Skip to content

Commit

Permalink
Update
Browse files Browse the repository at this point in the history
  • Loading branch information
aadi-bh committed Nov 28, 2023
1 parent a5e7269 commit 7f18900
Show file tree
Hide file tree
Showing 5 changed files with 38 additions and 19 deletions.
2 changes: 1 addition & 1 deletion .gitignore
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
**.DS_Store
**/.DS_Store
**__pycache__
Binary file removed _site/plots/.DS_Store
Binary file not shown.
41 changes: 23 additions & 18 deletions main.py
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,12 @@ def run(args):
tf = max(0, args.Tf)
cfl = min(1, args.cfl)

prefix = f"{args.pde}-{tf}-{args.ic}-{args.filter}-g{args.ggb}"
prefix = f"{args.pde}-{tf}-{args.ic}-{args.filter}-g"
if args.ggb:
prefix += str(args.ggblambda)
else:
prefix += "False"

print("PDE :", args.pde)
print("TF :", tf)
print("CFL :", cfl)
Expand All @@ -75,7 +80,6 @@ def run(args):

times = [0]
us = [u_hat_init]
print(tf)
if (tf > 0):
output = solve_ivp(fun = rhs,
t_span = [0, tf],
Expand All @@ -95,8 +99,8 @@ def run(args):
# Conveniently for us, the shock stays put at pi
left = np.where(x < pi, True, False)
right = np.where(x > pi, True, False)
ggbleft = ggb.recon_fft(x[left], uf, args.L)
ggbright = ggb.recon_fft(x[right], uf, args.L)
ggbleft = ggb.recon_fft(x[left], uf, args.ggblambda)
ggbright = ggb.recon_fft(x[right], uf, args.ggblambda)
v[left] = ggbleft
v[right] = ggbright
sols.append((times[-1], us[0], uf, v))
Expand All @@ -118,53 +122,54 @@ def run(args):
parser.add_argument('--Tf', type=float, default=1.0, help = "Final time")
parser.add_argument('--pde', choices=('linadv', 'burgers'), default='burgers', help = "PDE to solve")
# parser.add_argument('--add_visc', type=bool, default=False)
parser.add_argument('--filter', choices=('no_filter', 'exponential', 'cesaro', 'raisedcos', 'lanczos'),
default='no_filter', help = "Which filter, if any")
parser.add_argument('--filter', choices=('no_filter', 'exponential', 'cesaro', 'raisedcos', 'lanczos'), default='no_filter', help = "Which filter, if any")
# parser.add_argument('--filterp', type=int, default=1, "p value for the exponential")
parser.add_argument('--ggb', type=bool, default=False, action='store_true'
parser.add_argument('--ggb', action='store_true', default=False,
help = "Whether to reconstruct the analytic part of the solution")
parser.add_argument('-L', type=int, default=3, help = "Number of elements in the GGB basis")
parser.add_argument('---ggblambda', type=int, default=3, help = "Number of elements in the GGB basis")
# parser.add_argument('--max_lgN', type=int, default=7, help = "Largest power of 2 to calculate until")
#parser.add_argument('--integrator', choices=('solve_ivp', 'elrk4'), default='elrk4')
parser.add_argument('--show_markers', type=bool, default=False, action='store_true'
parser.add_argument('--show_markers',default=False, action='store_true',
help = "Show the original sample values in red crosses")
parser.add_argument('--exact', type=str, default=None, help = "Exact solution file")
args = parser.parse_args()
if args.N == None:
args.N = [16, 64]
args.N = [16, 64, 128]
sols, prefix, initial_condition = run(args)
plotname = prefix+'.svg'

# PLOT!
nn = 2048
fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 8), width_ratios=[3,1])
# Exact
if args.exact != None:
exact = np.loadtxt(args.exact).transpose()
exact[0] *= 2 * pi
ax[0].plot(exact[0], exact[1], 'ko', markersize=0.8, label="Exact")
# Initial
x = np.linspace(xmin, xmax, 1000)
ax[0].plot(x, initial_condition(x), linewidth=1, color='k', label="Init")
# All the modes
for (t, uhinit, uf, v) in sols:
n = len(uf)
label = str(n) + f", t={np.round(t, 3)}"
plots.smoothplot(uf, ax[0], label=label, linewidth=1)
plots.plot_resolution(uf, ax[1], linewidth=0.5, markersize=0.5)
# plots.smoothplot(uf, ax[0], label=label)
if args.ggb:
label += ",ggb"
if args.filter != 'no_filter':
label += "," + args.filter
ax[0].plot(cgrid(n), v, label=label)
elif args.filter != 'no_filter':
label += "," + args.filter
plots.smoothplot(uf, ax[0])
plots.smoothplot(uf, ax[0], label=label)
else:
# So no filter, no ggb
plots.smoothplot(uf, ax[0], label=label)
if args.show_markers:
ax[0].plot(cgrid(n), ifft(uf), "+", color= "red", markersize=5)
# Exact
if args.exact != None:
exact = np.loadtxt(args.exact).transpose()
ax[0].plot(exact[0] * 2 * np.pi, exact[1], 'ko', markersize=0.8, label="Exact")

ax[0].grid(visible=True)
ax[0].legend()
fig.tight_layout()
fig.savefig(plotname)
plt.close()
print("Saved plot to "+plotname)
13 changes: 13 additions & 0 deletions plots.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,19 @@ def plot_resolution(c, ax, **kwargs):
ax.semilogy(k, np.abs(fftshift(c)), **kwargs)
return

def plot_error_fft(uh, exact, ax, **kwargs):
x = exact[0]
if len(uh) < len(x):
u = ifft_at(x, uh)
e = np.abs(u - exact[1])
ax[1].plot(x, e, **kwargs)
else:
print("WARNING: Interpolating exact not supported yet.\nSupply exact solution on finer grid.")
return

def plot_error(c, exact, ax, **kwargs):
print("print_error not implemented yet.")

# Fourier interpolate the IFFT
def smoothplot(v, ax,nn=2048, **plotargs):
n = len(v)
Expand Down
1 change: 1 addition & 0 deletions presentation.mplstyle
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
axes.titlesize : 24
axes.labelsize : 20
axes.linewidth : 1.5
lines.linewidth : 3
lines.markersize : 10
xtick.labelsize : 16
Expand Down

0 comments on commit 7f18900

Please sign in to comment.