diff --git a/advection/advection.py b/advection/advection.py index 53a251f..5a69a60 100644 --- a/advection/advection.py +++ b/advection/advection.py @@ -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): diff --git a/advection/weno.py b/advection/weno.py index c73550b..e5a1a81 100644 --- a/advection/weno.py +++ b/advection/weno.py @@ -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 @@ -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) @@ -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 = [] @@ -458,18 +378,12 @@ 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) @@ -477,14 +391,9 @@ def rk_substep_scipy(t, y): 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) @@ -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() \ No newline at end of file