Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WENO convergence #2

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

WENO convergence #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
050390f
Try using 2-norm: doesn't help?
IanHawke Feb 27, 2019
e6d832a
Fix WENO convergence plot
IanHawke Feb 27, 2019
afb7f9f
Switch RK4 case back to Gaussian
IanHawke Feb 27, 2019
a480adb
Do convergence on 5 periods
IanHawke Mar 5, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 7 additions & 4 deletions advection/advection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,14 +86,17 @@ def fill_BCs(self):
# right boundary
self.a[self.ihi+1+n] = self.a[self.ilo+n]

def norm(self, e):
def norm(self, e, norm="inf"):
""" return the norm of quantity e which lives on the grid """
if len(e) != 2*self.ng + self.nx:
return None

#return np.sqrt(self.dx*np.sum(e[self.ilo:self.ihi+1]**2))
return np.max(abs(e[self.ilo:self.ihi+1]))

if norm==2:
return np.sqrt(self.dx*np.sum(e[self.ilo:self.ihi+1]**2))
elif norm=="inf":
return np.max(abs(e[self.ilo:self.ihi+1]))
else:
raise NotImplementedError("Only norms implemented are '2' and 'inf'")

class Simulation(object):

Expand Down
140 changes: 25 additions & 115 deletions advection/weno.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -287,95 +287,16 @@ def rk_substep_scipy(t, y):
label="WENO3")


# #-------------------------------------------------------------------------
# # convergence test
# # Note that WENO schemes with standard weights lose convergence at
# # critical points. For high degree critical points they lose more orders.
# # The suggestion in Gerolymos is that you may expect to drop down to
# # order r-1 in the limit.
# # The Gaussian has all odd derivatives vanishing at the origin, so
# # the higher order schemes will lose accuracy.
# # For the Gaussian:
# # This shows clean 5th order convergence for r=3
# # But for r=4-6 the best you get is ~6th order, and 5th order is more
# # realistic
# # For sin(x - sin(x)) type data Gerolymos expects better results
# # But the problem actually appears to be the time integrator
# # Switching to Dormand-Price 8th order from scipy (a hack) will make it
# # work for all cases. With sin(.. sin) data you get 2r - 2 thanks to
# # the one critical point.
#
# problem = "sine_sine"
#
# xmin =-1.0
# xmax = 1.0
## orders = [4]
# orders = [3, 4, 5, 6]
## N1 = [2**4*3**i//2**i for i in range(5)]
## N2 = [2**5*3**i//2**i for i in range(6)]
## N3 = [3**4*4**i//3**i for i in range(5)]
## N4 = [2**(4+i) for i in range(4)]
## N = numpy.unique(numpy.array(N1+N2+N3+N4, dtype=numpy.int))
## N.sort()
## N = [32, 64, 128, 256, 512]
## N = [32, 64, 128]
# N = [24, 32, 54, 64, 81, 108, 128]
#
# errs = []
# errsM = []
#
# u = 1.0
#
# colors="bygrc"
#
# for order in orders:
# ng = order+1
# errs.append([])
# errsM.append([])
# for nx in N:
# print(order, nx)
# gu = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax)
# su = WENOSimulation(gu, u, C=0.5, weno_order=order)
## guM = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax)
## suM = WENOMSimulation(guM, u, C=0.5, weno_order=order)
#
# su.init_cond("sine_sine")
## suM.init_cond("sine_sine")
# ainit = su.grid.a.copy()
#
# su.evolve_scipy(num_periods=1)
## suM.evolve_scipy(num_periods=1)
#
# errs[-1].append(gu.norm(gu.a - ainit))
## errsM[-1].append(guM.norm(guM.a - ainit))
#
# pyplot.clf()
# N = numpy.array(N, dtype=numpy.float64)
# for n_order, order in enumerate(orders):
# pyplot.scatter(N, errs[n_order],
# color=colors[n_order],
# label=r"WENO, $r={}$".format(order))
## pyplot.scatter(N, errsM[n_order],
## color=colors[n_order],
## label=r"WENOM, $r={}$".format(order))
# pyplot.plot(N, errs[n_order][0]*(N[0]/N)**(2*order-2),
# linestyle="--", color=colors[n_order],
# label=r"$\mathcal{{O}}(\Delta x^{{{}}})$".format(2*order-2))
## pyplot.plot(N, errs[n_order][len(N)-1]*(N[len(N)-1]/N)**4,
## color="k", label=r"$\mathcal{O}(\Delta x^4)$")
#
# ax = pyplot.gca()
# ax.set_ylim(numpy.min(errs)/5, numpy.max(errs)*5)
# ax.set_xscale('log')
# ax.set_yscale('log')
#
# pyplot.xlabel("N")
# pyplot.ylabel(r"$\| a^\mathrm{final} - a^\mathrm{init} \|_2$",
# fontsize=16)
#
# pyplot.legend(frameon=False)
# pyplot.savefig("weno-converge-sine-sine.pdf")
## pyplot.show()
#-------------------------------------------------------------------------
# convergence test
# Note that WENO schemes with standard weights lose convergence at
# critical points. For high degree critical points they lose more orders.
# The suggestion in Gerolymos is that you may expect to drop down to
# order r-1 in the limit.
#
# For the odd r values and using sine initial data we can get optimal
# convergence using 8th order time integration. For other cases the
# results are not so nice.

#-------------- RK4

Expand Down Expand Up @@ -405,7 +326,7 @@ def rk_substep_scipy(t, y):

su.evolve(num_periods=5)

errs[-1].append(gu.norm(gu.a - ainit))
errs[-1].append(gu.norm(gu.a - ainit, norm=2))

pyplot.clf()
N = numpy.array(N, dtype=numpy.float64)
Expand All @@ -431,17 +352,16 @@ def rk_substep_scipy(t, y):

pyplot.legend(frameon=False)
pyplot.savefig("weno-converge-gaussian-rk4.pdf")
# pyplot.show()
pyplot.show()

#-------------- Gaussian
#-------------- Sine wave, 8th order time integrator

problem = "gaussian"
problem = "sine"

xmin = 0.0
xmax = 1.0
orders = [3, 4, 5, 6]
orders = [3, 5, 7]
N = [24, 32, 54, 64, 81, 108, 128]
# N = [32, 64, 108, 128]

errs = []
errsM = []
Expand All @@ -458,33 +378,22 @@ def rk_substep_scipy(t, y):
print(order, nx)
gu = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax)
su = WENOSimulation(gu, u, C=0.5, weno_order=order)
# guM = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax)
# suM = WENOMSimulation(guM, u, C=0.5, weno_order=order)

su.init_cond("gaussian")
# suM.init_cond("gaussian")
su.init_cond("sine")
ainit = su.grid.a.copy()

su.evolve_scipy(num_periods=1)
# suM.evolve_scipy(num_periods=1)

errs[-1].append(gu.norm(gu.a - ainit))
# errsM[-1].append(guM.norm(guM.a - ainit))
su.evolve_scipy(num_periods=5)
errs[-1].append(gu.norm(gu.a - ainit, norm=2))

pyplot.clf()
N = numpy.array(N, dtype=numpy.float64)
for n_order, order in enumerate(orders):
pyplot.scatter(N, errs[n_order],
color=colors[n_order],
label=r"WENO, $r={}$".format(order))
# pyplot.scatter(N, errsM[n_order],
# color=colors[n_order],
# label=r"WENOM, $r={}$".format(order))
pyplot.plot(N, errs[n_order][0]*(N[0]/N)**(2*order-2),
pyplot.plot(N, errs[n_order][0]*(N[0]/N)**(2*order-1),
linestyle="--", color=colors[n_order],
label=r"$\mathcal{{O}}(\Delta x^{{{}}})$".format(2*order-2))
# pyplot.plot(N, errs[n_order][len(N)-1]*(N[len(N)-1]/N)**4,
# color="k", label=r"$\mathcal{O}(\Delta x^4)$")
label=r"$\mathcal{{O}}(\Delta x^{{{}}})$".format(2*order-1))

ax = pyplot.gca()
ax.set_ylim(numpy.min(errs)/5, numpy.max(errs)*5)
Expand All @@ -494,9 +403,10 @@ def rk_substep_scipy(t, y):
pyplot.xlabel("N")
pyplot.ylabel(r"$\| a^\mathrm{final} - a^\mathrm{init} \|_2$",
fontsize=16)
pyplot.title("Convergence of Gaussian, DOPRK8")
pyplot.title("Convergence of sine wave, DOPRK8")

pyplot.legend(frameon=False)
pyplot.savefig("weno-converge-gaussian.pdf")
# pyplot.show()
lgd = ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
pyplot.savefig("weno-converge-sine.pdf",
bbox_extra_artists=(lgd,), bbox_inches='tight')
pyplot.show()