doc/misc.rst

@@ -17,7 +17,11 @@ MacOS support is in the works.
     Everywhere else, just making sure you have the ``g++`` package should be
     enough.
 
-#.  Install `miniforge <https://github.com/conda-forge/miniforge>`__.
+#.  Install `miniforge <https://github.com/conda-forge/miniforge>`_::
+
+        curl -L -O https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh
+        # then run
+        bash ./Miniforge3-*.sh
 
 #.  ``export CONDA=/WHERE/YOU/INSTALLED/miniforge3``
 
@@ -33,7 +37,7 @@ MacOS support is in the works.
 
 #.  Type the following command::
 
-        hash -r; for i in pymbolic cgen genpy modepy pyvisfile loopy meshmode dagrt leap grudge; do python -m pip install --editable "git+https://gitlab.tiker.net/inducer/$i.git#egg=$i"; done
+        hash -r; for i in pymbolic cgen genpy modepy pyvisfile loopy arraycontext meshmode dagrt leap grudge; do python -m pip install --editable "git+https://github.com/inducer/$i.git#egg=$i"; done
 
 .. note::
 
@@ -46,8 +50,8 @@ Next time you want to use `grudge`, just run the following command::
 
 You may also like to add this to a startup file (like :file:`$HOME/.bashrc`) or create an alias for it.
 
-After this, you should be able to run the `tests <https://gitlab.tiker.net/inducer/grudge/tree/master/test>`_
-or `examples <https://gitlab.tiker.net/inducer/grudge/tree/master/examples>`_.
+After this, you should be able to run the `tests <https://gitlab.tiker.net/inducer/grudge/tree/main/test>`_
+or `examples <https://gitlab.tiker.net/inducer/grudge/tree/main/examples>`_.
 
 Troubleshooting the Installation
 --------------------------------
@@ -100,7 +104,7 @@ Licensing
 
 :mod:`grudge` is licensed to you under the MIT/X Consortium license:
 
-Copyright (c) 2014-16 Andreas Klöckner and Contributors.
+Copyright (c) 2014-21 Andreas Klöckner and Contributors.
 
 Permission is hereby granted, free of charge, to any person
 obtaining a copy of this software and associated documentation
@@ -138,4 +142,3 @@ AK also gratefully acknowledges a hardware gift from Nvidia Corporation.
 
 The views and opinions expressed herein do not necessarily reflect those of the
 funding agencies.
-

doc/models.rst

+Discontinuous Galerkin Models
+=============================
+
+Compressible Euler Equations
+----------------------------
+
+.. automodule:: grudge.models.euler

doc/operators.rst

+Discontinuous Galerkin operators
+================================
+
+.. automodule:: grudge.op
+.. automodule:: grudge.trace_pair
+
+
+Transferring data between discretizations
+=========================================
+
+.. automodule:: grudge.projection
+
+Reductions
+==========
+
+.. automodule:: grudge.reductions

doc/references.rst

+References
+==========
+
+..
+    When adding references here, please use the demonstrated format:
+    [FirstAuthor_pubyear]
+
+.. [Shewchuk_2002] J. R. Shewchuk (2002),
+    *What Is a Good Linear Element? - Interpolation, Conditioning, Quality Measures*,
+    In Proceedings of the 11th International Meshing Roundtable.
+    `(link) <https://hal.science/hal-04614934/>`__
+    `(slides) <https://people.eecs.berkeley.edu/~jrs/papers/elemtalk.pdf>`__
+
+.. [Hesthaven_2008] J. S. Hesthaven and T. Warburton (2008),
+    *Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications*,
+    Springer.
+    `(doi) <https://doi.org/10.1007/978-0-387-72067-8>`__
+
+.. [Chan_2016] J. Chan, R. J. Hewett, and T. Warburton (2016),
+    *Weight-Adjusted Discontinuous Galerkin Methods: Curvilinear Meshes*,
+    SIAM Journal on Scientific Computing.
+    `(doi) <https://doi.org/10.1137/16M1089198>`__

doc/symbolic.rst

-Symbolic Operator Representation
-================================
-
-Based on :mod:`pymbolic`.
-
-Basic Objects
--------------
-
-.. automodule:: grudge.symbolic.primitives
-
-Operators
----------
-
-.. automodule:: grudge.symbolic.operators

doc/upload-docs.sh

 #! /bin/sh
 
-rsync --verbose --archive --delete _build/html/* doc-upload:doc/grudge
+rsync --verbose --archive --delete _build/html/ doc-upload:doc/grudge

doc/utils.rst

+Helper functions
+================
+
+Estimating stable time-steps
+----------------------------
+
+.. automodule:: grudge.dt_utils
+
+Array contexts
+--------------
+
+.. automodule:: grudge.array_context
+
+Miscellaneous tools
+-------------------
+
+.. automodule:: grudge.tools

examples/advection/surface.py

-__copyright__ = "Copyright (C) 2020 Alexandru Fikl"
+"""Minimal example of a grudge driver for DG on surfaces."""
+
+__copyright__ = """
+Copyright (C) 2020 Alexandru Fikl
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -20,28 +25,33 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+import logging
 import os
 
 import numpy as np
+
 import pyopencl as cl
+import pyopencl.tools as cl_tools
+from arraycontext import flatten
+from meshmode.discretization.connection import FACE_RESTR_INTERIOR
+from pytools.obj_array import make_obj_array
 
-from meshmode.array_context import PyOpenCLArrayContext
-from meshmode.dof_array import thaw, flatten
+import grudge.dof_desc as dof_desc
+import grudge.geometry as geo
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext
 
-from grudge import bind, sym
-from pytools.obj_array import make_obj_array
 
-import logging
 logger = logging.getLogger(__name__)
 
 
 # {{{ plotting (keep in sync with `var-velocity.py`)
 
 class Plotter:
-    def __init__(self, actx, discr, order, visualize=True):
+    def __init__(self, actx, dcoll, order, visualize=True):
         self.actx = actx
-        self.ambient_dim = discr.ambient_dim
-        self.dim = discr.dim
+        self.ambient_dim = dcoll.ambient_dim
+        self.dim = dcoll.dim
 
         self.visualize = visualize
         if not self.visualize:
@@ -51,11 +61,11 @@ class Plotter:
             import matplotlib.pyplot as pt
             self.fig = pt.figure(figsize=(8, 8), dpi=300)
 
-            x = thaw(actx, discr.discr_from_dd(sym.DD_VOLUME).nodes())
-            self.x = actx.to_numpy(flatten(actx.np.atan2(x[1], x[0])))
+            x = actx.thaw(dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL).nodes())
+            self.x = actx.to_numpy(flatten(actx.np.arctan2(x[1], x[0]), self.actx))
         elif self.ambient_dim == 3:
             from grudge.shortcuts import make_visualizer
-            self.vis = make_visualizer(discr, vis_order=order)
+            self.vis = make_visualizer(dcoll)
         else:
             raise ValueError("unsupported dimension")
 
@@ -64,12 +74,12 @@ class Plotter:
             return
 
         if self.ambient_dim == 2:
-            u = self.actx.to_numpy(flatten(evt.state_component))
+            u = self.actx.to_numpy(flatten(evt.state_component, self.actx))
 
-            filename = "%s.png" % basename
+            filename = f"{basename}.png"
             if not overwrite and os.path.exists(filename):
                 from meshmode import FileExistsError
-                raise FileExistsError("output file '%s' already exists" % filename)
+                raise FileExistsError(f"output file '{filename}' already exists")
 
             ax = self.fig.gca()
             ax.grid()
@@ -82,7 +92,7 @@ class Plotter:
             self.fig.savefig(filename)
             self.fig.clf()
         elif self.ambient_dim == 3:
-            self.vis.write_vtk_file("%s.vtu" % basename, [
+            self.vis.write_vtk_file(f"{basename}.vtu", [
                 ("u", evt.state_component)
                 ], overwrite=overwrite)
         else:
@@ -91,10 +101,12 @@ class Plotter:
 # }}}
 
 
-def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
+def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
+
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
 
     # {{{ parameters
 
@@ -103,16 +115,9 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
     # sphere resolution
     resolution = 64 if dim == 2 else 1
 
-    # cfl
-    dt_factor = 2.0
     # final time
     final_time = np.pi
 
-    # velocity field
-    sym_x = sym.nodes(dim)
-    c = make_obj_array([
-        -sym_x[1], sym_x[0], 0.0
-        ])[:dim]
     # flux
     flux_type = "lf"
 
@@ -121,7 +126,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
     # {{{ discretization
 
     if dim == 2:
-        from meshmode.mesh.generation import make_curve_mesh, ellipse
+        from meshmode.mesh.generation import ellipse, make_curve_mesh
         mesh = make_curve_mesh(
                 lambda t: radius * ellipse(1.0, t),
                 np.linspace(0.0, 1.0, resolution + 1),
@@ -133,94 +138,108 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
     else:
         raise ValueError("unsupported dimension")
 
-    quad_tag_to_group_factory = {}
-    if product_tag == "none":
-        product_tag = None
+    discr_tag_to_group_factory = {}
+    if use_quad:
+        qtag = dof_desc.DISCR_TAG_QUAD
+    else:
+        qtag = None
+
+    from meshmode.discretization.poly_element import (
+        QuadratureSimplexGroupFactory,
+        default_simplex_group_factory,
+    )
 
-    from meshmode.discretization.poly_element import \
-            PolynomialWarpAndBlendGroupFactory, \
-            QuadratureSimplexGroupFactory
+    discr_tag_to_group_factory[dof_desc.DISCR_TAG_BASE] = \
+        default_simplex_group_factory(base_dim=dim-1, order=order)
 
-    quad_tag_to_group_factory[sym.QTAG_NONE] = \
-            PolynomialWarpAndBlendGroupFactory(order)
+    if use_quad:
+        discr_tag_to_group_factory[qtag] = \
+            QuadratureSimplexGroupFactory(order=4*order)
 
-    if product_tag:
-        quad_tag_to_group_factory[product_tag] = \
-                QuadratureSimplexGroupFactory(order=4*order)
+    from grudge.discretization import make_discretization_collection
 
-    from grudge import DGDiscretizationWithBoundaries
-    discr = DGDiscretizationWithBoundaries(actx, mesh,
-            quad_tag_to_group_factory=quad_tag_to_group_factory)
+    dcoll = make_discretization_collection(
+        actx, mesh,
+        discr_tag_to_group_factory=discr_tag_to_group_factory
+    )
 
-    volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
+    volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
     logger.info("ndofs:     %d", volume_discr.ndofs)
     logger.info("nelements: %d", volume_discr.mesh.nelements)
 
     # }}}
 
-    # {{{ symbolic operators
+    # {{{ Surface advection operator
+
+    # velocity field
+    x = actx.thaw(dcoll.nodes())
+    c = make_obj_array([-x[1], x[0], 0.0])[:dim]
 
     def f_initial_condition(x):
         return x[0]
 
     from grudge.models.advection import SurfaceAdvectionOperator
-    op = SurfaceAdvectionOperator(c,
+    adv_operator = SurfaceAdvectionOperator(
+        dcoll,
+        c,
         flux_type=flux_type,
-        quad_tag=product_tag)
+        quad_tag=qtag
+    )
 
-    bound_op = bind(discr, op.sym_operator())
-    u0 = bind(discr, f_initial_condition(sym_x))(actx, t=0)
+    u0 = f_initial_condition(x)
 
     def rhs(t, u):
-        return bound_op(actx, t=t, u=u)
+        return adv_operator.operator(t, u)
 
     # check velocity is tangential
-    sym_normal = sym.surface_normal(dim, dim=dim - 1, dd=sym.DD_VOLUME).as_vector()
-    error = bind(discr, sym.norm(2, c.dot(sym_normal)))(actx)
+    from grudge.geometry import normal
+
+    surf_normal = normal(actx, dcoll, dd=dof_desc.DD_VOLUME_ALL)
+
+    error = op.norm(dcoll, c.dot(surf_normal), 2)
     logger.info("u_dot_n:   %.5e", error)
 
     # }}}
 
     # {{{ time stepping
 
-    # compute time step
-    h_min = bind(discr,
-            sym.h_max_from_volume(discr.ambient_dim, dim=discr.dim))(actx)
-    dt = dt_factor * h_min/order**2
+    # FIXME: dt estimate is not necessarily valid for surfaces
+    dt = actx.to_numpy(
+        0.45 * adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u0))
     nsteps = int(final_time // dt) + 1
-    dt = final_time/nsteps + 1.0e-15
 
     logger.info("dt:        %.5e", dt)
     logger.info("nsteps:    %d", nsteps)
 
     from grudge.shortcuts import set_up_rk4
-    dt_stepper = set_up_rk4("u", dt, u0, rhs)
-    plot = Plotter(actx, discr, order, visualize=visualize)
+    dt_stepper = set_up_rk4("u", float(dt), u0, rhs)
+    plot = Plotter(actx, dcoll, order, visualize=visualize)
 
-    norm = bind(discr, sym.norm(2, sym.var("u")))
-    norm_u = norm(actx, u=u0)
+    norm_u = actx.to_numpy(op.norm(dcoll, u0, 2))
 
     step = 0
 
     event = dt_stepper.StateComputed(0.0, 0, 0, u0)
-    plot(event, "fld-surface-%04d" % 0)
+    plot(event, f"fld-surface-{0:04d}")
 
     if visualize and dim == 3:
         from grudge.shortcuts import make_visualizer
-        vis = make_visualizer(discr, vis_order=order)
-        vis.write_vtk_file("fld-surface-velocity.vtu", [
-            ("u", bind(discr, c)(actx)),
-            ("n", bind(discr, sym_normal)(actx))
-            ], overwrite=True)
-
-        df = sym.DOFDesc(sym.FACE_RESTR_INTERIOR)
-        face_discr = discr.connection_from_dds(sym.DD_VOLUME, df).to_discr
-
-        face_normal = bind(discr, sym.normal(
-            df, face_discr.ambient_dim, dim=face_discr.dim))(actx)
+        vis = make_visualizer(dcoll)
+        vis.write_vtk_file(
+            "fld-surface-velocity.vtu",
+            [
+                ("u", c),
+                ("n", surf_normal)
+            ],
+            overwrite=True
+        )
+
+        df = dof_desc.as_dofdesc(FACE_RESTR_INTERIOR)
+        face_discr = dcoll.discr_from_dd(df)
+        face_normal = geo.normal(actx, dcoll, dd=df)
 
         from meshmode.discretization.visualization import make_visualizer
-        vis = make_visualizer(actx, face_discr, vis_order=order)
+        vis = make_visualizer(actx, face_discr)
         vis.write_vtk_file("fld-surface-face-normals.vtu", [
             ("n", face_normal)
             ], overwrite=True)
@@ -231,12 +250,14 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
 
         step += 1
         if step % 10 == 0:
-            norm_u = norm(actx, u=event.state_component)
-            plot(event, "fld-surface-%04d" % step)
+            norm_u = actx.to_numpy(op.norm(dcoll, event.state_component, 2))
+            plot(event, f"fld-surface-{step:04d}")
 
         logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
 
-    plot(event, "fld-surface-%04d" % step)
+        # NOTE: These are here to ensure the solution is bounded for the
+        # time interval specified
+        assert norm_u < 3
 
     # }}}
 
@@ -246,10 +267,14 @@ if __name__ == "__main__":
 
     parser = argparse.ArgumentParser()
     parser.add_argument("--dim", choices=[2, 3], default=2, type=int)
-    parser.add_argument("--qtag", choices=["none", "product"], default="none")
+    parser.add_argument("--order", default=4, type=int)
+    parser.add_argument("--use-quad", action="store_true")
+    parser.add_argument("--visualize", action="store_true")
     args = parser.parse_args()
 
     logging.basicConfig(level=logging.INFO)
     main(cl.create_some_context,
             dim=args.dim,
-            product_tag=args.qtag)
+            order=args.order,
+            use_quad=args.use_quad,
+            visualize=args.visualize)

examples/advection/var-velocity.py

-__copyright__ = "Copyright (C) 2017 Bogdan Enache"
+__copyright__ = """
+Copyright (C) 2017 Bogdan Enache
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -20,27 +23,31 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+import logging
 import os
+
 import numpy as np
 
 import pyopencl as cl
+import pyopencl.tools as cl_tools
+from arraycontext import flatten
+from meshmode.mesh import BTAG_ALL
+from pytools.obj_array import flat_obj_array
 
-from meshmode.array_context import PyOpenCLArrayContext
-from meshmode.dof_array import thaw, flatten
+import grudge.dof_desc as dof_desc
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext
 
-from grudge import bind, sym
-from pytools.obj_array import flat_obj_array
 
-import logging
 logger = logging.getLogger(__name__)
 
 
 # {{{ plotting (keep in sync with `weak.py`)
 
 class Plotter:
-    def __init__(self, actx, discr, order, visualize=True, ylim=None):
+    def __init__(self, actx, dcoll, order, visualize=True, ylim=None):
         self.actx = actx
-        self.dim = discr.ambient_dim
+        self.dim = dcoll.ambient_dim
 
         self.visualize = visualize
         if not self.visualize:
@@ -51,23 +58,23 @@ class Plotter:
             self.fig = pt.figure(figsize=(8, 8), dpi=300)
             self.ylim = ylim
 
-            volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
-            self.x = actx.to_numpy(flatten(thaw(actx, volume_discr.nodes()[0])))
+            volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
+            self.x = actx.to_numpy(flatten(volume_discr.nodes()[0], self.actx))
         else:
             from grudge.shortcuts import make_visualizer
-            self.vis = make_visualizer(discr, vis_order=order)
+            self.vis = make_visualizer(dcoll)
 
     def __call__(self, evt, basename, overwrite=True):
         if not self.visualize:
             return
 
         if self.dim == 1:
-            u = self.actx.to_numpy(flatten(evt.state_component))
+            u = self.actx.to_numpy(flatten(evt.state_component, self.actx))
 
-            filename = "%s.png" % basename
+            filename = f"{basename}.png"
             if not overwrite and os.path.exists(filename):
                 from meshmode import FileExistsError
-                raise FileExistsError("output file '%s' already exists" % filename)
+                raise FileExistsError(f"output file '{filename}' already exists")
 
             ax = self.fig.gca()
             ax.plot(self.x, u, "-")
@@ -81,17 +88,20 @@ class Plotter:
             self.fig.savefig(filename)
             self.fig.clf()
         else:
-            self.vis.write_vtk_file("%s.vtu" % basename, [
+            self.vis.write_vtk_file(f"{basename}.vtu", [
                 ("u", evt.state_component)
                 ], overwrite=overwrite)
 
 # }}}
 
 
-def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
+def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False,
+        flux_type="upwind"):
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
+
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
 
     # {{{ parameters
 
@@ -99,30 +109,14 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
     d = 1.0
     # number of points in each dimension
     npoints = 25
-    # grid spacing
-    h = d / npoints
-
-    # cfl
-    dt_factor = 1.0
-    # finale time
-    final_time = 0.5
-    # time steps
-    dt = dt_factor * h/order**2
-    nsteps = int(final_time // dt) + 1
-    dt = final_time/nsteps + 1.0e-15
-
-    # flux
-    flux_type = "upwind"
-    # velocity field
-    sym_x = sym.nodes(dim)
-    if dim == 1:
-        c = sym_x
+
+    # final time
+    final_time = 1
+
+    if use_quad:
+        qtag = dof_desc.DISCR_TAG_QUAD
     else:
-        # solid body rotation
-        c = flat_obj_array(
-                np.pi * (d/2 - sym_x[1]),
-                np.pi * (sym_x[0] - d/2),
-                0)[:dim]
+        qtag = None
 
     # }}}
 
@@ -131,77 +125,98 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False):
     from meshmode.mesh.generation import generate_regular_rect_mesh
     mesh = generate_regular_rect_mesh(
             a=(0,)*dim, b=(d,)*dim,
-            n=(npoints,)*dim,
+            npoints_per_axis=(npoints,)*dim,
             order=order)
 
-    from meshmode.discretization.poly_element import \
-            QuadratureSimplexGroupFactory
+    from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
 
-    if product_tag:
-        quad_tag_to_group_factory = {
-                product_tag: QuadratureSimplexGroupFactory(order=4*order)
-                }
+    if use_quad:
+        discr_tag_to_group_factory = {
+            qtag: QuadratureSimplexGroupFactory(order=4*order)
+        }
     else:
-        quad_tag_to_group_factory = {}
+        discr_tag_to_group_factory = {}
 
-    from grudge import DGDiscretizationWithBoundaries
-    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order,
-            quad_tag_to_group_factory=quad_tag_to_group_factory)
+    from grudge.discretization import make_discretization_collection
+
+    dcoll = make_discretization_collection(
+        actx, mesh, order=order,
+        discr_tag_to_group_factory=discr_tag_to_group_factory
+    )
 
     # }}}
 
-    # {{{ symbolic operators
+    # {{{ advection operator
 
     # gaussian parameters
-    source_center = np.array([0.5, 0.75, 0.0])[:dim]
-    source_width = 0.05
-    dist_squared = np.dot(sym_x - source_center, sym_x - source_center)
-
-    def f_gaussian(x):
-        return sym.exp(-dist_squared / source_width**2)
 
-    def f_step(x):
-        return sym.If(sym.Comparison(
-            dist_squared, "<", (4*source_width) ** 2),
-            1, 0)
+    def f_halfcircle(x):
+        source_center = np.array([d/2, d/2, d/2])[:dim]
+        dist = x - source_center
+        return (
+                (0.5+0.5*actx.np.tanh(500*(-np.dot(dist, dist) + 0.4**2)))
+                * (0.5+0.5*actx.np.tanh(500*(dist[0]))))
 
-    def u_bc(x):
-        return 0.0
+    def zero_inflow_bc(dtag, t=0):
+        dd = dof_desc.as_dofdesc(dtag, qtag)
+        return dcoll.discr_from_dd(dd).zeros(actx)
 
     from grudge.models.advection import VariableCoefficientAdvectionOperator
-    op = VariableCoefficientAdvectionOperator(
-            c,
-            u_bc(sym.nodes(dim, sym.BTAG_ALL)),
-            quad_tag=product_tag,
-            flux_type=flux_type)
 
-    bound_op = bind(discr, op.sym_operator())
-    u = bind(discr, f_gaussian(sym.nodes(dim)))(actx, t=0)
+    x = actx.thaw(dcoll.nodes())
+
+    # velocity field
+    if dim == 1:
+        c = x
+    else:
+        # solid body rotation
+        c = flat_obj_array(
+            np.pi * (d/2 - x[1]),
+            np.pi * (x[0] - d/2),
+            0
+        )[:dim]
+
+    adv_operator = VariableCoefficientAdvectionOperator(
+        dcoll,
+        c,
+        inflow_u=lambda t: zero_inflow_bc(BTAG_ALL, t),
+        quad_tag=qtag,
+        flux_type=flux_type
+    )
+
+    u = f_halfcircle(x)
 
     def rhs(t, u):
-        return bound_op(t=t, u=u)
+        return adv_operator.operator(t, u)
+
+    dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u))
+
+    logger.info("Timestep size: %g", dt)
 
     # }}}
 
     # {{{ time stepping
 
     from grudge.shortcuts import set_up_rk4
-    dt_stepper = set_up_rk4("u", dt, u, rhs)
-    plot = Plotter(actx, discr, order, visualize=visualize,
+    dt_stepper = set_up_rk4("u", float(dt), u, rhs)
+    plot = Plotter(actx, dcoll, order, visualize=visualize,
             ylim=[-0.1, 1.1])
 
     step = 0
-    norm = bind(discr, sym.norm(2, sym.var("u")))
     for event in dt_stepper.run(t_end=final_time):
         if not isinstance(event, dt_stepper.StateComputed):
             continue
 
         if step % 10 == 0:
-            norm_u = norm(u=event.state_component)
-            plot(event, "fld-var-velocity-%04d" % step)
+            norm_u = actx.to_numpy(op.norm(dcoll, event.state_component, 2))
+            plot(event, f"fld-var-velocity-{step:04d}")
+
+            logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
+            # NOTE: These are here to ensure the solution is bounded for the
+            # time interval specified
+            assert norm_u < 1
 
         step += 1
-        logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
 
     # }}}
 
@@ -211,10 +226,18 @@ if __name__ == "__main__":
 
     parser = argparse.ArgumentParser()
     parser.add_argument("--dim", default=2, type=int)
-    parser.add_argument("--qtag", default="product")
+    parser.add_argument("--order", default=4, type=int)
+    parser.add_argument("--use-quad", action="store_true")
+    parser.add_argument("--visualize", action="store_true")
+    parser.add_argument("--flux", default="upwind",
+            help="'central' or 'upwind'. Run with central to observe aliasing "
+            "instability. Add --use-quad to fix that instability.")
     args = parser.parse_args()
 
     logging.basicConfig(level=logging.INFO)
     main(cl.create_some_context,
-            dim=args.dim,
-            product_tag=args.qtag)
+         dim=args.dim,
+         order=args.order,
+         use_quad=args.use_quad,
+         visualize=args.visualize,
+         flux_type=args.flux)

examples/advection/weak.py

-__copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
+__copyright__ = """
+Copyright (C) 2007 Andreas Kloeckner
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -20,27 +23,31 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+import logging
 import os
+
 import numpy as np
 import numpy.linalg as la
 
 import pyopencl as cl
+import pyopencl.tools as cl_tools
+from arraycontext import flatten
+from meshmode.mesh import BTAG_ALL
 
-from meshmode.array_context import PyOpenCLArrayContext
-from meshmode.dof_array import thaw, flatten
+import grudge.dof_desc as dof_desc
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext
 
-from grudge import bind, sym
 
-import logging
 logger = logging.getLogger(__name__)
 
 
 # {{{ plotting (keep in sync with `var-velocity.py`)
 
 class Plotter:
-    def __init__(self, actx, discr, order, visualize=True, ylim=None):
+    def __init__(self, actx, dcoll, order, visualize=True, ylim=None):
         self.actx = actx
-        self.dim = discr.ambient_dim
+        self.dim = dcoll.ambient_dim
 
         self.visualize = visualize
         if not self.visualize:
@@ -51,23 +58,23 @@ class Plotter:
             self.fig = pt.figure(figsize=(8, 8), dpi=300)
             self.ylim = ylim
 
-            volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
-            self.x = actx.to_numpy(flatten(thaw(actx, volume_discr.nodes()[0])))
+            volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
+            self.x = actx.to_numpy(flatten(volume_discr.nodes()[0], self.actx))
         else:
             from grudge.shortcuts import make_visualizer
-            self.vis = make_visualizer(discr, vis_order=order)
+            self.vis = make_visualizer(dcoll)
 
     def __call__(self, evt, basename, overwrite=True):
         if not self.visualize:
             return
 
         if self.dim == 1:
-            u = self.actx.to_numpy(flatten(evt.state_component))
+            u = self.actx.to_numpy(flatten(evt.state_component, self.actx))
 
-            filename = "%s.png" % basename
+            filename = f"{basename}.png"
             if not overwrite and os.path.exists(filename):
                 from meshmode import FileExistsError
-                raise FileExistsError("output file '%s' already exists" % filename)
+                raise FileExistsError(f"output file '{filename}' already exists")
 
             ax = self.fig.gca()
             ax.plot(self.x, u, "-")
@@ -82,7 +89,7 @@ class Plotter:
             self.fig.savefig(filename)
             self.fig.clf()
         else:
-            self.vis.write_vtk_file("%s.vtu" % basename, [
+            self.vis.write_vtk_file(f"{basename}.vtu", [
                 ("u", evt.state_component)
                 ], overwrite=overwrite)
 
@@ -92,7 +99,9 @@ class Plotter:
 def main(ctx_factory, dim=2, order=4, visualize=False):
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
+
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
 
     # {{{ parameters
 
@@ -100,21 +109,14 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
     d = 1.0
     # number of points in each dimension
     npoints = 20
-    # grid spacing
-    h = d / npoints
 
-    # cfl
-    dt_factor = 2.0
     # final time
     final_time = 1.0
-    # compute number of steps
-    dt = dt_factor * h/order**2
-    nsteps = int(final_time // dt) + 1
-    dt = final_time/nsteps + 1.0e-15
 
     # velocity field
     c = np.array([0.5] * dim)
     norm_c = la.norm(c)
+
     # flux
     flux_type = "central"
 
@@ -127,42 +129,51 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
             [np.linspace(-d/2, d/2, npoints) for _ in range(dim)],
             order=order)
 
-    from grudge import DGDiscretizationWithBoundaries
-    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
+    from grudge.discretization import make_discretization_collection
+
+    dcoll = make_discretization_collection(actx, mesh, order=order)
 
     # }}}
 
-    # {{{ symbolic operators
+    # {{{ weak advection operator
 
     def f(x):
-        return sym.sin(3 * x)
+        return actx.np.sin(3 * x)
 
-    def u_analytic(x):
-        t = sym.var("t", sym.DD_SCALAR)
+    def u_analytic(x, t=0):
         return f(-np.dot(c, x) / norm_c + t * norm_c)
 
     from grudge.models.advection import WeakAdvectionOperator
-    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)))(actx, t=0)
+    adv_operator = WeakAdvectionOperator(
+        dcoll,
+        c,
+        inflow_u=lambda t: u_analytic(
+            actx.thaw(dcoll.nodes(dd=BTAG_ALL)),
+            t=t
+        ),
+        flux_type=flux_type
+    )
+
+    nodes = actx.thaw(dcoll.nodes())
+    u = u_analytic(nodes, t=0)
 
     def rhs(t, u):
-        return bound_op(t=t, u=u)
+        return adv_operator.operator(t, u)
+
+    dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u))
+
+    logger.info("Timestep size: %g", dt)
 
     # }}}
 
     # {{{ time stepping
 
     from grudge.shortcuts import set_up_rk4
-    dt_stepper = set_up_rk4("u", dt, u, rhs)
-    plot = Plotter(actx, discr, order, visualize=visualize,
+    dt_stepper = set_up_rk4("u", float(dt), u, rhs)
+    plot = Plotter(actx, dcoll, order, visualize=visualize,
             ylim=[-1.1, 1.1])
 
-    norm = bind(discr, sym.norm(2, sym.var("u")))
-
     step = 0
     norm_u = 0.0
     for event in dt_stepper.run(t_end=final_time):
@@ -170,12 +181,16 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
             continue
 
         if step % 10 == 0:
-            norm_u = norm(u=event.state_component)
-            plot(event, "fld-weak-%04d" % step)
+            norm_u = actx.to_numpy(op.norm(dcoll, event.state_component, 2))
+            plot(event, f"fld-weak-{step:04d}")
 
         step += 1
         logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
 
+        # NOTE: These are here to ensure the solution is bounded for the
+        # time interval specified
+        assert norm_u < 1
+
     # }}}
 
 
@@ -184,8 +199,12 @@ if __name__ == "__main__":
 
     parser = argparse.ArgumentParser()
     parser.add_argument("--dim", default=2, type=int)
+    parser.add_argument("--order", default=4, type=int)
+    parser.add_argument("--visualize", action="store_true")
     args = parser.parse_args()
 
     logging.basicConfig(level=logging.INFO)
     main(cl.create_some_context,
-            dim=args.dim)
+         dim=args.dim,
+         order=args.order,
+         visualize=args.visualize)

examples/euler/acoustic_pulse.py

+__copyright__ = """
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import logging
+
+import numpy as np
+
+import pyopencl as cl
+import pyopencl.tools as cl_tools
+from arraycontext import ArrayContext
+from meshmode.mesh import BTAG_ALL
+from pytools.obj_array import make_obj_array
+
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext, PytatoPyOpenCLArrayContext
+from grudge.models.euler import ConservedEulerField, EulerOperator, InviscidWallBC
+from grudge.shortcuts import rk4_step
+
+
+logger = logging.getLogger(__name__)
+
+
+def gaussian_profile(
+        actx: ArrayContext,
+        x_vec, t=0, rho0=1.0, rhoamp=1.0, p0=1.0, gamma=1.4,
+        center=None, velocity=None):
+
+    dim = len(x_vec)
+    if center is None:
+        center = np.zeros(shape=(dim,))
+    if velocity is None:
+        velocity = np.zeros(shape=(dim,))
+
+    lump_loc = center + t * velocity
+
+    # coordinates relative to lump center
+    rel_center = make_obj_array(
+        [x_vec[i] - lump_loc[i] for i in range(dim)]
+    )
+    r = actx.np.sqrt(np.dot(rel_center, rel_center))
+    expterm = rhoamp * actx.np.exp(1 - r ** 2)
+
+    mass = expterm + rho0
+    mom = velocity.astype(object) * mass
+    energy = (p0 / (gamma - 1.0)) + np.dot(mom, mom) / (2.0 * mass)
+
+    return ConservedEulerField(mass=mass, energy=energy, momentum=mom)
+
+
+def make_pulse(amplitude, r0, w, r):
+    dim = len(r)
+    r_0 = np.zeros(dim)
+    r_0 = r_0 + r0
+    rel_center = make_obj_array(
+        [r[i] - r_0[i] for i in range(dim)]
+    )
+    actx = r[0].array_context
+    rms2 = w * w
+    r2 = np.dot(rel_center, rel_center) / rms2
+    return amplitude * actx.np.exp(-.5 * r2)
+
+
+def acoustic_pulse_condition(actx: ArrayContext, x_vec, t=0):
+    dim = len(x_vec)
+    vel = np.zeros(shape=(dim,))
+    orig = np.zeros(shape=(dim,))
+    uniform_gaussian = gaussian_profile(
+        actx,
+        x_vec, t=t, center=orig, velocity=vel, rhoamp=0.0)
+
+    amplitude = 1.0
+    width = 0.1
+    pulse = make_pulse(amplitude, orig, width, x_vec)
+
+    return ConservedEulerField(
+        mass=uniform_gaussian.mass,
+        energy=uniform_gaussian.energy + pulse,
+        momentum=uniform_gaussian.momentum
+    )
+
+
+def run_acoustic_pulse(actx,
+                       order=3,
+                       final_time=1,
+                       resolution=16,
+                       overintegration=False,
+                       visualize=False):
+
+    # eos-related parameters
+    gamma = 1.4
+
+    # {{{ discretization
+
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+
+    dim = 2
+    box_ll = -0.5
+    box_ur = 0.5
+    mesh = generate_regular_rect_mesh(
+        a=(box_ll,)*dim,
+        b=(box_ur,)*dim,
+        nelements_per_axis=(resolution,)*dim)
+
+    from meshmode.discretization.poly_element import (
+        QuadratureSimplexGroupFactory,
+        default_simplex_group_factory,
+    )
+
+    from grudge.discretization import make_discretization_collection
+    from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD
+
+    exp_name = f"fld-acoustic-pulse-N{order}-K{resolution}"
+    if overintegration:
+        exp_name += "-overintegrated"
+        quad_tag = DISCR_TAG_QUAD
+    else:
+        quad_tag = None
+
+    dcoll = make_discretization_collection(
+        actx, mesh,
+        discr_tag_to_group_factory={
+            DISCR_TAG_BASE: default_simplex_group_factory(
+                base_dim=mesh.dim, order=order),
+            DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order)
+        }
+    )
+
+    # }}}
+
+    # {{{ Euler operator
+
+    euler_operator = EulerOperator(
+        dcoll,
+        bdry_conditions={BTAG_ALL: InviscidWallBC()},
+        flux_type="lf",
+        gamma=gamma,
+        quadrature_tag=quad_tag
+    )
+
+    def rhs(t, q):
+        return euler_operator.operator(actx, t, q)
+
+    compiled_rhs = actx.compile(rhs)
+
+    from grudge.dt_utils import h_min_from_volume
+
+    cfl = 0.125
+    cn = 0.5*(order + 1)**2
+    dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn
+
+    fields = acoustic_pulse_condition(actx, actx.thaw(dcoll.nodes()))
+
+    logger.info("Timestep size: %g", dt)
+
+    # }}}
+
+    from grudge.shortcuts import make_visualizer
+
+    vis = make_visualizer(dcoll)
+
+    # {{{ time stepping
+
+    step = 0
+    t = 0.0
+    while t < final_time:
+        if step % 10 == 0:
+            norm_q = actx.to_numpy(op.norm(dcoll, fields, 2))
+            logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q)
+            if visualize:
+                vis.write_vtk_file(
+                    f"{exp_name}-{step:04d}.vtu",
+                    [
+                        ("rho", fields.mass),
+                        ("energy", fields.energy),
+                        ("momentum", fields.momentum)
+                    ]
+                )
+            assert norm_q < 5
+
+        fields = actx.thaw(actx.freeze(fields))
+        fields = rk4_step(fields, t, dt, compiled_rhs)
+        t += dt
+        step += 1
+
+    # }}}
+
+
+def main(ctx_factory, order=3, final_time=1, resolution=16,
+         overintegration=False, visualize=False, lazy=False):
+    cl_ctx = ctx_factory()
+    queue = cl.CommandQueue(cl_ctx)
+
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    if lazy:
+        actx = PytatoPyOpenCLArrayContext(queue, allocator=allocator)
+    else:
+        actx = PyOpenCLArrayContext(queue, allocator=allocator)
+
+    run_acoustic_pulse(
+        actx,
+        order=order,
+        resolution=resolution,
+        overintegration=overintegration,
+        final_time=final_time,
+        visualize=visualize
+    )
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--order", default=3, type=int)
+    parser.add_argument("--tfinal", default=0.1, type=float)
+    parser.add_argument("--resolution", default=16, type=int)
+    parser.add_argument("--oi", action="store_true",
+                        help="use overintegration")
+    parser.add_argument("--visualize", action="store_true",
+                        help="write out vtk output")
+    parser.add_argument("--lazy", action="store_true",
+                        help="switch to a lazy computation mode")
+    args = parser.parse_args()
+
+    logging.basicConfig(level=logging.INFO)
+    main(cl.create_some_context,
+         order=args.order,
+         final_time=args.tfinal,
+         resolution=args.resolution,
+         overintegration=args.oi,
+         visualize=args.visualize,
+         lazy=args.lazy)

examples/euler/vortex.py

+__copyright__ = """
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import logging
+
+import pyopencl as cl
+import pyopencl.tools as cl_tools
+
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext, PytatoPyOpenCLArrayContext
+from grudge.models.euler import EulerOperator, vortex_initial_condition
+from grudge.shortcuts import rk4_step
+
+
+logger = logging.getLogger(__name__)
+
+
+def run_vortex(actx, order=3, resolution=8, final_time=5,
+               overintegration=False,
+               flux_type="central",
+               visualize=False):
+
+    logger.info(
+        """
+        Isentropic vortex parameters:\n
+        order: %s\n
+        final_time: %s\n
+        resolution: %s\n
+        overintegration: %s\n
+        flux_type: %s\n
+        visualize: %s
+        """,
+        order, final_time, resolution,
+        overintegration, flux_type, visualize
+    )
+
+    # eos-related parameters
+    gamma = 1.4
+
+    # {{{ discretization
+
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+
+    mesh = generate_regular_rect_mesh(
+        a=(0, -5),
+        b=(20, 5),
+        nelements_per_axis=(2*resolution, resolution),
+        periodic=(True, True))
+
+    from meshmode.discretization.poly_element import (
+        QuadratureSimplexGroupFactory,
+        default_simplex_group_factory,
+    )
+
+    from grudge.discretization import make_discretization_collection
+    from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD
+
+    exp_name = f"fld-vortex-N{order}-K{resolution}-{flux_type}"
+
+    if overintegration:
+        exp_name += "-overintegrated"
+        quad_tag = DISCR_TAG_QUAD
+    else:
+        quad_tag = None
+
+    dcoll = make_discretization_collection(
+        actx, mesh,
+        discr_tag_to_group_factory={
+            DISCR_TAG_BASE: default_simplex_group_factory(
+                base_dim=mesh.dim, order=order),
+            DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order)
+        }
+    )
+
+    # }}}
+
+    # {{{ Euler operator
+
+    euler_operator = EulerOperator(
+        dcoll,
+        flux_type=flux_type,
+        gamma=gamma,
+        quadrature_tag=quad_tag
+    )
+
+    def rhs(t, q):
+        return euler_operator.operator(actx, t, q)
+
+    compiled_rhs = actx.compile(rhs)
+
+    fields = vortex_initial_condition(actx.thaw(dcoll.nodes()))
+
+    from grudge.dt_utils import h_min_from_volume
+
+    cfl = 0.01
+    cn = 0.5*(order + 1)**2
+    dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn
+
+    logger.info("Timestep size: %g", dt)
+
+    # }}}
+
+    from grudge.shortcuts import make_visualizer
+
+    vis = make_visualizer(dcoll)
+
+    # {{{ time stepping
+
+    step = 0
+    t = 0.0
+    while t < final_time:
+        if step % 10 == 0:
+            norm_q = actx.to_numpy(op.norm(dcoll, fields, 2))
+            logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q)
+
+            if visualize:
+                vis.write_vtk_file(
+                    f"{exp_name}-{step:04d}.vtu",
+                    [
+                        ("rho", fields.mass),
+                        ("energy", fields.energy),
+                        ("momentum", fields.momentum)
+                    ]
+                )
+            assert norm_q < 200
+
+        fields = actx.thaw(actx.freeze(fields))
+        fields = rk4_step(fields, t, dt, compiled_rhs)
+        t += dt
+        step += 1
+
+    # }}}
+
+
+def main(ctx_factory, order=3, final_time=5, resolution=8,
+         overintegration=False,
+         lf_stabilization=False,
+         visualize=False,
+         lazy=False):
+    cl_ctx = ctx_factory()
+    queue = cl.CommandQueue(cl_ctx)
+
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    if lazy:
+        actx = PytatoPyOpenCLArrayContext(queue, allocator=allocator)
+    else:
+        actx = PyOpenCLArrayContext(queue, allocator=allocator)
+
+    if lf_stabilization:
+        flux_type = "lf"
+    else:
+        flux_type = "central"
+
+    run_vortex(
+        actx,
+        order=order,
+        resolution=resolution,
+        final_time=final_time,
+        overintegration=overintegration,
+        flux_type=flux_type,
+        visualize=visualize
+    )
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--order", default=3, type=int)
+    parser.add_argument("--tfinal", default=0.015, type=float)
+    parser.add_argument("--resolution", default=8, type=int)
+    parser.add_argument("--oi", action="store_true",
+                        help="use overintegration")
+    parser.add_argument("--lf", action="store_true",
+                        help="turn on lax-friedrichs dissipation")
+    parser.add_argument("--visualize", action="store_true",
+                        help="write out vtk output")
+    parser.add_argument("--lazy", action="store_true",
+                        help="switch to a lazy computation mode")
+    args = parser.parse_args()
+
+    logging.basicConfig(level=logging.INFO)
+    main(cl.create_some_context,
+         order=args.order,
+         final_time=args.tfinal,
+         resolution=args.resolution,
+         overintegration=args.oi,
+         lf_stabilization=args.lf,
+         visualize=args.visualize,
+         lazy=args.lazy)

examples/geometry.py

-"""Minimal example of a grudge driver."""
+"""Minimal example of viewing geometric quantities."""
 
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = """
+Copyright (C) 2015 Andreas Kloeckner
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -22,41 +25,38 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-import numpy as np  # noqa
 import pyopencl as cl
-from grudge import sym, bind, DGDiscretizationWithBoundaries, shortcuts
+import pyopencl.tools as cl_tools
 
-from meshmode.array_context import PyOpenCLArrayContext
+from grudge import geometry, shortcuts
+from grudge.array_context import PyOpenCLArrayContext
+from grudge.discretization import make_discretization_collection
 
 
-def main(write_output=True):
+def main(write_output: bool = True) -> None:
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
 
-    from meshmode.mesh.generation import generate_warped_rect_mesh
-    mesh = generate_warped_rect_mesh(dim=2, order=4, n=6)
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
 
-    discr = DGDiscretizationWithBoundaries(actx, mesh, order=4)
+    from meshmode.mesh import BTAG_ALL
+    from meshmode.mesh.generation import generate_warped_rect_mesh
 
-    sym_op = sym.normal(sym.BTAG_ALL, mesh.dim)
-    #sym_op = sym.nodes(mesh.dim, where=sym.BTAG_ALL)
-    print(sym.pretty(sym_op))
-    op = bind(discr, sym_op)
-    print()
-    print(op.eval_code)
+    mesh = generate_warped_rect_mesh(dim=2, order=4, nelements_side=6)
+    dcoll = make_discretization_collection(actx, mesh, order=4)
 
-    vec = op(actx)
+    nodes = actx.thaw(dcoll.nodes())
+    bdry_nodes = actx.thaw(dcoll.nodes(dd=BTAG_ALL))
+    bdry_normals = geometry.normal(actx, dcoll, dd=BTAG_ALL)
 
-    vis = shortcuts.make_visualizer(discr, 4)
-    vis.write_vtk_file("geo.vtu", [
-        ])
+    if write_output:
+        vis = shortcuts.make_visualizer(dcoll)
+        vis.write_vtk_file("geo.vtu", [("nodes", nodes)])
 
-    bvis = shortcuts.make_boundary_visualizer(discr, 4)
-    bvis.write_vtk_file("bgeo.vtu", [
-        ("normals", vec)
-        ])
+        bvis = shortcuts.make_boundary_visualizer(dcoll)
+        bvis.write_vtk_file("bgeo.vtu", [("bdry normals", bdry_normals),
+                                         ("bdry nodes", bdry_nodes)])
 
 
 if __name__ == "__main__":

examples/hello-grudge.py

+# Solves the PDE:
+# \begin{cases}
+#   u_t + 2\pi u_x = 0, \\
+#   u(0, t) = -\sin(2\pi t), \\
+#   u(x, 0) = \sin(x),
+# \end{cases}
+# on the domain $x \in [0, 2\pi]$. We closely follow Chapter 3 of
+# "Nodal Discontinuous Galerkin Methods" by Hesthaven & Warburton.
+
+# BEGINEXAMPLE
+import numpy as np
+
+import pyopencl as cl
+from meshmode.array_context import PyOpenCLArrayContext
+from meshmode.mesh.generation import generate_box_mesh
+
+import grudge.geometry as geo
+import grudge.op as op
+from grudge.discretization import make_discretization_collection
+from grudge.dof_desc import FACE_RESTR_INTERIOR, BoundaryDomainTag, as_dofdesc
+
+
+ctx = cl.create_some_context()
+queue = cl.CommandQueue(ctx)
+actx = PyOpenCLArrayContext(queue)
+
+nel = 10
+coords = np.linspace(0, 2*np.pi, nel)
+mesh = generate_box_mesh((coords,),
+                         boundary_tag_to_face={"left": ["-x"],
+                                               "right": ["+x"]})
+dcoll = make_discretization_collection(actx, mesh, order=1)
+
+
+def initial_condition(x):
+    # 'x' contains ndim arrays.
+    # 'x[0]' gets the first coordinate value of all the nodes
+    return actx.np.sin(x[0])
+
+
+def left_boundary_condition(x, t):
+    return actx.np.sin(x[0] - 2 * np.pi * t)
+
+
+def flux(dcoll, u_tpair):
+    dd = u_tpair.dd
+    velocity = np.array([2 * np.pi])
+    normal = geo.normal(actx, dcoll, dd)
+
+    v_dot_n = np.dot(velocity, normal)
+    u_upwind = actx.np.where(v_dot_n > 0,
+                             u_tpair.int, u_tpair.ext)
+    return u_upwind * v_dot_n
+
+
+vol_discr = dcoll.discr_from_dd("vol")
+left_bndry = as_dofdesc(BoundaryDomainTag("left"))
+right_bndry = as_dofdesc(BoundaryDomainTag("right"))
+
+x_vol = actx.thaw(dcoll.nodes())
+x_bndry = actx.thaw(dcoll.discr_from_dd(left_bndry).nodes())
+
+uh = initial_condition(x_vol)
+
+dt = 0.001
+t = 0
+t_final = 0.5
+
+# timestepper loop
+while t < t_final:
+    # extract the left boundary trace pair
+    lbnd_tpair = op.bv_trace_pair(dcoll,
+                                  dd=left_bndry,
+                                  interior=uh,
+                                  exterior=left_boundary_condition(x_bndry, t))
+    # extract the right boundary trace pair
+    rbnd_tpair = op.bv_trace_pair(dcoll,
+                                  dd=right_bndry,
+                                  interior=uh,
+                                  exterior=op.project(dcoll, "vol",
+                                                      right_bndry, uh))
+    # extract the trace pairs on the interior faces
+    interior_tpair = op.local_interior_trace_pair(dcoll, uh)
+    Su = op.weak_local_grad(dcoll, uh)
+
+    lift = op.face_mass(dcoll,
+                        # left boundary weak-flux terms
+                        op.project(dcoll,
+                                   left_bndry, "all_faces",
+                                   flux(dcoll, lbnd_tpair))
+                        # right boundary weak-flux terms
+                        + op.project(dcoll,
+                                     right_bndry, "all_faces",
+                                     flux(dcoll, rbnd_tpair))
+                        # interior weak-flux terms
+                        + op.project(dcoll,
+                                     FACE_RESTR_INTERIOR, "all_faces",
+                                     flux(dcoll, interior_tpair)))
+
+    duh_by_dt = op.inverse_mass(dcoll,
+                                np.dot([2 * np.pi], Su) - lift)
+
+    # forward euler time step
+    uh = uh + dt * duh_by_dt
+    t += dt
+# ENDEXAMPLE
+
+
+# Plot the solution:
+def u_exact(x, t):
+    return actx.np.sin(x[0] - 2 * np.pi * t)
+
+
+assert op.norm(dcoll,
+               uh - u_exact(x_vol, t_final),
+               p=2) <= 0.1
+import matplotlib.pyplot as plt
+
+
+plt.plot(actx.to_numpy(actx.np.ravel(x_vol[0][0])),
+         actx.to_numpy(actx.np.ravel(uh[0])), label="Numerical")
+plt.plot(actx.to_numpy(actx.np.ravel(x_vol[0][0])),
+         actx.to_numpy(actx.np.ravel(u_exact(x_vol, t_final)[0])), label="Exact")
+plt.xlabel("$x$")
+plt.ylabel("$u$")
+plt.legend()
+plt.show()

examples/maxwell/cavities.py

-"""Minimal example of a grudge driver."""
-
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = """
+Copyright (C) 2015 Andreas Kloeckner
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -23,32 +24,37 @@ THE SOFTWARE.
 """
 
 
-import numpy as np
-import pyopencl as cl
+import logging
 
-from meshmode.array_context import PyOpenCLArrayContext
+import numpy as np
 
-from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, DGDiscretizationWithBoundaries
+import pyopencl as cl
+import pyopencl.tools as cl_tools
 
+from grudge import op
+from grudge.array_context import PyOpenCLArrayContext
+from grudge.discretization import make_discretization_collection
 from grudge.models.em import get_rectangular_cavity_mode
+from grudge.shortcuts import set_up_rk4
 
 
-STEPS = 60
+logger = logging.getLogger(__name__)
 
 
-def main(dims, write_output=True, order=4):
-    cl_ctx = cl.create_some_context()
+def main(ctx_factory, dim=3, order=4, visualize=False):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
+
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
 
     from meshmode.mesh.generation import generate_regular_rect_mesh
     mesh = generate_regular_rect_mesh(
-            a=(0.0,)*dims,
-            b=(1.0,)*dims,
-            n=(5,)*dims)
+            a=(0.0,)*dim,
+            b=(1.0,)*dim,
+            nelements_per_axis=(4,)*dim)
 
-    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
+    dcoll = make_discretization_collection(actx, mesh, order=order)
 
     if 0:
         epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
@@ -60,75 +66,105 @@ def main(dims, write_output=True, order=4):
         mu = 1
 
     from grudge.models.em import MaxwellOperator
-    op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
 
-    if dims == 3:
-        sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
-        fields = bind(discr, sym_mode)(actx, t=0, epsilon=epsilon, mu=mu)
-    else:
-        sym_mode = get_rectangular_cavity_mode(1, (2, 3))
-        fields = bind(discr, sym_mode)(actx, t=0)
+    maxwell_operator = MaxwellOperator(
+        dcoll,
+        epsilon,
+        mu,
+        flux_type=0.5,
+        dimensions=dim
+    )
 
-    # FIXME
-    #dt = op.estimate_rk4_timestep(discr, fields=fields)
+    def cavity_mode(x, t=0):
+        if dim == 3:
+            return get_rectangular_cavity_mode(actx, x, t, 1, (1, 2, 2))
+        else:
+            return get_rectangular_cavity_mode(actx, x, t, 1, (2, 3))
 
-    op.check_bc_coverage(mesh)
+    fields = cavity_mode(actx.thaw(dcoll.nodes()), t=0)
 
-    # print(sym.pretty(op.sym_operator()))
-    bound_op = bind(discr, op.sym_operator())
+    maxwell_operator.check_bc_coverage(mesh)
 
     def rhs(t, w):
-        return bound_op(t=t, w=w)
+        return maxwell_operator.operator(t, w)
 
-    if mesh.dim == 2:
-        dt = 0.004
-    elif mesh.dim == 3:
-        dt = 0.002
+    dt = actx.to_numpy(
+        maxwell_operator.estimate_rk4_timestep(actx, dcoll, fields=fields))
 
     dt_stepper = set_up_rk4("w", dt, fields, rhs)
 
-    final_t = dt * STEPS
-    nsteps = int(final_t/dt)
+    target_steps = 60
+    final_t = dt * target_steps
+    nsteps = int(final_t/dt) + 1
 
-    print("dt=%g nsteps=%d" % (dt, nsteps))
+    logger.info("dt = %g nsteps = %d", dt, nsteps)
 
     from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
+    vis = make_visualizer(dcoll)
 
     step = 0
 
-    norm = bind(discr, sym.norm(2, sym.var("u")))
-
-    from time import time
-    t_last_step = time()
+    def norm(u):
+        return op.norm(dcoll, u, 2)
 
-    e, h = op.split_eh(fields)
+    e, h = maxwell_operator.split_eh(fields)
 
-    if 1:
-        vis.write_vtk_file("fld-cavities-%04d.vtu" % step,
-                [
-                    ("e", e),
-                    ("h", h),
-                    ])
+    if visualize:
+        vis.write_vtk_file(
+            f"fld-cavities-{step:04d}.vtu",
+            [
+                ("e", e),
+                ("h", h),
+            ]
+        )
 
     for event in dt_stepper.run(t_end=final_t):
         if isinstance(event, dt_stepper.StateComputed):
             assert event.component_id == "w"
 
             step += 1
+            e, h = maxwell_operator.split_eh(event.state_component)
+
+            norm_e0 = actx.to_numpy(norm(u=e[0]))
+            norm_e1 = actx.to_numpy(norm(u=e[1]))
+            norm_h0 = actx.to_numpy(norm(u=h[0]))
+            norm_h1 = actx.to_numpy(norm(u=h[1]))
+
+            logger.info(
+                "[%04d] t = %.5f |e0| = %.5e, |e1| = %.5e, |h0| = %.5e, |h1| = %.5e",
+                step, event.t,
+                norm_e0, norm_e1, norm_h0, norm_h1
+            )
 
-            print(step, event.t, norm(u=e[0]), norm(u=e[1]),
-                    norm(u=h[0]), norm(u=h[1]),
-                    time()-t_last_step)
             if step % 10 == 0:
-                e, h = op.split_eh(event.state_component)
-                vis.write_vtk_file("fld-cavities-%04d.vtu" % step,
+                if visualize:
+                    vis.write_vtk_file(
+                        f"fld-cavities-{step:04d}.vtu",
                         [
                             ("e", e),
                             ("h", h),
-                            ])
-            t_last_step = time()
+                        ]
+                    )
+
+            # NOTE: These are here to ensure the solution is bounded for the
+            # time interval specified
+            assert norm_e0 < 0.5
+            assert norm_e1 < 0.5
+            assert norm_h0 < 0.5
+            assert norm_h1 < 0.5
 
 
 if __name__ == "__main__":
-    main(3)
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--dim", default=3, type=int)
+    parser.add_argument("--order", default=4, type=int)
+    parser.add_argument("--visualize", action="store_true")
+    args = parser.parse_args()
+
+    logging.basicConfig(level=logging.INFO)
+    main(cl.create_some_context,
+         dim=args.dim,
+         order=args.order,
+         visualize=args.visualize)

examples/wave/var-propagation-speed.py

-"""Minimal example of a grudge driver."""
-
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = """
+Copyright (C) 2015 Andreas Kloeckner
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -23,104 +24,151 @@ THE SOFTWARE.
 """
 
 
+import logging
+
 import numpy as np
+
 import pyopencl as cl
+import pyopencl.tools as cl_tools
+from pytools.obj_array import flat_obj_array
+
+from grudge import op
+from grudge.array_context import PyOpenCLArrayContext
+from grudge.discretization import make_discretization_collection
 from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, DGDiscretizationWithBoundaries
 
-from meshmode.array_context import PyOpenCLArrayContext
+
+logger = logging.getLogger(__name__)
 
 
-def main(write_output=True, order=4):
-    cl_ctx = cl.create_some_context()
+def main(ctx_factory, dim=2, order=4, visualize=False):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
 
-    dims = 2
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
+
     from meshmode.mesh.generation import generate_regular_rect_mesh
     mesh = generate_regular_rect_mesh(
-            a=(-0.5,)*dims,
-            b=(0.5,)*dims,
-            n=(20,)*dims)
-
-    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
+            a=(-0.5,)*dim,
+            b=(0.5,)*dim,
+            nelements_per_axis=(20,)*dim)
+
+    dcoll = make_discretization_collection(actx, mesh, order=order)
+
+    def source_f(actx, dcoll, t=0):
+        source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
+        source_width = 0.05
+        source_omega = 3
+        nodes = actx.thaw(dcoll.nodes())
+        source_center_dist = flat_obj_array(
+            [nodes[i] - source_center[i] for i in range(dcoll.dim)]
+        )
+        return (
+            np.sin(source_omega*t)
+            * actx.np.exp(
+                -np.dot(source_center_dist, source_center_dist)
+                / source_width**2
+            )
+        )
+
+    x = actx.thaw(dcoll.nodes())
+    ones = dcoll.zeros(actx) + 1
+    c = actx.np.where(np.dot(x, x) < 0.15, 0.1 * ones, 0.2 * ones)
 
-    source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
-    source_width = 0.05
-    source_omega = 3
-
-    sym_x = sym.nodes(mesh.dim)
-    sym_source_center_dist = sym_x - source_center
-    sym_t = sym.ScalarVariable("t")
-    c = sym.If(sym.Comparison(
-                np.dot(sym_x, sym_x), "<", 0.15),
-                np.float32(0.1), np.float32(0.2))
-
-    from grudge.models.wave import VariableCoefficientWeakWaveOperator
     from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    op = VariableCoefficientWeakWaveOperator(c,
-            discr.dim,
-            source_f=(
-                sym.sin(source_omega*sym_t)
-                * sym.exp(
-                    -np.dot(sym_source_center_dist, sym_source_center_dist)
-                    / source_width**2)),
-            dirichlet_tag=BTAG_NONE,
-            neumann_tag=BTAG_NONE,
-            radiation_tag=BTAG_ALL,
-            flux_type="upwind")
 
-    from pytools.obj_array import flat_obj_array
-    fields = flat_obj_array(discr.zeros(actx),
-            [discr.zeros(actx) for i in range(discr.dim)])
+    from grudge.models.wave import VariableCoefficientWeakWaveOperator
 
-    op.check_bc_coverage(mesh)
+    wave_op = VariableCoefficientWeakWaveOperator(
+        dcoll,
+        actx.freeze(c),
+        source_f=source_f,
+        dirichlet_tag=BTAG_NONE,
+        neumann_tag=BTAG_NONE,
+        radiation_tag=BTAG_ALL,
+        flux_type="upwind"
+    )
 
-    c_eval = bind(discr, c)(actx)
+    fields = flat_obj_array(
+        dcoll.zeros(actx),
+        [dcoll.zeros(actx) for i in range(dcoll.dim)]
+    )
 
-    bound_op = bind(discr, op.sym_operator())
+    wave_op.check_bc_coverage(mesh)
 
     def rhs(t, w):
-        return bound_op(t=t, w=w)
-
-    if mesh.dim == 2:
-        dt = 0.04 * 0.3
-    elif mesh.dim == 3:
-        dt = 0.02 * 0.1
+        return wave_op.operator(t, w)
 
+    dt = actx.to_numpy(
+        2/3 * wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields))
     dt_stepper = set_up_rk4("w", dt, fields, rhs)
 
     final_t = 1
-    nsteps = int(final_t/dt)
-    print("dt=%g nsteps=%d" % (dt, nsteps))
+    nsteps = int(final_t/dt) + 1
+
+    logger.info("dt=%g nsteps=%d", dt, nsteps)
 
     from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
+    vis = make_visualizer(dcoll)
 
     step = 0
 
-    norm = bind(discr, sym.norm(2, sym.var("u")))
+    def norm(u):
+        return op.norm(dcoll, u, 2)
 
     from time import time
     t_last_step = time()
 
+    if visualize:
+        u = fields[0]
+        v = fields[1:]
+        vis.write_vtk_file(
+            f"fld-var-propagation-speed-{step:04d}.vtu",
+            [
+                ("u", u),
+                ("v", v),
+                ("c", c),
+            ]
+        )
+
     for event in dt_stepper.run(t_end=final_t):
         if isinstance(event, dt_stepper.StateComputed):
             assert event.component_id == "w"
 
             step += 1
 
-            print(step, event.t, norm(u=event.state_component[0]),
-                    time()-t_last_step)
             if step % 10 == 0:
-                vis.write_vtk_file("fld-var-propogation-speed-%04d.vtu" % step,
+                logger.info("step: %d t: %.8e L2: %.8e",
+                            step, time() - t_last_step,
+                            actx.to_numpy(norm(event.state_component[0])))
+                if visualize:
+                    vis.write_vtk_file(
+                        f"fld-var-propagation-speed-{step:04d}.vtu",
                         [
                             ("u", event.state_component[0]),
                             ("v", event.state_component[1:]),
-                            ("c", c_eval),
-                            ])
+                            ("c", c),
+                        ]
+                    )
             t_last_step = time()
 
+            # NOTE: These are here to ensure the solution is bounded for the
+            # time interval specified
+            assert norm(u=event.state_component[0]) < 1
+
 
 if __name__ == "__main__":
-    main()
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--dim", default=2, type=int)
+    parser.add_argument("--order", default=4, type=int)
+    parser.add_argument("--visualize", action="store_true")
+    args = parser.parse_args()
+
+    logging.basicConfig(level=logging.INFO)
+    main(cl.create_some_context,
+         dim=args.dim,
+         order=args.order,
+         visualize=args.visualize)

examples/wave/wave-min-mpi.py

-"""Minimal example of a grudge driver."""
-
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = """
+Copyright (C) 2015 Andreas Kloeckner
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -23,123 +24,175 @@ THE SOFTWARE.
 """
 
 
+import logging
+
 import numpy as np
+from mpi4py import MPI
+
 import pyopencl as cl
-from meshmode.array_context import PyOpenCLArrayContext
+import pyopencl.tools as cl_tools
+from pytools.obj_array import flat_obj_array
+
+import grudge.op as op
+from grudge import make_discretization_collection
+from grudge.array_context import MPIPyOpenCLArrayContext
 from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-from mpi4py import MPI
 
 
-def main(write_output=True, order=4):
+logger = logging.getLogger(__name__)
+
+
+class WaveTag:
+    pass
+
+
+def main(dim=2, order=4, visualize=True):
+    comm = MPI.COMM_WORLD
+    num_parts = comm.size
+
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
 
-    comm = MPI.COMM_WORLD
-    num_parts = comm.Get_size()
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = MPIPyOpenCLArrayContext(comm, queue, allocator=allocator)
 
-    from meshmode.distributed import MPIMeshDistributor, get_partition_by_pymetis
-    mesh_dist = MPIMeshDistributor(comm)
+    from meshmode.distributed import get_partition_by_pymetis, membership_list_to_map
+    from meshmode.mesh.processing import partition_mesh
 
-    if mesh_dist.is_mananger_rank():
-        dims = 2
+    if comm.rank == 0:
         from meshmode.mesh.generation import generate_regular_rect_mesh
         mesh = generate_regular_rect_mesh(
-                a=(-0.5,)*dims,
-                b=(0.5,)*dims,
-                n=(16,)*dims)
-
-        print("%d elements" % mesh.nelements)
+                a=(-0.5,)*dim,
+                b=(0.5,)*dim,
+                nelements_per_axis=(16,)*dim)
 
-        part_per_element = get_partition_by_pymetis(mesh, num_parts)
+        logger.info("%d elements", mesh.nelements)
 
-        local_mesh = mesh_dist.send_mesh_parts(mesh, part_per_element, num_parts)
+        part_id_to_part = partition_mesh(mesh,
+                       membership_list_to_map(
+                           get_partition_by_pymetis(mesh, num_parts)))
+        parts = [part_id_to_part[i] for i in range(num_parts)]
+        local_mesh = comm.scatter(parts)
 
         del mesh
 
     else:
-        local_mesh = mesh_dist.receive_mesh_part()
+        local_mesh = comm.scatter(None)
+
+    dcoll = make_discretization_collection(actx, local_mesh, order=order)
+
+    def source_f(actx, dcoll, t=0):
+        source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
+        source_width = 0.05
+        source_omega = 3
+        nodes = actx.thaw(dcoll.nodes())
+        source_center_dist = flat_obj_array(
+            [nodes[i] - source_center[i] for i in range(dcoll.dim)]
+        )
+        return (
+            np.sin(source_omega*t)
+            * actx.np.exp(
+                -np.dot(source_center_dist, source_center_dist)
+                / source_width**2
+            )
+        )
 
-    discr = DGDiscretizationWithBoundaries(actx, local_mesh, order=order,
-            mpi_communicator=comm)
-
-    if local_mesh.dim == 2:
-        dt = 0.04
-    elif local_mesh.dim == 3:
-        dt = 0.02
+    from meshmode.mesh import BTAG_ALL, BTAG_NONE
 
-    source_center = np.array([0.1, 0.22, 0.33])[:local_mesh.dim]
-    source_width = 0.05
-    source_omega = 3
+    from grudge.models.wave import WeakWaveOperator
 
-    sym_x = sym.nodes(local_mesh.dim)
-    sym_source_center_dist = sym_x - source_center
-    sym_t = sym.ScalarVariable("t")
+    wave_op = WeakWaveOperator(
+        dcoll,
+        0.1,
+        source_f=source_f,
+        dirichlet_tag=BTAG_NONE,
+        neumann_tag=BTAG_NONE,
+        radiation_tag=BTAG_ALL,
+        flux_type="upwind"
+    )
 
-    from grudge.models.wave import WeakWaveOperator
-    from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    op = WeakWaveOperator(0.1, discr.dim,
-            source_f=(
-                sym.sin(source_omega*sym_t)
-                * sym.exp(
-                    -np.dot(sym_source_center_dist, sym_source_center_dist)
-                    / source_width**2)),
-            dirichlet_tag=BTAG_NONE,
-            neumann_tag=BTAG_NONE,
-            radiation_tag=BTAG_ALL,
-            flux_type="upwind")
-
-    from pytools.obj_array import flat_obj_array
     fields = flat_obj_array(
-            discr.zeros(actx),
-            [discr.zeros(actx) for i in range(discr.dim)])
-
-    # FIXME
-    #dt = op.estimate_rk4_timestep(discr, fields=fields)
+        dcoll.zeros(actx),
+        [dcoll.zeros(actx) for i in range(dcoll.dim)]
+    )
 
-    op.check_bc_coverage(local_mesh)
+    dt = actx.to_numpy(
+        2/3 * wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields))
 
-    # print(sym.pretty(op.sym_operator()))
-    bound_op = bind(discr, op.sym_operator())
+    wave_op.check_bc_coverage(local_mesh)
 
     def rhs(t, w):
-        return bound_op(t=t, w=w)
+        return wave_op.operator(t, w)
 
     dt_stepper = set_up_rk4("w", dt, fields, rhs)
 
     final_t = 10
-    nsteps = int(final_t/dt)
-    print("dt=%g nsteps=%d" % (dt, nsteps))
+    nsteps = int(final_t/dt) + 1
+
+    if comm.rank == 0:
+        logger.info("dt=%g nsteps=%d", dt, nsteps)
 
     from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
+    vis = make_visualizer(dcoll)
 
     step = 0
 
-    norm = bind(discr, sym.norm(2, sym.var("u")))
+    def norm(u):
+        return op.norm(dcoll, u, 2)
 
     from time import time
     t_last_step = time()
 
+    if visualize:
+        u = fields[0]
+        v = fields[1:]
+        vis.write_parallel_vtk_file(
+            comm,
+            f"fld-wave-min-mpi-{{rank:03d}}-{step:04d}.vtu",
+            [
+                ("u", u),
+                ("v", v),
+            ]
+        )
+
     for event in dt_stepper.run(t_end=final_t):
         if isinstance(event, dt_stepper.StateComputed):
             assert event.component_id == "w"
 
             step += 1
+            l2norm = norm(u=event.state_component[0])
 
-            print(step, event.t, norm(u=event.state_component[0]),
-                    time()-t_last_step)
             if step % 10 == 0:
-                vis.write_parallel_vtk_file(
+                if comm.rank == 0:
+                    logger.info("step: %d t: %.8e L2: %.8e",
+                                step, time() - t_last_step, l2norm)
+                if visualize:
+                    vis.write_parallel_vtk_file(
                         comm,
                         f"fld-wave-min-mpi-{{rank:03d}}-{step:04d}.vtu",
                         [
                             ("u", event.state_component[0]),
                             ("v", event.state_component[1:]),
-                            ])
+                        ]
+                    )
             t_last_step = time()
 
+            # NOTE: These are here to ensure the solution is bounded for the
+            # time interval specified
+            assert l2norm < 1
+
 
 if __name__ == "__main__":
-    main()
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--dim", default=2, type=int)
+    parser.add_argument("--order", default=4, type=int)
+    parser.add_argument("--visualize", action="store_true")
+    args = parser.parse_args()
+
+    logging.basicConfig(level=logging.INFO)
+    main(
+         dim=args.dim,
+         order=args.order,
+         visualize=args.visualize)

examples/wave/wave-min.py

-"""Minimal example of a grudge driver."""
-
-
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
-
-__license__ = """
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in
-all copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-THE SOFTWARE.
-"""
-
-
-import numpy as np
-import pyopencl as cl
-from meshmode.array_context import PyOpenCLArrayContext
-from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-
-
-def main(write_output=True, order=4):
-    cl_ctx = cl.create_some_context()
-    queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue)
-
-    dims = 2
-    from meshmode.mesh.generation import generate_regular_rect_mesh
-    mesh = generate_regular_rect_mesh(
-            a=(-0.5,)*dims,
-            b=(0.5,)*dims,
-            n=(16,)*dims)
-
-    if mesh.dim == 2:
-        dt = 0.04
-    elif mesh.dim == 3:
-        dt = 0.02
-
-    print("%d elements" % mesh.nelements)
-
-    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
-
-    source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
-    source_width = 0.05
-    source_omega = 3
-
-    sym_x = sym.nodes(mesh.dim)
-    sym_source_center_dist = sym_x - source_center
-    sym_t = sym.ScalarVariable("t")
-
-    from grudge.models.wave import WeakWaveOperator
-    from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    op = WeakWaveOperator(0.1, discr.dim,
-            source_f=(
-                sym.sin(source_omega*sym_t)
-                * sym.exp(
-                    -np.dot(sym_source_center_dist, sym_source_center_dist)
-                    / source_width**2)),
-            dirichlet_tag=BTAG_NONE,
-            neumann_tag=BTAG_NONE,
-            radiation_tag=BTAG_ALL,
-            flux_type="upwind")
-
-    from pytools.obj_array import flat_obj_array
-    fields = flat_obj_array(discr.zeros(actx),
-            [discr.zeros(actx) for i in range(discr.dim)])
-
-    # FIXME
-    #dt = op.estimate_rk4_timestep(discr, fields=fields)
-
-    op.check_bc_coverage(mesh)
-
-    # print(sym.pretty(op.sym_operator()))
-    bound_op = bind(discr, op.sym_operator())
-
-    def rhs(t, w):
-        return bound_op(t=t, w=w)
-
-    dt_stepper = set_up_rk4("w", dt, fields, rhs)
-
-    final_t = 10
-    nsteps = int(final_t/dt)
-    print("dt=%g nsteps=%d" % (dt, nsteps))
-
-    from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
-
-    step = 0
-
-    norm = bind(discr, sym.norm(2, sym.var("u")))
-
-    from time import time
-    t_last_step = time()
-
-    for event in dt_stepper.run(t_end=final_t):
-        if isinstance(event, dt_stepper.StateComputed):
-            assert event.component_id == "w"
-
-            step += 1
-
-            print(step, event.t, norm(u=event.state_component[0]),
-                    time()-t_last_step)
-            if step % 10 == 0:
-                vis.write_vtk_file("fld-wave-min-%04d.vtu" % step,
-                        [
-                            ("u", event.state_component[0]),
-                            ("v", event.state_component[1:]),
-                            ])
-            t_last_step = time()
-
-
-if __name__ == "__main__":
-    main()