Skip to content

Commit

Permalink
More bugs
Browse files Browse the repository at this point in the history
  • Loading branch information
@aadi-bh
aadi-bh committed Nov 24, 2023
1 parent 148496d commit 5cfe38d
Show file tree
Hide file tree
Showing 2 changed files with 103 additions and 96 deletions.
197 changes: 102 additions & 95 deletions burgers_filtered.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,13 +18,15 @@
xmax = 2 * pi

def saw_tooth(x):
left = np.where(x <= pi, 1, 0);
return x * left + (x - 2 * pi) * (1-left)
left = np.where(x < pi, 1, 0);
return x * left + (x - 2 * pi) * (1-left) / (2*pi)
def sin(x):
return np.sin(x)
def step(x):
y = np.where(x < np.pi, 1, 0) * np.where(x > np.pi/2, 1, 0)
y = np.where(x < pi, 1, 0) * np.where(x > pi/2, 1, 0)
return y
def bump(x):
return np.exp(-6 * (x-pi)**2)

def plot_resolution(c, ax, kwargs):
k = fftshift(freqs(len(c)))
Expand All @@ -38,15 +40,15 @@ def pad(c, m):
r = fftshift(c)
r = np.r_[np.zeros(m), r, np.zeros(m)]
r *= 1 # N / newN
r = fftshift(r)
r = ifftshift(r)
return r
def unpad(c, m):
N = len(c)
newN = N - 2 * m
r = fftshift(c)
r = r[m:m + newN]
r *= 1 # N / newN
r = fftshift(r)
r = ifftshift(r)
return r
def freqs(n, x1=xmin, x2=xmax):
return fftfreq(n, 1./ (n))
Expand All @@ -55,6 +57,7 @@ def mask(c, tol=1e-14):
return c * np.where(np.abs(c)<tol, 0, 1)

def sigma(eta, p=1):
return 1 - eta
return np.exp(-15* np.power(eta, 2*p))

def create_filter(k, sigma, args):
Expand All @@ -67,25 +70,26 @@ def no_filter(k, sigma, args):
def zero(t, u_hat, N, M, eps, filter, a=1):
return np.zeros(len(u_hat));

def linadv(t, u_hat, N, M, filter, a=1):
def linadv(t, u_hat, N, M, filter, a = 2 * pi):
# because fft scales these, and we want the ifft to work as usual
NN = (2 * M//2) + N
u_hat = pad(u_hat, M//2) * NN / N
kk = freqs(NN)
nonlinear = -1j * kk * u_hat
nonlinear *= a
nonlinear *= filter
return unpad(mask(nonlinear), M//2)
return unpad(mask(nonlinear), M//2) * N / NN

def burgers(t, u_hat, N, M, filter, a=1):
# because fft scales these, and we want the ifft to work as usual
NN = (2 * M//2) + N
u_hat = pad(u_hat, M//2) * NN / N
u = ifft(u_hat)
kk = freqs(NN)
nonlinear = -1 * fft(0.5 * u * u) * 1j * kk
nonlinear = -1j * kk * fft(0.5 * u * u)
nonlinear *= a
nonlinear *= filter
return unpad(mask(nonlinear), M//2)
return unpad(mask(nonlinear), M//2) * N / NN

def semigroup_heat(dt, k, eps):
S_half = sparse.diags(np.exp(-1 * eps * k**2 * dt / 2))
Expand All @@ -103,7 +107,7 @@ def semigroup_none(dt, k, eps):
# choices=('C','LF','GLF','LLF','LW','ROE','EROE','GOD'),
# help='Scheme', default='LF')
parser.add_argument('--ic',
choices=('saw_tooth','sin','step'),
choices=('saw_tooth','sin','step', 'bump'),
help='Initial condition', default='sin')
parser.add_argument('--Tf', type=float, help='Final time', default=1.0)
parser.add_argument('--pde', choices=('linadv', 'burgers'), default='burgers')
Expand All @@ -113,89 +117,92 @@ def semigroup_none(dt, k, eps):
parser.add_argument('--integrator', choices=('solve_ivp', 'elrk4'), default='elrk4')
args = parser.parse_args()

a = 2 * pi
rhs = linadv
initial_condition = sin
visc = semigroup_none
USE_FILTER = False
tf = 2
cfl = 0.5
if args.pde == 'burgers':
a = 1
rhs = burgers
if args.ic == 'saw_tooth':
initial_condition = saw_tooth
elif args.ic == 'sin':
if __name__ == '__main__':
rhs = linadv
initial_condition = sin
elif args.ic == 'step':
initial_condition = step
if args.add_visc == True:
visc = semigroup_heat
tf = max(0, args.Tf)
cfl = min(1, args.cfl)
USE_FILTER = args.use_filter

dt = 0.02
fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 10), width_ratios=[3,1])
fig.tight_layout()
x = np.linspace(xmin, xmax, 1000)
ax[0].plot(x, initial_condition(x), linewidth=0.1, color='k', label="init")

for i in range(0, args.max_lgN - 4 + 1):
N = np.power(2, i + 4);
print(f"N={N}")
M = 10 * N // 2;
m = M // 2;
NN = (2 * m) + N
dx = (xmax-xmin) / N
dt = min(dt, cfl * dx)
x = np.arange(xmin, xmax, dx)
k = freqs(N)
kk = freqs(NN)
filter = np.ones(len(kk))
p = i
if USE_FILTER == True:
filter = create_filter(kk, sigma, {'p':p})
args = (N, M, filter, a)
u_hat_init = fft(initial_condition(x))
S_half, S = visc(dt, k, eps = 1e-2)
'''
output = solve_ivp(fun = rhs,
t_span = [0, tf],
t_eval = [0, tf],
y0 = u_hat_init,
args = args)
times, u = output.t, output.y.transpose()
'''
times, u = elrk4([S_half, S], rhs, u_hat_init, (0, tf), dt, args)
ax[0].plot(x, ifft(u[-1]).real, label=str(N)+f", t={np.round(times[-1], 3)}")
plot_resolution(u[-1], ax[1], {'linewidth':0.1, 'markersize':0.1})
'''
nr = 11
f, a = plt.subplots(ncols=2, nrows = nr,
figsize=(10, 2 * nr),width_ratios=[3,1])
f.tight_layout()
for k in range(nr):
j = int(np.linspace(0, len(times), nr, endpoint=True)[k])
print(j)
a[k][0].plot(x, ifft(u[j]).real, label=str(np.round(times[j], 3)))
a[k][0].legend()
plot_resolution(u[j], a[k][1], {'color':'k', 'linewidth':0.1, 'markersize':0.1})
f.savefig(f"{N}-f{USE_FILTER}-f{str(p)}.png")
'''
np.savetxt(f"{N}.txt", np.vstack((k, u[-1].real)))
print("Saved file.")

if rhs == burgers:
if tf == 1:
god = np.loadtxt("burg3_GOD_5.txt").transpose()
if tf == 0.6:
god = np.loadtxt("burg3_GOD_3.txt").transpose()
if tf == 0.2:
god = np.loadtxt("burg3_GOD_1.txt").transpose()
ax[0].plot(god[0] * 2 * np.pi, god[1], 'ko', markersize=0.1, label="Godunov flux")
ax[0].legend()
fig.savefig(f"f{rhs}-f{str(USE_FILTER)}.png")

plt.close()
print("Done.")
visc = semigroup_none
USE_FILTER = False
tf = 2
cfl = 0.5
if args.pde == 'burgers':
rhs = burgers
if args.ic == 'saw_tooth':
initial_condition = saw_tooth
elif args.ic == 'sin':
initial_condition = sin
elif args.ic == 'step':
initial_condition = step
elif args.ic == 'bump':
initial_condition = bump
if args.add_visc == True:
visc = semigroup_heat
tf = max(0, args.Tf)
cfl = min(1, args.cfl)
USE_FILTER = args.use_filter

dt = 0.02
fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 8), width_ratios=[3,1])
fig.tight_layout()
x = np.linspace(xmin, xmax, 1000)
ax[0].plot(x, initial_condition(x), linewidth=0.1, color='k', label="init")

sols = []
for i in range(0, args.max_lgN - 4 + 1):
N = np.power(2, i + 4);
print(f"N={N}")
M = 3 * N // 2;
m = M // 2;
NN = (2 * m) + N
dx = (xmax-xmin) / N
dt = min(dt, cfl * dx)
x = np.arange(xmin, xmax, dx)
k = freqs(N)
kk = freqs(NN)
filter = np.ones(len(kk))
p = i
if USE_FILTER == True:
filter = create_filter(kk, sigma, {'p':p})
args = (N, M, filter)
u_hat_init = fft(initial_condition(x))
S_half, S = visc(dt, k, eps = 1e-2)
'''
output = solve_ivp(fun = rhs,
t_span = [0, tf],
t_eval = [0, tf],
y0 = u_hat_init,
args = args)
times, u = output.t, output.y.transpose()
'''
times, u = elrk4([S_half, S], rhs, u_hat_init, (0, tf), dt, args)
sols.append((times, u))
ax[0].plot(x, ifft(u[-1]).real, label=str(N)+f", t={np.round(times[-1], 3)}")
plot_resolution(u[-1], ax[1], {'linewidth':0.5, 'markersize':0.5})
'''
nr = 11
f, a = plt.subplots(ncols=2, nrows = nr,
figsize=(10, 2 * nr),width_ratios=[3,1])
f.tight_layout()
for k in range(nr):
j = int(np.linspace(0, len(times), nr, endpoint=True)[k])
print(j)
a[k][0].plot(x, ifft(u[j]).real, label=str(np.round(times[j], 3)))
a[k][0].legend()
plot_resolution(u[j], a[k][1], {'color':'k', 'linewidth':0.1, 'markersize':0.1})
f.savefig(f"{N}-{USE_FILTER}-{str(p)}.png")
'''
np.savetxt(f"{N}.txt", np.vstack((k, u[-1].real)))
print("Saved file.")

if rhs == burgers and initial_condition == sin:
if tf >= 1 or True:
god = np.loadtxt("burg3_GOD_5.txt").transpose()
elif tf == 0.6:
god = np.loadtxt("burg3_GOD_3.txt").transpose()
elif tf == 0.2:
god = np.loadtxt("burg3_GOD_1.txt").transpose()
ax[0].plot(god[0] * 2 * np.pi, god[1], 'ko', markersize=0.1, label="Godunov flux")
ax[0].legend()
fig.savefig(f"{rhs}-{str(USE_FILTER)}.png")

plt.close()
print("Done.")
2 changes: 1 addition & 1 deletion common.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ def elrk4(SemiGroup,Nonlinear,y0,tinterval,dt,args):
flag = False
elif np.any(np.isnan(y)):
flag = False
print("t, dt = ", t, dt)
# print("t, dt = ", t, dt)

times = np.array(time)
solution.append(y)
Expand Down

0 comments on commit 5cfe38d

Please sign in to comment.