clean up advection example and doing it in 1d

examples/advection/weak.py

@@ -13,106 +13,127 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-
 import numpy as np
-import pyopencl as cl  # noqa
-import pyopencl.array  # noqa
-import pyopencl.clmath  # noqa
+import numpy.linalg as la
 
-import pytest  # noqa
+import pyopencl as cl
+import pyopencl.array
+import pyopencl.clmath
 
-from pyopencl.tools import (  # noqa
-                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+from grudge import bind, sym
 
 import logging
+
 logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.INFO)
 
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-
-import numpy.linalg as la
 
+def main(ctx_factory, dim=1, order=4, visualize=True):
+    cl_ctx = ctx_factory()
+    queue = cl.CommandQueue(cl_ctx)
 
+    # {{{ parameters
 
+    # domain side [-d/2, d/2]^dim
+    d = 1.0
+    # number of points in each dimension
+    npoints = 20
+    # grid spacing
+    h = d / npoints
+    # cfl?
+    dt_factor = 4
+    # final time
+    final_time = 1.0
 
-def main(write_output=True, order=4):
-    cl_ctx = cl.create_some_context()
-    queue = cl.CommandQueue(cl_ctx)
+    # velocity field
+    c = np.array([0.1] * dim)
+    norm_c = la.norm(c)
+    # flux
+    flux_type = "central"
 
-    dim = 2
+    # compute number of steps
+    dt = dt_factor * h / order**2
+    nsteps = int(final_time // dt) + 1
+    dt = final_time/nsteps + 1.0e-15
 
+    # }}}
 
-    from meshmode.mesh.generation import generate_regular_rect_mesh
-    mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
-            n=(20, 20), order=order)
+    # {{{ discretization
 
-    dt_factor = 4
-    h = 1/20
+    from meshmode.mesh.generation import generate_box_mesh
+    mesh = generate_box_mesh(
+            [np.linspace(-d/2, d/2, npoints) for _ in range(dim)],
+            order=order)
 
+    from grudge import DGDiscretizationWithBoundaries
     discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
 
-    c = np.array([0.1,0.1])
-    norm_c = la.norm(c)
-
+    # }}}
 
-    flux_type = "central"
-         
+    # {{{ solve advection
 
     def f(x):
-        return sym.sin(3*x)
+        return sym.sin(3 * x)
 
     def u_analytic(x):
-        return f(-np.dot(c, x)/norm_c+sym.var("t", sym.DD_SCALAR)*norm_c)
+        t = sym.var("t", sym.DD_SCALAR)
+        return f(-np.dot(c, x) / norm_c + t * norm_c)
 
     from grudge.models.advection import WeakAdvectionOperator
-    from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
     op = WeakAdvectionOperator(c,
         inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
         flux_type=flux_type)
 
     bound_op = bind(discr, op.sym_operator())
-
     u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0)
 
+    return
+
     def rhs(t, u):
         return bound_op(queue, t=t, u=u)
 
-    final_time = 0.3
-
-    dt = dt_factor * h/order**2
-    nsteps = (final_time // dt) + 1
-    dt = final_time/nsteps + 1e-15
-
-
     from grudge.shortcuts import set_up_rk4
     dt_stepper = set_up_rk4("u", dt, u, rhs)
 
-    last_u = None
-
-    from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
+    if dim == 1:
+        import matplotlib.pyplot as pt
+        pt.figure(figsize=(8, 8), dpi=300)
+    else:
+        from grudge.shortcuts import make_visualizer
+        vis = make_visualizer(discr, vis_order=order)
 
     step = 0
-
     norm = bind(discr, sym.norm(2, sym.var("u")))
-
     for event in dt_stepper.run(t_end=final_time):
-        if isinstance(event, dt_stepper.StateComputed):
-
-            step += 1
-
-            #print(step, event.t, norm(queue, u=event.state_component[0]))
-            vis.write_vtk_file("fld-weak-%04d.vtu" % step,
-                    [  ("u", event.state_component) ])
-
-
+        if not isinstance(event, dt_stepper.StateComputed):
+            continue
 
+        step += 1
+        norm_u = norm(queue, u=event.state_component)
+        logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
 
+        if dim == 1:
+            x = discr.discr_from_dd(sym.DTAG_VOLUME_ALL).nodes().get(queue)
+            u = event.state_component.get(queue)
 
+            pt.plot(x, u)
+            pt.xlabel("$x$")
+            pt.ylabel("$u$")
+            pt.title("t = {:.2f}".format(event.t))
+            pt.savefig("fld-weak-%04d.png" % step)
+            pt.clf()
+        else:
+            vis.write_vtk_file("fld-weak-%04d.vtu" % step, [
+                ("u", event.state_component)
+                ], overwrite=True)
 
 
 if __name__ == "__main__":
-    main()
+    import argparse
 
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--dim", default=2, type=int)
+    args = parser.parse_args()
 
+    main(cl.create_some_context,
+            dim=args.dim)

grudge/models/advection.py

@@ -74,7 +74,7 @@ class AdvectionOperatorBase(HyperbolicOperator):
 
 class StrongAdvectionOperator(AdvectionOperatorBase):
     def flux(self, u):
-        normal = sym.normal(u. dd, self.ambient_dim)
+        normal = sym.normal(u.dd, self.ambient_dim)
         v_dot_normal = sym.cse(self.v.dot(normal), "v_dot_normal")
 
         return u.int * v_dot_normal - self.weak_flux(u)
@@ -142,7 +142,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
 
     def flux(self, u):
 
-        normal = sym.normal(u. dd, self.ambient_dim)
+        normal = sym.normal(u.dd, self.ambient_dim)
 
         surf_v = sym.interp("vol", u.dd)(self.v)
 

grudge/symbolic/compiler.py

@@ -1047,6 +1047,11 @@ class ToLoopyInstructionMapper(object):
 
         self.insn_count += 1
 
+        for expr in insn.exprs:
+            print(expr)
+            print(self.dd_inference_mapper(expr))
+        print()
+
         from pytools import single_valued
         governing_dd = single_valued(
                 self.dd_inference_mapper(expr)