add norm and restrict methods and a convergence plot

Author: zingale

multiphysics/burgersvisc.py

@@ -43,6 +43,10 @@ def scratch_array(self):
         return np.zeros((2*self.ng+self.nx), dtype=np.float64)
 
 
+    def norm(self, e):
+        return np.sqrt(self.dx*np.sum(e[self.ilo:self.ihi+1]**2))
+
+        
     def fill_BCs(self, var):
 
         if not var in self.data.keys():
@@ -54,6 +58,21 @@ def fill_BCs(self, var):
         vp[0:self.ilo+1] = vp[self.ilo]
         vp[self.ihi+1:] = vp[self.ihi]
 
+    def restrict(self, fac=2):
+        """ restrict the data q that lives on this grid by a 
+        factor fac and return the new data"""
+        
+        # we require fac to be a multiple of 2
+
+        gnew = Grid1d(self.nx//fac, ng=self.ng, 
+                      xmin=self.xmin, xmax=self.xmax, vars=self.data.keys())
+
+        for v in self.data:
+            b = self.data[v][self.ilo:self.ihi+1]
+            gnew.data[v][gnew.ilo:gnew.ihi+1] = \
+                        np.mean(b.reshape(-1, fac), axis=1)
+        return gnew
+
 
 class Simulation(object):
 
@@ -120,7 +139,7 @@ def diffuse(self, eps, S, dt):
 
         return unew
 
-    def advect(self, S, dt):
+    def advect(self, S, dt, limit=1):
         """ compute the advective term that updates u in time.  Here, S is
         a source term """
 
@@ -143,8 +162,11 @@ def advect(self, S, dt):
         dl[ib:ie+1] = np.fabs(u[ib+1:ie+2] - u[ib  :ie+1])
         dr[ib:ie+1] = np.fabs(u[ib  :ie+1] - u[ib-1:ie ])
 
-        minslope = np.minimum(np.fabs(dc), np.minimum(2.0*dl, 2.0*dr))
-        ldeltau = np.where(test > 0.0, minslope, 0.0)*np.sign(dc)
+        if limit:
+            minslope = np.minimum(np.fabs(dc), np.minimum(2.0*dl, 2.0*dr))
+            ldeltau = np.where(test > 0.0, minslope, 0.0)*np.sign(dc)
+        else:
+            ldeltau = dc
 
         # construct the interface states, to second order in space and
         # time
@@ -193,7 +215,7 @@ def lap(self):
 
         return lapu
 
-    def evolve(self, eps, cfl, tmax, dovis=0):
+    def evolve(self, eps, cfl, tmax, limit=1, dovis=0):
         """ 
         the main evolution loop.  Evolve 
   
@@ -225,7 +247,7 @@ def evolve(self, eps, cfl, tmax, dovis=0):
             S = eps*self.lap()
 
             # construct the advective update
-            A = self.advect(S, dt)
+            A = self.advect(S, dt, limit=limit)
 
             # diffuse for dt
             unew = self.diffuse(eps, A, dt)
@@ -247,43 +269,43 @@ def evolve(self, eps, cfl, tmax, dovis=0):
                 plt.draw()
 
 
-nx = 256
-cfl = 0.8
-tmax = 0.2
-
-g1 = Grid1d(nx, ng=2, vars=["u"])
-g2 = Grid1d(nx, ng=2, vars=["u"])
-g3 = Grid1d(nx, ng=2, vars=["u"])
 
-eps1 = 0.005
-s1 = Simulation(g1)
-s1.init_cond()
-s1.evolve(eps1, cfl, tmax, dovis=0)
+if __name__ == "__main__":
+    nx = 256
+    cfl = 0.8
+    tmax = 0.2
 
+    g1 = Grid1d(nx, ng=2, vars=["u"])
+    g2 = Grid1d(nx, ng=2, vars=["u"])
+    g3 = Grid1d(nx, ng=2, vars=["u"])
 
-eps2 = 0.0005
-s2 = Simulation(g2)
-s2.init_cond()
-s2.evolve(eps2, cfl, tmax, dovis=0)
+    eps1 = 0.005
+    s1 = Simulation(g1)
+    s1.init_cond()
+    s1.evolve(eps1, cfl, tmax, dovis=0)
 
+    eps2 = 0.0005
+    s2 = Simulation(g2)
+    s2.init_cond()
+    s2.evolve(eps2, cfl, tmax, dovis=0)
 
-eps3 = 0.00005
-s3 = Simulation(g3)
-s3.init_cond()
-s3.evolve(eps3, cfl, tmax, dovis=0)
+    eps3 = 0.00005
+    s3 = Simulation(g3)
+    s3.init_cond()
+    s3.evolve(eps3, cfl, tmax, dovis=0)
 
-u1 = s1.grid.data["u"]
-plt.plot(g1.x, u1, label=r"$\epsilon = %f$" % (eps1))
+    u1 = s1.grid.data["u"]
+    plt.plot(g1.x, u1, label=r"$\epsilon = %f$" % (eps1))
 
-u2 = s2.grid.data["u"]
-plt.plot(g2.x, u2, label=r"$\epsilon = %f$" % (eps2), ls="--", color="k")
+    u2 = s2.grid.data["u"]
+    plt.plot(g2.x, u2, label=r"$\epsilon = %f$" % (eps2), ls="--", color="k")
 
-u3 = s3.grid.data["u"]
-plt.plot(g3.x, u3, label=r"$\epsilon = %f$" % (eps3), ls=":", color="k")
+    u3 = s3.grid.data["u"]
+    plt.plot(g3.x, u3, label=r"$\epsilon = %f$" % (eps3), ls=":", color="k")
 
-plt.legend(frameon=False)
+    plt.legend(frameon=False)
 
-plt.xlim(0.0, 1.0)
+    plt.xlim(0.0, 1.0)
 
-plt.tight_layout()
-plt.savefig("burgersvisc.pdf")
+    plt.tight_layout()
+    plt.savefig("burgersvisc.pdf")

multiphysics/burgersvisc_converge.py

@@ -0,0 +1,54 @@
+# test convergence of our inviscid Burgers problem by comparing against
+# a high-resolution test case
+
+import burgersvisc as bv
+import matplotlib.pyplot as plt
+import numpy as np
+
+cfl = 0.8
+tmax = 0.2
+eps = 0.005
+
+# high-resolution benchmark
+nxb = 4096
+gb = bv.Grid1d(nxb, ng=2, vars=["u"])
+
+sb = bv.Simulation(gb)
+sb.init_cond()
+sb.evolve(eps, cfl, tmax, dovis=0)
+
+
+# comparison simulations
+N = [64, 128, 256, 512, 1024]
+
+err = []
+for nx in N:
+    # our comparison -- a coarsen version of the high-resolution
+    # solution
+    gc = gb.restrict(fac=nxb//nx)
+
+    g = bv.Grid1d(nx, ng=2, vars=["u"])
+    s = bv.Simulation(g)
+    s.init_cond()
+    s.evolve(eps, cfl, tmax, dovis=0)
+
+    uc = gc.data["u"]
+    u = g.data["u"]
+
+    e = uc - u
+    err.append(g.norm(e))
+
+N = np.array(N, dtype=np.float64)
+err = np.array(err)
+
+plt.scatter(N, err, color="C1", label="Burgers' solution")
+plt.loglog(N, err[len(N)-1]*(N[len(N)-1]/N)**2, color="C0", 
+           label="$\mathcal{O}(\Delta x^2)$")
+
+plt.xlabel(r"$N$", fontsize="large")
+plt.ylabel(r"L2 norm of absolute error")
+
+plt.legend(frameon=False)
+
+plt.tight_layout()
+plt.savefig("burgersvisc_converge.pdf")

multiphysics/diffusion-reaction.py

@@ -280,7 +280,7 @@ def evolve(nx, kappa, tau, tmax, dovis=0, return_initial=0):
 
 plt.xlabel("$x$")
 plt.ylabel("$\phi$")
-plt.title(r"Diffusion-Reaction, $N = {}, \, \kappa = {:3.2f}, \, \tau = {:3.2f}$".format(nx, kappa, tau, dt))
+#plt.title(r"Diffusion-Reaction, $N = {}, \, \kappa = {:3.2f}, \, \tau = {:3.2f}$".format(nx, kappa, tau, dt))
 #, \, \delta t = {:3.2f}$ (between lines)
 plt.tight_layout()