doc/geometry.rst

+Metric terms and transformations
+================================
+
+.. automodule:: grudge.geometry.metrics

doc/index.rst

-.. grudge documentation master file, created by
-   sphinx-quickstart on Sun Sep 27 13:08:30 2015.
-   You can adapt this file completely to your liking, but it should at least
-   contain the root `toctree` directive.
-
-Welcome to grudge's documentation!
+Welcome to grudge's Documentation!
 ==================================
 
+Here's an example to solve the PDE
+
+.. math::
+    \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 :math:`x \in [0, 2\pi]`. We closely follow Chapter 3 of
+[Hesthaven_2008]_.
+
+
+.. literalinclude:: ../examples/hello-grudge.py
+   :start-after: BEGINEXAMPLE
+   :end-before: ENDEXAMPLE
+
+Plotting numerical solution ``uh`` in results in
+
+.. plot:: ../examples/hello-grudge.py
+
 Contents:
 
 .. toctree::
     :maxdepth: 2
 
+    discretization
+    dof_desc
+    geometry
+    operators
+    utils
+    models
+    references
     misc
+    🚀 Github <https://github.com/inducer/grudge>
+    💾 Download Releases <https://pypi.org/project/grudge>
 
-
-
-Indices and tables
+Indices and Tables
 ==================
 
 * :ref:`genindex`

doc/misc.rst

@@ -17,34 +17,41 @@ MacOS support is in the works.
     Everywhere else, just making sure you have the ``g++`` package should be
     enough.
 
-#.  Installing `miniconda for Python 3 on 64-bit Linux <https://conda.io/miniconda.html>`_.
+#.  Install `miniforge <https://github.com/conda-forge/miniforge>`_::
 
-#.  ``export CONDA=/WHERE/YOU/INSTALLED/miniconda3``
+        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``
 
     If you accepted the default location, this should work:
 
-    ``export CONDA=$HOME/miniconda3``
+    ``export CONDA=$HOME/miniforge3``
 
 #.  ``$CONDA/bin/conda create -n dgfem``
 
 #.  ``source $CONDA/bin/activate dgfem``
 
-#.  ``conda config --add channels conda-forge``
-
 #.  ``conda install git pip pocl islpy pyopencl``
 
 #.  Type the following command::
 
-        hash -r; for i in pymbolic cgen genpy modepy pyvisfile loopy meshmode dagrt leap grudge; do python -m pip install git+https://gitlab.tiker.net/inducer/$i.git; 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::
+
+    In each case, you may leave out the ``--editable`` flag if you would not like
+    an editable copy of the source code checked out in a subfolder.
 
 Next time you want to use `grudge`, just run the following command::
 
-    source /WHERE/YOU/INSTALLED/miniconda3/bin/activate dgfem
+    source /WHERE/YOU/INSTALLED/miniforge3/bin/activate dgfem
 
 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
 --------------------------------
@@ -97,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
@@ -123,12 +130,15 @@ OTHER DEALINGS IN THE SOFTWARE.
 Acknowledgments
 ===============
 
-Andreas Klöckner's work on :mod:`grudge` was supported in part by
+Work on grudge was supported in part by
 
-* US Navy ONR grant number N00014-14-1-0117
-* the US National Science Foundation under grant number CCF-1524433.
+* the Department of Energy, National Nuclear Security Administration,
+  under Award Number DE-NA0003963,
+* the US Navy ONR, under grant number N00014-14-1-0117, and
+* the US National Science Foundation under grant numbers CCF-1524433,
+  and OAC-1931577.
 
-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.
+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/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/advection.py

-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-import numpy.linalg as la
-
-
-
-
-def main(write_output=True, flux_type_arg="upwind"):
-    from grudge.tools import mem_checkpoint
-    from math import sin, cos, pi, sqrt
-    from math import floor
-
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context()
-
-    def f(x):
-        return sin(pi*x)
-
-    def u_analytic(x, el, t):
-        return f((-numpy.dot(v, x)/norm_v+t*norm_v))
-
-    def boundary_tagger(vertices, el, face_nr, all_v):
-        if numpy.dot(el.face_normals[face_nr], v) < 0:
-            return ["inflow"]
-        else:
-            return ["outflow"]
-
-    dim = 2
-
-    if dim == 1:
-        v = numpy.array([1])
-        if rcon.is_head_rank:
-            from grudge.mesh.generator import make_uniform_1d_mesh
-            mesh = make_uniform_1d_mesh(0, 2, 10, periodic=True)
-    elif dim == 2:
-        v = numpy.array([2,0])
-        if rcon.is_head_rank:
-            from grudge.mesh.generator import make_disk_mesh
-            mesh = make_disk_mesh(boundary_tagger=boundary_tagger)
-    elif dim == 3:
-        v = numpy.array([0,0,1])
-        if rcon.is_head_rank:
-            from grudge.mesh.generator import make_cylinder_mesh, make_ball_mesh, make_box_mesh
-
-            mesh = make_cylinder_mesh(max_volume=0.04, height=2, boundary_tagger=boundary_tagger,
-                    periodic=False, radial_subdivisions=32)
-    else:
-        raise RuntimeError("bad number of dimensions")
-
-    norm_v = la.norm(v)
-
-    if rcon.is_head_rank:
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    if dim != 1:
-        mesh_data = mesh_data.reordered_by("cuthill")
-
-    discr = rcon.make_discretization(mesh_data, order=4)
-    vis_discr = discr
-
-    from grudge.visualization import VtkVisualizer
-    if write_output:
-        vis = VtkVisualizer(vis_discr, rcon, "fld")
-
-    # operator setup ----------------------------------------------------------
-    from grudge.data import \
-            ConstantGivenFunction, \
-            TimeConstantGivenFunction, \
-            TimeDependentGivenFunction
-    from grudge.models.advection import StrongAdvectionOperator, WeakAdvectionOperator
-    op = WeakAdvectionOperator(v, 
-            inflow_u=TimeDependentGivenFunction(u_analytic),
-            flux_type=flux_type_arg)
-
-    u = discr.interpolate_volume_function(lambda x, el: u_analytic(x, el, 0))
-
-    # timestep setup ----------------------------------------------------------
-    from grudge.timestep.runge_kutta import LSRK4TimeStepper
-    stepper = LSRK4TimeStepper()
-
-    if rcon.is_head_rank:
-        print("%d elements" % len(discr.mesh.elements))
-
-    # diagnostics setup -------------------------------------------------------
-    from pytools.log import LogManager, \
-            add_general_quantities, \
-            add_simulation_quantities, \
-            add_run_info
-
-    if write_output:
-        log_file_name = "advection.dat"
-    else:
-        log_file_name = None
-
-    logmgr = LogManager(log_file_name, "w", rcon.communicator)
-    add_run_info(logmgr)
-    add_general_quantities(logmgr)
-    add_simulation_quantities(logmgr)
-    discr.add_instrumentation(logmgr)
-
-    stepper.add_instrumentation(logmgr)
-
-    from grudge.log import Integral, LpNorm
-    u_getter = lambda: u
-    logmgr.add_quantity(Integral(u_getter, discr, name="int_u"))
-    logmgr.add_quantity(LpNorm(u_getter, discr, p=1, name="l1_u"))
-    logmgr.add_quantity(LpNorm(u_getter, discr, name="l2_u"))
-
-    logmgr.add_watches(["step.max", "t_sim.max", "l2_u", "t_step.max"])
-
-    # timestep loop -----------------------------------------------------------
-    rhs = op.bind(discr)
-
-    try:
-        from grudge.timestep import times_and_steps
-        step_it = times_and_steps(
-                final_time=3, logmgr=logmgr,
-                max_dt_getter=lambda t: op.estimate_timestep(discr,
-                    stepper=stepper, t=t, fields=u))
-
-        for step, t, dt in step_it:
-            if step % 5 == 0 and write_output:
-                visf = vis.make_file("fld-%04d" % step)
-                vis.add_data(visf, [ 
-                    ("u", discr.convert_volume(u, kind="numpy")), 
-                    ], time=t, step=step)
-                visf.close()
-
-            u = stepper(u, t, dt, rhs)
-
-        true_u = discr.interpolate_volume_function(lambda x, el: u_analytic(x, el, t))
-        print(discr.norm(u-true_u))
-        assert discr.norm(u-true_u) < 1e-2
-    finally:
-        if write_output:
-            vis.close()
-
-        logmgr.close()
-        discr.close()
-
-
-
-
-if __name__ == "__main__":
-    main()
-
-
-
-
-# entry points for py.test ----------------------------------------------------
-def test_advection():
-    from pytools.test import mark_test
-    mark_long = mark_test.long
-
-    for flux_type in ["upwind", "central", "lf"]:
-        yield "advection with %s flux" % flux_type, \
-                mark_long(main), False, flux_type

examples/advection/surface.py

+"""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
+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 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
+
+import grudge.dof_desc as dof_desc
+import grudge.geometry as geo
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext
+
+
+logger = logging.getLogger(__name__)
+
+
+# {{{ plotting (keep in sync with `var-velocity.py`)
+
+class Plotter:
+    def __init__(self, actx, dcoll, order, visualize=True):
+        self.actx = actx
+        self.ambient_dim = dcoll.ambient_dim
+        self.dim = dcoll.dim
+
+        self.visualize = visualize
+        if not self.visualize:
+            return
+
+        if self.ambient_dim == 2:
+            import matplotlib.pyplot as pt
+            self.fig = pt.figure(figsize=(8, 8), dpi=300)
+
+            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(dcoll)
+        else:
+            raise ValueError("unsupported dimension")
+
+    def __call__(self, evt, basename, overwrite=True):
+        if not self.visualize:
+            return
+
+        if self.ambient_dim == 2:
+            u = self.actx.to_numpy(flatten(evt.state_component, self.actx))
+
+            filename = f"{basename}.png"
+            if not overwrite and os.path.exists(filename):
+                from meshmode import FileExistsError
+                raise FileExistsError(f"output file '{filename}' already exists")
+
+            ax = self.fig.gca()
+            ax.grid()
+            ax.plot(self.x, u, "-")
+            ax.plot(self.x, u, "k.")
+            ax.set_xlabel(r"$\theta$")
+            ax.set_ylabel("$u$")
+            ax.set_title(f"t = {evt.t:.2f}")
+
+            self.fig.savefig(filename)
+            self.fig.clf()
+        elif self.ambient_dim == 3:
+            self.vis.write_vtk_file(f"{basename}.vtu", [
+                ("u", evt.state_component)
+                ], overwrite=overwrite)
+        else:
+            raise ValueError("unsupported dimension")
+
+# }}}
+
+
+def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
+    cl_ctx = ctx_factory()
+    queue = cl.CommandQueue(cl_ctx)
+
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
+
+    # {{{ parameters
+
+    # sphere radius
+    radius = 1.0
+    # sphere resolution
+    resolution = 64 if dim == 2 else 1
+
+    # final time
+    final_time = np.pi
+
+    # flux
+    flux_type = "lf"
+
+    # }}}
+
+    # {{{ discretization
+
+    if dim == 2:
+        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),
+                order)
+    elif dim == 3:
+        from meshmode.mesh.generation import generate_icosphere
+        mesh = generate_icosphere(radius, order=4 * order,
+                uniform_refinement_rounds=resolution)
+    else:
+        raise ValueError("unsupported dimension")
+
+    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,
+    )
+
+    discr_tag_to_group_factory[dof_desc.DISCR_TAG_BASE] = \
+        default_simplex_group_factory(base_dim=dim-1, order=order)
+
+    if use_quad:
+        discr_tag_to_group_factory[qtag] = \
+            QuadratureSimplexGroupFactory(order=4*order)
+
+    from grudge.discretization import make_discretization_collection
+
+    dcoll = make_discretization_collection(
+        actx, mesh,
+        discr_tag_to_group_factory=discr_tag_to_group_factory
+    )
+
+    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)
+
+    # }}}
+
+    # {{{ 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
+    adv_operator = SurfaceAdvectionOperator(
+        dcoll,
+        c,
+        flux_type=flux_type,
+        quad_tag=qtag
+    )
+
+    u0 = f_initial_condition(x)
+
+    def rhs(t, u):
+        return adv_operator.operator(t, u)
+
+    # check velocity is tangential
+    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
+
+    # 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
+
+    logger.info("dt:        %.5e", dt)
+    logger.info("nsteps:    %d", nsteps)
+
+    from grudge.shortcuts import set_up_rk4
+    dt_stepper = set_up_rk4("u", float(dt), u0, rhs)
+    plot = Plotter(actx, dcoll, order, visualize=visualize)
+
+    norm_u = actx.to_numpy(op.norm(dcoll, u0, 2))
+
+    step = 0
+
+    event = dt_stepper.StateComputed(0.0, 0, 0, u0)
+    plot(event, f"fld-surface-{0:04d}")
+
+    if visualize and dim == 3:
+        from grudge.shortcuts import make_visualizer
+        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.write_vtk_file("fld-surface-face-normals.vtu", [
+            ("n", face_normal)
+            ], overwrite=True)
+
+    for event in dt_stepper.run(t_end=final_time, max_steps=nsteps + 1):
+        if not isinstance(event, dt_stepper.StateComputed):
+            continue
+
+        step += 1
+        if step % 10 == 0:
+            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)
+
+        # NOTE: These are here to ensure the solution is bounded for the
+        # time interval specified
+        assert norm_u < 3
+
+    # }}}
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--dim", choices=[2, 3], default=2, type=int)
+    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,
+            order=args.order,
+            use_quad=args.use_quad,
+            visualize=args.visualize)

examples/advection/var-velocity.py

-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+__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
+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 os
 
 import numpy as np
-import pyopencl as cl  # noqa
-import pyopencl.array  # noqa
-import pyopencl.clmath  # noqa
 
-import pytest  # noqa
+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
+
+import grudge.dof_desc as dof_desc
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext
 
-from pyopencl.tools import (  # noqa
-                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
 
-import logging
 logger = logging.getLogger(__name__)
 
-from grudge import sym, bind, Discretization
-from pytools.obj_array import join_fields
 
+# {{{ plotting (keep in sync with `weak.py`)
+
+class Plotter:
+    def __init__(self, actx, dcoll, order, visualize=True, ylim=None):
+        self.actx = actx
+        self.dim = dcoll.ambient_dim
+
+        self.visualize = visualize
+        if not self.visualize:
+            return
+
+        if self.dim == 1:
+            import matplotlib.pyplot as pt
+            self.fig = pt.figure(figsize=(8, 8), dpi=300)
+            self.ylim = ylim
+
+            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(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, self.actx))
+
+            filename = f"{basename}.png"
+            if not overwrite and os.path.exists(filename):
+                from meshmode import FileExistsError
+                raise FileExistsError(f"output file '{filename}' already exists")
+
+            ax = self.fig.gca()
+            ax.plot(self.x, u, "-")
+            ax.plot(self.x, u, "k.")
+            if self.ylim is not None:
+                ax.set_ylim(self.ylim)
+
+            ax.set_xlabel("$x$")
+            ax.set_ylabel("$u$")
+            ax.set_title(f"t = {evt.t:.2f}")
+            self.fig.savefig(filename)
+            self.fig.clf()
+        else:
+            self.vis.write_vtk_file(f"{basename}.vtu", [
+                ("u", evt.state_component)
+                ], overwrite=overwrite)
+
+# }}}
 
-def main(write_output=True, order=4):
-    cl_ctx = cl.create_some_context()
+
+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)
 
-    dim = 2
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
+
+    # {{{ parameters
+
+    # domain [0, d]^dim
+    d = 1.0
+    # number of points in each dimension
+    npoints = 25
+
+    # final time
+    final_time = 1
+
+    if use_quad:
+        qtag = dof_desc.DISCR_TAG_QUAD
+    else:
+        qtag = None
+
+    # }}}
+
+    # {{{ discretization
 
     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=(15, 15), order=order)
+    mesh = generate_regular_rect_mesh(
+            a=(0,)*dim, b=(d,)*dim,
+            npoints_per_axis=(npoints,)*dim,
+            order=order)
 
-    dt_factor = 4
-    h = 1/20
+    from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
 
-    discr = Discretization(cl_ctx, mesh, order=order)
+    if use_quad:
+        discr_tag_to_group_factory = {
+            qtag: QuadratureSimplexGroupFactory(order=4*order)
+        }
+    else:
+        discr_tag_to_group_factory = {}
 
-    sym_x = sym.nodes(2)
+    from grudge.discretization import make_discretization_collection
 
-    advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2
+    dcoll = make_discretization_collection(
+        actx, mesh, order=order,
+        discr_tag_to_group_factory=discr_tag_to_group_factory
+    )
 
-    flux_type = "central"
+    # }}}
 
-    source_center = np.array([0.1, 0.1])
-    source_width = 0.05
+    # {{{ advection operator
 
-    sym_x = sym.nodes(2)
-    sym_source_center_dist = sym_x - source_center
+    # gaussian parameters
 
-    def f(x):
-        return sym.exp(
-                    -np.dot(sym_source_center_dist, sym_source_center_dist)
-                    / source_width**2)
+    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_analytic(x):
-        return 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
 
-    discr = Discretization(cl_ctx, mesh, order=order)
-    op = VariableCoefficientAdvectionOperator(2, advec_v,
-        inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
-        flux_type=flux_type)
+    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)
 
-    bound_op = bind(discr, op.sym_operator())
+    def rhs(t, u):
+        return adv_operator.operator(t, u)
 
-    u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
+    dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u))
 
-    def rhs(t, u):
-        return bound_op(queue, t=t, u=u)
+    logger.info("Timestep size: %g", dt)
 
-    final_time = 2
+    # }}}
 
-    dt = dt_factor * h/order**2
-    nsteps = (final_time // dt) + 1
-    dt = final_time/nsteps + 1e-15
+    # {{{ time stepping
 
     from grudge.shortcuts import set_up_rk4
-    dt_stepper = set_up_rk4("u", dt, u, rhs)
-
-    from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
+    dt_stepper = set_up_rk4("u", float(dt), u, rhs)
+    plot = Plotter(actx, dcoll, order, visualize=visualize,
+            ylim=[-0.1, 1.1])
 
     step = 0
-
     for event in dt_stepper.run(t_end=final_time):
-        if isinstance(event, dt_stepper.StateComputed):
+        if not isinstance(event, dt_stepper.StateComputed):
+            continue
+
+        if step % 10 == 0:
+            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
+        step += 1
 
-            #print(step, event.t, norm(queue, u=event.state_component[0]))
-            vis.write_vtk_file("fld-%04d.vtu" % step,
-                    [("u", event.state_component)])
+    # }}}
 
 
 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("--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,
+         order=args.order,
+         use_quad=args.use_quad,
+         visualize=args.visualize,
+         flux_type=args.flux)

examples/advection/weak.py

-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+__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
+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 os
 
 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.tools as cl_tools
+from arraycontext import flatten
+from meshmode.mesh import BTAG_ALL
+
+import grudge.dof_desc as dof_desc
+import grudge.op as op
+from grudge.array_context import PyOpenCLArrayContext
 
-from pyopencl.tools import (  # noqa
-                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
 
-import logging
 logger = logging.getLogger(__name__)
 
-from grudge import sym, bind, Discretization
 
-import numpy.linalg as la
+# {{{ plotting (keep in sync with `var-velocity.py`)
 
+class Plotter:
+    def __init__(self, actx, dcoll, order, visualize=True, ylim=None):
+        self.actx = actx
+        self.dim = dcoll.ambient_dim
 
+        self.visualize = visualize
+        if not self.visualize:
+            return
 
+        if self.dim == 1:
+            import matplotlib.pyplot as pt
+            self.fig = pt.figure(figsize=(8, 8), dpi=300)
+            self.ylim = ylim
 
-def main(write_output=True, order=4):
-    cl_ctx = cl.create_some_context()
-    queue = cl.CommandQueue(cl_ctx)
+            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(dcoll)
 
-    dim = 2
+    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, self.actx))
 
-    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)
+            filename = f"{basename}.png"
+            if not overwrite and os.path.exists(filename):
+                from meshmode import FileExistsError
+                raise FileExistsError(f"output file '{filename}' already exists")
 
-    dt_factor = 4
-    h = 1/20
+            ax = self.fig.gca()
+            ax.plot(self.x, u, "-")
+            ax.plot(self.x, u, "k.")
+            if self.ylim is not None:
+                ax.set_ylim(self.ylim)
 
-    discr = Discretization(cl_ctx, mesh, order=order)
+            ax.set_xlabel("$x$")
+            ax.set_ylabel("$u$")
+            ax.set_title(f"t = {evt.t:.2f}")
 
-    c = np.array([0.1,0.1])
-    norm_c = la.norm(c)
+            self.fig.savefig(filename)
+            self.fig.clf()
+        else:
+            self.vis.write_vtk_file(f"{basename}.vtu", [
+                ("u", evt.state_component)
+                ], overwrite=overwrite)
 
+# }}}
 
-    flux_type = "central"
-         
 
-    def f(x):
-        return sym.sin(3*x)
+def main(ctx_factory, dim=2, order=4, visualize=False):
+    cl_ctx = ctx_factory()
+    queue = cl.CommandQueue(cl_ctx)
 
-    def u_analytic(x):
-        return f(-np.dot(c, x)/norm_c+sym.var("t", sym.DD_SCALAR)*norm_c)
+    allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
+    actx = PyOpenCLArrayContext(queue, allocator=allocator)
 
-    from grudge.models.advection import WeakAdvectionOperator
-    from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    
-    discr = Discretization(cl_ctx, mesh, order=order)
-    op = WeakAdvectionOperator(c,
-        inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
-        flux_type=flux_type)
+    # {{{ parameters
 
-    bound_op = bind(discr, op.sym_operator())
+    # domain [-d/2, d/2]^dim
+    d = 1.0
+    # number of points in each dimension
+    npoints = 20
 
-    u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0)
+    # final time
+    final_time = 1.0
 
-    def rhs(t, u):
-        return bound_op(queue, t=t, u=u)
+    # velocity field
+    c = np.array([0.5] * dim)
+    norm_c = la.norm(c)
 
-    final_time = 0.3
+    # flux
+    flux_type = "central"
 
-    dt = dt_factor * h/order**2
-    nsteps = (final_time // dt) + 1
-    dt = final_time/nsteps + 1e-15
+    # }}}
 
+    # {{{ discretization
 
-    from grudge.shortcuts import set_up_rk4
-    dt_stepper = set_up_rk4("u", dt, u, rhs)
+    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)
 
-    last_u = None
+    from grudge.discretization import make_discretization_collection
 
-    from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
+    dcoll = make_discretization_collection(actx, mesh, order=order)
 
-    step = 0
+    # }}}
 
-    norm = bind(discr, sym.norm(2, sym.var("u")))
+    # {{{ weak advection operator
 
-    for event in dt_stepper.run(t_end=final_time):
-        if isinstance(event, dt_stepper.StateComputed):
+    def f(x):
+        return actx.np.sin(3 * x)
 
-            step += 1
+    def u_analytic(x, t=0):
+        return f(-np.dot(c, x) / norm_c + t * norm_c)
 
-            #print(step, event.t, norm(queue, u=event.state_component[0]))
-            vis.write_vtk_file("fld-%04d.vtu" % step,
-                    [  ("u", event.state_component) ])
+    from grudge.models.advection import WeakAdvectionOperator
 
+    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 adv_operator.operator(t, u)
 
+    dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u))
 
+    logger.info("Timestep size: %g", dt)
 
+    # }}}
 
-if __name__ == "__main__":
-    main()
+    # {{{ time stepping
 
+    from grudge.shortcuts import set_up_rk4
+    dt_stepper = set_up_rk4("u", float(dt), u, rhs)
+    plot = Plotter(actx, dcoll, order, visualize=visualize,
+            ylim=[-1.1, 1.1])
 
+    step = 0
+    norm_u = 0.0
+    for event in dt_stepper.run(t_end=final_time):
+        if not isinstance(event, dt_stepper.StateComputed):
+            continue
+
+        if step % 10 == 0:
+            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
+
+    # }}}
+
+
+if __name__ == "__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/burgers/burgers.py

-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-import numpy.linalg as la
-from math import sin, cos, pi, sqrt
-from pytools.test import mark_test
-
-
-
-
-class ExactTestCase:
-    a = 0
-    b = 150
-    final_time = 1000
-
-    def u0(self, x):
-        return self.u_exact(x, 0)
-
-    def u_exact(self, x, t):
-        # CAUTION: This gets the shock speed wrong as soon as the pulse
-        # starts interacting with itself.
-
-        def f(x, shock_loc):
-            if x < (t-40)/4:
-                return 1/4
-            else:
-                if t < 40:
-                    if x < (3*t)/4:
-                        return (x+15)/(t+20)
-                    elif x < (t+80)/4:
-                        return (x-30)/(t-40)
-                    else:
-                        return 1/4
-                else:
-                    if x < shock_loc:
-                        return (x+15)/(t+20)
-                    else:
-                        return 1/4
-
-        from math import sqrt
-
-        shock_loc = 30*sqrt(2*t+40)/sqrt(120) + t/4 - 10
-        shock_win = (shock_loc + 20) // self.b
-        x += shock_win * 150
-
-        x -= 20
-
-        return max(f(x, shock_loc), f(x-self.b, shock_loc-self.b))
-
-class OffCenterMigratingTestCase:
-    a = -pi
-    b = pi
-    final_time = 10
-
-    def u0(self, x):
-        return -0.4+sin(x+0.1)
-
-
-class CenteredStationaryTestCase:
-    # does funny things to P-P
-    a = -pi
-    b = pi
-    final_time = 10
-
-    def u0(self, x):
-        return -sin(x)
-
-class OffCenterStationaryTestCase:
-    # does funny things to P-P
-    a = -pi
-    b = pi
-    final_time = 10
-
-    def u0(self, x):
-        return -sin(x+0.3)
-
-
-
-def main(write_output=True, flux_type_arg="upwind", 
-        #case = CenteredStationaryTestCase(),
-        #case = OffCenterStationaryTestCase(),
-        #case = OffCenterMigratingTestCase(),
-        case = ExactTestCase(),
-        ):
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context()
-
-    order = 3
-    if rcon.is_head_rank:
-        if True:
-            from grudge.mesh.generator import make_uniform_1d_mesh
-            mesh = make_uniform_1d_mesh(case.a, case.b, 20, periodic=True)
-        else:
-            from grudge.mesh.generator import make_rect_mesh
-            print((pi*2)/(11*5*2))
-            mesh = make_rect_mesh((-pi, -1), (pi, 1),
-                    periodicity=(True, True),
-                    subdivisions=(11,5),
-                    max_area=(pi*2)/(11*5*2)
-                    )
-
-    if rcon.is_head_rank:
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    discr = rcon.make_discretization(mesh_data, order=order,
-            quad_min_degrees={"quad": 3*order})
-
-    if write_output:
-        from grudge.visualization import VtkVisualizer
-        vis = VtkVisualizer(discr, rcon, "fld")
-
-    # operator setup ----------------------------------------------------------
-    from grudge.second_order import IPDGSecondDerivative
-
-    from grudge.models.burgers import BurgersOperator
-    op = BurgersOperator(mesh.dimensions,
-            viscosity_scheme=IPDGSecondDerivative())
-
-    if rcon.is_head_rank:
-        print("%d elements" % len(discr.mesh.elements))
-
-    # exact solution ----------------------------------------------------------
-    import pymbolic
-    var = pymbolic.var
-
-    u = discr.interpolate_volume_function(lambda x, el: case.u0(x[0]))
-
-    # diagnostics setup -------------------------------------------------------
-    from pytools.log import LogManager, \
-            add_general_quantities, \
-            add_simulation_quantities, \
-            add_run_info
-
-    if write_output:
-        log_file_name = "burgers.dat"
-    else:
-        log_file_name = None
-
-    logmgr = LogManager(log_file_name, "w", rcon.communicator)
-    add_run_info(logmgr)
-    add_general_quantities(logmgr)
-    add_simulation_quantities(logmgr)
-    discr.add_instrumentation(logmgr)
-
-    from grudge.log import LpNorm
-    u_getter = lambda: u
-    logmgr.add_quantity(LpNorm(u_getter, discr, p=1, name="l1_u"))
-
-    logmgr.add_watches(["step.max", "t_sim.max", "l1_u", "t_step.max"])
-
-    # timestep loop -----------------------------------------------------------
-    rhs = op.bind(discr)
-
-    from grudge.timestep.runge_kutta import ODE45TimeStepper, LSRK4TimeStepper
-    stepper = ODE45TimeStepper()
-
-    stepper.add_instrumentation(logmgr)
-
-    try:
-        from grudge.timestep import times_and_steps
-        # for visc=0.01
-        #stab_fac = 0.1 # RK4
-        #stab_fac = 1.6 # dumka3(3), central
-        #stab_fac = 3 # dumka3(4), central
-
-        #stab_fac = 0.01 # RK4
-        stab_fac = 0.2 # dumka3(3), central
-        #stab_fac = 3 # dumka3(4), central
-
-        dt = stab_fac*op.estimate_timestep(discr,
-                stepper=LSRK4TimeStepper(), t=0, fields=u)
-
-        step_it = times_and_steps(
-                final_time=case.final_time, logmgr=logmgr, max_dt_getter=lambda t: dt)
-        from grudge.symbolic import  InverseVandermondeOperator
-        inv_vdm = InverseVandermondeOperator().bind(discr)
-
-        for step, t, dt in step_it:
-            if step % 3 == 0 and write_output:
-                if hasattr(case, "u_exact"):
-                    extra_fields = [
-                            ("u_exact",
-                                discr.interpolate_volume_function(
-                                    lambda x, el: case.u_exact(x[0], t)))]
-                else:
-                    extra_fields = []
-
-                visf = vis.make_file("fld-%04d" % step)
-                vis.add_data(visf, [
-                    ("u", u),
-                    ] + extra_fields,
-                    time=t,
-                    step=step)
-                visf.close()
-
-            u = stepper(u, t, dt, rhs)
-
-        if isinstance(case, ExactTestCase):
-            assert discr.norm(u, 1) < 50
-
-    finally:
-        if write_output:
-            vis.close()
-
-        logmgr.save()
-
-
-
-
-if __name__ == "__main__":
-    main()
-
-
-
-
-# entry points for py.test ----------------------------------------------------
-@mark_test.long
-def test_stability():
-    main(write_output=False)
-

examples/conftest.py

-# This file tells py.test to scan for test_xxx() functions in
-# any file below here, not just those named test_xxxx.py.
-
-
-def pytest_collect_file(path, parent):
-    if "grudge/examples" in str(path.dirpath()) and path.ext == ".py":
-        return parent.Module(path, parent=parent)

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/gas_dynamics/box-in-box.py

-# Copyright (C) 2008 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-
-
-def make_boxmesh():
-    from meshpy.tet import MeshInfo, build
-    from meshpy.geometry import GeometryBuilder, Marker, make_box
-
-    geob = GeometryBuilder()
-
-    box_marker = Marker.FIRST_USER_MARKER
-    extent_small = 0.1*numpy.ones(3, dtype=numpy.float64)
-
-    geob.add_geometry(*make_box(-extent_small, extent_small))
-
-    # make small "separator box" for region attribute
-
-    geob.add_geometry(
-            *make_box(
-                -extent_small*4,
-                numpy.array([4, 0.4, 0.4], dtype=numpy.float64)))
-
-    geob.add_geometry(
-            *make_box(
-                numpy.array([-1, -1, -1], dtype=numpy.float64),
-                numpy.array([5, 1, 1], dtype=numpy.float64)))
-
-    mesh_info = MeshInfo()
-    geob.set(mesh_info)
-    mesh_info.set_holes([(0, 0, 0)])
-
-    # region attributes
-    mesh_info.regions.resize(1)
-    mesh_info.regions[0] = (
-            # point in region
-            list(extent_small*2) + [
-            # region number
-            1,
-            # max volume in region
-            #0.0001
-            0.005
-            ])
-
-    mesh = build(mesh_info, max_volume=0.02,
-            volume_constraints=True, attributes=True)
-    print("%d elements" % len(mesh.elements))
-    #mesh.write_vtk("box-in-box.vtk")
-    #print "done writing"
-
-    fvi2fm = mesh.face_vertex_indices_to_face_marker
-
-    face_marker_to_tag = {
-            box_marker: "noslip",
-            Marker.MINUS_X: "inflow",
-            Marker.PLUS_X: "outflow",
-            Marker.MINUS_Y: "inflow",
-            Marker.PLUS_Y: "inflow",
-            Marker.PLUS_Z: "inflow",
-            Marker.MINUS_Z: "inflow"
-            }
-
-    def bdry_tagger(fvi, el, fn, all_v):
-        face_marker = fvi2fm[fvi]
-        return [face_marker_to_tag[face_marker]]
-
-    from grudge.mesh import make_conformal_mesh
-    return make_conformal_mesh(
-            mesh.points, mesh.elements, bdry_tagger)
-
-
-def main():
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context(["cuda"])
-
-    if rcon.is_head_rank:
-        mesh = make_boxmesh()
-        #from grudge.mesh import make_rect_mesh
-        #mesh = make_rect_mesh(
-        #       boundary_tagger=lambda fvi, el, fn, all_v: ["inflow"])
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    for order in [3]:
-        from pytools import add_python_path_relative_to_script
-        add_python_path_relative_to_script("..")
-
-        from gas_dynamics_initials import UniformMachFlow
-        box = UniformMachFlow(angle_of_attack=0)
-
-        from grudge.models.gas_dynamics import GasDynamicsOperator
-        op = GasDynamicsOperator(dimensions=3,
-                gamma=box.gamma, mu=box.mu,
-                prandtl=box.prandtl, spec_gas_const=box.spec_gas_const,
-                bc_inflow=box, bc_outflow=box, bc_noslip=box,
-                inflow_tag="inflow", outflow_tag="outflow", noslip_tag="noslip")
-
-        discr = rcon.make_discretization(mesh_data, order=order,
-                        debug=[
-                            #"cuda_no_plan",
-                            #"cuda_dump_kernels",
-                            #"dump_dataflow_graph",
-                            #"dump_optemplate_stages",
-                            #"dump_dataflow_graph",
-                            #"print_op_code",
-                            "cuda_no_plan_el_local",
-                            ],
-                        default_scalar_type=numpy.float32,
-                        tune_for=op.sym_operator())
-
-        from grudge.visualization import SiloVisualizer, VtkVisualizer  # noqa
-        #vis = VtkVisualizer(discr, rcon, "shearflow-%d" % order)
-        vis = SiloVisualizer(discr, rcon)
-
-        fields = box.volume_interpolant(0, discr)
-
-        navierstokes_ex = op.bind(discr)
-
-        max_eigval = [0]
-
-        def rhs(t, q):
-            ode_rhs, speed = navierstokes_ex(t, q)
-            max_eigval[0] = speed
-            return ode_rhs
-
-        rhs(0, fields)
-
-        if rcon.is_head_rank:
-            print("---------------------------------------------")
-            print("order %d" % order)
-            print("---------------------------------------------")
-            print("#elements=", len(mesh.elements))
-
-        from grudge.timestep import RK4TimeStepper
-        stepper = RK4TimeStepper()
-
-        # diagnostics setup ---------------------------------------------------
-        from pytools.log import LogManager, add_general_quantities, \
-                add_simulation_quantities, add_run_info
-
-        logmgr = LogManager("navierstokes-%d.dat" % order, "w", rcon.communicator)
-        add_run_info(logmgr)
-        add_general_quantities(logmgr)
-        add_simulation_quantities(logmgr)
-        discr.add_instrumentation(logmgr)
-        stepper.add_instrumentation(logmgr)
-
-        logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
-
-        from pytools.log import LogQuantity
-
-        class ChangeSinceLastStep(LogQuantity):
-            """Records the change of a variable between a time step and the previous
-               one"""
-
-            def __init__(self, name="change"):
-                LogQuantity.__init__(self, name, "1", "Change since last time step")
-
-                self.old_fields = 0
-
-            def __call__(self):
-                result = discr.norm(fields - self.old_fields)
-                self.old_fields = fields
-                return result
-
-        logmgr.add_quantity(ChangeSinceLastStep())
-
-        # timestep loop -------------------------------------------------------
-        try:
-            from grudge.timestep import times_and_steps
-            step_it = times_and_steps(
-                    final_time=200,
-                    #max_steps=500,
-                    logmgr=logmgr,
-                    max_dt_getter=lambda t: op.estimate_timestep(discr,
-                        stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
-
-            for step, t, dt in step_it:
-                if step % 200 == 0:
-                #if False:
-                    visf = vis.make_file("box-%d-%06d" % (order, step))
-
-                    #rhs_fields = rhs(t, fields)
-
-                    vis.add_data(visf,
-                            [
-                                ("rho", discr.convert_volume(
-                                    op.rho(fields), kind="numpy")),
-                                ("e", discr.convert_volume(
-                                    op.e(fields), kind="numpy")),
-                                ("rho_u", discr.convert_volume(
-                                    op.rho_u(fields), kind="numpy")),
-                                ("u", discr.convert_volume(
-                                    op.u(fields), kind="numpy")),
-
-                                # ("rhs_rho", discr.convert_volume(
-                                #     op.rho(rhs_fields), kind="numpy")),
-                                # ("rhs_e", discr.convert_volume(
-                                #     op.e(rhs_fields), kind="numpy")),
-                                # ("rhs_rho_u", discr.convert_volume(
-                                #     op.rho_u(rhs_fields), kind="numpy")),
-                                ],
-                            expressions=[
-                                ("p", "(0.4)*(e- 0.5*(rho_u*u))"),
-                                ],
-                            time=t, step=step
-                            )
-                    visf.close()
-
-                fields = stepper(fields, t, dt, rhs)
-
-        finally:
-            vis.close()
-            logmgr.save()
-            discr.close()
-
-if __name__ == "__main__":
-    main()

examples/gas_dynamics/euler/sine-wave.py

-# Copyright (C) 2008 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-import numpy.linalg as la
-
-
-
-
-class SineWave:
-    def __init__(self):
-        self.gamma = 1.4
-        self.mu = 0
-        self.prandtl = 0.72
-        self.spec_gas_const = 287.1
-
-    def __call__(self, t, x_vec):
-        rho = 2 + numpy.sin(x_vec[0] + x_vec[1] + x_vec[2] - 2 * t)
-        velocity = numpy.array([1, 1, 0])
-        p = 1
-        e = p/(self.gamma-1) + rho/2 * numpy.dot(velocity, velocity)
-        rho_u = rho * velocity[0]
-        rho_v = rho * velocity[1]
-        rho_w = rho * velocity[2]
-
-        from grudge.tools import join_fields
-        return join_fields(rho, e, rho_u, rho_v, rho_w)
-
-    def properties(self):
-        return(self.gamma, self.mu, self.prandtl, self.spec_gas_const)
-
-    def volume_interpolant(self, t, discr):
-        return discr.convert_volume(
-                        self(t, discr.nodes.T),
-                        kind=discr.compute_kind)
-
-    def boundary_interpolant(self, t, discr, tag):
-        return discr.convert_boundary(
-                        self(t, discr.get_boundary(tag).nodes.T),
-                         tag=tag, kind=discr.compute_kind)
-
-
-
-
-def main(final_time=1, write_output=False):
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context()
-
-    from grudge.tools import EOCRecorder, to_obj_array
-    eoc_rec = EOCRecorder()
-
-    if rcon.is_head_rank:
-        from grudge.mesh import make_box_mesh
-        mesh = make_box_mesh((0,0,0), (10,10,10), max_volume=0.5)
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    for order in [3, 4, 5]:
-        discr = rcon.make_discretization(mesh_data, order=order,
-                        default_scalar_type=numpy.float64)
-
-        from grudge.visualization import SiloVisualizer, VtkVisualizer
-        vis = VtkVisualizer(discr, rcon, "sinewave-%d" % order)
-        #vis = SiloVisualizer(discr, rcon)
-
-        sinewave = SineWave()
-        fields = sinewave.volume_interpolant(0, discr)
-        gamma, mu, prandtl, spec_gas_const = sinewave.properties()
-
-        from grudge.mesh import BTAG_ALL
-        from grudge.models.gas_dynamics import GasDynamicsOperator
-        op = GasDynamicsOperator(dimensions=mesh.dimensions, gamma=gamma, mu=mu,
-                prandtl=prandtl, spec_gas_const=spec_gas_const,
-                bc_inflow=sinewave, bc_outflow=sinewave, bc_noslip=sinewave,
-                inflow_tag=BTAG_ALL, source=None)
-
-        euler_ex = op.bind(discr)
-
-        max_eigval = [0]
-        def rhs(t, q):
-            ode_rhs, speed = euler_ex(t, q)
-            max_eigval[0] = speed
-            return ode_rhs
-        rhs(0, fields)
-
-        if rcon.is_head_rank:
-            print("---------------------------------------------")
-            print("order %d" % order)
-            print("---------------------------------------------")
-            print("#elements=", len(mesh.elements))
-
-        from grudge.timestep import RK4TimeStepper
-        stepper = RK4TimeStepper()
-
-        # diagnostics setup ---------------------------------------------------
-        from pytools.log import LogManager, add_general_quantities, \
-                add_simulation_quantities, add_run_info
-
-        if write_output:
-            log_name = ("euler-sinewave-%(order)d-%(els)d.dat"
-                    % {"order":order, "els":len(mesh.elements)})
-        else:
-            log_name = False
-        logmgr = LogManager(log_name, "w", rcon.communicator)
-        add_run_info(logmgr)
-        add_general_quantities(logmgr)
-        add_simulation_quantities(logmgr)
-        discr.add_instrumentation(logmgr)
-        stepper.add_instrumentation(logmgr)
-
-        logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
-
-        # timestep loop -------------------------------------------------------
-        try:
-            from grudge.timestep import times_and_steps
-            step_it = times_and_steps(
-                    final_time=final_time, logmgr=logmgr,
-                    max_dt_getter=lambda t: op.estimate_timestep(discr,
-                        stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
-
-            for step, t, dt in step_it:
-                #if step % 10 == 0:
-                if write_output:
-                    visf = vis.make_file("sinewave-%d-%04d" % (order, step))
-
-                    #from pyvisfile.silo import DB_VARTYPE_VECTOR
-                    vis.add_data(visf,
-                            [
-                                ("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
-                                ("e", discr.convert_volume(op.e(fields), kind="numpy")),
-                                ("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
-                                ("u", discr.convert_volume(op.u(fields), kind="numpy")),
-
-                                #("true_rho", op.rho(true_fields)),
-                                #("true_e", op.e(true_fields)),
-                                #("true_rho_u", op.rho_u(true_fields)),
-                                #("true_u", op.u(true_fields)),
-
-                                #("rhs_rho", op.rho(rhs_fields)),
-                                #("rhs_e", op.e(rhs_fields)),
-                                #("rhs_rho_u", op.rho_u(rhs_fields)),
-                                ],
-                            #expressions=[
-                                #("diff_rho", "rho-true_rho"),
-                                #("diff_e", "e-true_e"),
-                                #("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
-
-                                #("p", "0.4*(e- 0.5*(rho_u*u))"),
-                                #],
-                            time=t, step=step
-                            )
-                    visf.close()
-
-                fields = stepper(fields, t, dt, rhs)
-
-        finally:
-            vis.close()
-            logmgr.close()
-            discr.close()
-
-        true_fields = sinewave.volume_interpolant(t, discr)
-        eoc_rec.add_data_point(order, discr.norm(fields-true_fields))
-        print()
-        print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
-
-
-
-
-if __name__ == "__main__":
-    main()
-
-
-
-
-# entry points for py.test ----------------------------------------------------
-from pytools.test import mark_test
-@mark_test.long
-def test_euler_sine_wave():
-    main(final_time=0.1, write_output=False)
-

examples/gas_dynamics/euler/sod-2d.py

-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-import numpy.linalg as la
-
-
-
-
-class Sod:
-    def __init__(self, gamma):
-        self.gamma = gamma
-        self.prandtl = 0.72
-
-    def __call__(self, t, x_vec):
-
-        from grudge.tools import heaviside
-        from grudge.tools import heaviside_a
-
-        x_rel = x_vec[0]
-        y_rel = x_vec[1]
-
-        from math import pi
-        r = numpy.sqrt(x_rel**2+y_rel**2)
-        r_shift=r-3.0
-        u = 0.0
-        v = 0.0
-        from numpy import sign
-        rho = heaviside(-r_shift)+.125*heaviside_a(r_shift,1.0)
-        e = (1.0/(self.gamma-1.0))*(heaviside(-r_shift)+.1*heaviside_a(r_shift,1.0))
-        p = (self.gamma-1.0)*e
-
-        from grudge.tools import join_fields
-        return join_fields(rho, e, rho*u, rho*v)
-
-
-    def volume_interpolant(self, t, discr):
-        return discr.convert_volume(
-                        self(t, discr.nodes.T),
-                        kind=discr.compute_kind)
-
-    def boundary_interpolant(self, t, discr, tag):
-        return discr.convert_boundary(
-                        self(t, discr.get_boundary(tag).nodes.T),
-                         tag=tag, kind=discr.compute_kind)
-
-
-
-
-def main():
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context()
-
-    from grudge.tools import to_obj_array
-
-    if rcon.is_head_rank:
-        from grudge.mesh.generator import make_rect_mesh
-        mesh = make_rect_mesh((-5,-5), (5,5), max_area=0.01)
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    for order in [1]:
-        discr = rcon.make_discretization(mesh_data, order=order,
-                        default_scalar_type=numpy.float64)
-
-        from grudge.visualization import SiloVisualizer, VtkVisualizer
-        vis = VtkVisualizer(discr, rcon, "Sod2D-%d" % order)
-        #vis = SiloVisualizer(discr, rcon)
-
-        sod_field = Sod(gamma=1.4)
-        fields = sod_field.volume_interpolant(0, discr)
-
-        from grudge.models.gas_dynamics import GasDynamicsOperator
-        from grudge.mesh import BTAG_ALL
-        op = GasDynamicsOperator(dimensions=2, gamma=sod_field.gamma, mu=0.0,
-                prandtl=sod_field.prandtl,
-                bc_inflow=sod_field,
-                bc_outflow=sod_field,
-                bc_noslip=sod_field,
-                inflow_tag=BTAG_ALL,
-                source=None)
-
-        euler_ex = op.bind(discr)
-
-        max_eigval = [0]
-        def rhs(t, q):
-            ode_rhs, speed = euler_ex(t, q)
-            max_eigval[0] = speed
-            return ode_rhs
-        rhs(0, fields)
-
-        if rcon.is_head_rank:
-            print("---------------------------------------------")
-            print("order %d" % order)
-            print("---------------------------------------------")
-            print("#elements=", len(mesh.elements))
-
-        # limiter setup ------------------------------------------------------------
-        from grudge.models.gas_dynamics import SlopeLimiter1NEuler
-        limiter =  SlopeLimiter1NEuler(discr, sod_field.gamma, 2, op)
-
-        # integrator setup---------------------------------------------------------
-        from grudge.timestep import SSPRK3TimeStepper, RK4TimeStepper
-        stepper = SSPRK3TimeStepper(limiter=limiter)
-        #stepper = SSPRK3TimeStepper()
-        #stepper = RK4TimeStepper()
-
-        # diagnostics setup ---------------------------------------------------
-        from pytools.log import LogManager, add_general_quantities, \
-                add_simulation_quantities, add_run_info
-
-        logmgr = LogManager("euler-%d.dat" % order, "w", rcon.communicator)
-        add_run_info(logmgr)
-        add_general_quantities(logmgr)
-        add_simulation_quantities(logmgr)
-        discr.add_instrumentation(logmgr)
-        stepper.add_instrumentation(logmgr)
-
-        logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
-
-        # filter setup-------------------------------------------------------------
-        from grudge.discretization import Filter, ExponentialFilterResponseFunction
-        mode_filter = Filter(discr,
-                ExponentialFilterResponseFunction(min_amplification=0.9,order=4))
-
-        # timestep loop -------------------------------------------------------
-        try:
-            from grudge.timestep import times_and_steps
-            step_it = times_and_steps(
-                    final_time=1.0, logmgr=logmgr,
-                    max_dt_getter=lambda t: op.estimate_timestep(discr,
-                        stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
-
-            for step, t, dt in step_it:
-                if step % 5 == 0:
-                #if False:
-                    visf = vis.make_file("vortex-%d-%04d" % (order, step))
-
-                    #true_fields = vortex.volume_interpolant(t, discr)
-
-                    #from pyvisfile.silo import DB_VARTYPE_VECTOR
-                    vis.add_data(visf,
-                            [
-                                ("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
-                                ("e", discr.convert_volume(op.e(fields), kind="numpy")),
-                                ("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
-                                ("u", discr.convert_volume(op.u(fields), kind="numpy")),
-
-                                #("true_rho", op.rho(true_fields)),
-                                #("true_e", op.e(true_fields)),
-                                #("true_rho_u", op.rho_u(true_fields)),
-                                #("true_u", op.u(true_fields)),
-
-                                #("rhs_rho", op.rho(rhs_fields)),
-                                #("rhs_e", op.e(rhs_fields)),
-                                #("rhs_rho_u", op.rho_u(rhs_fields)),
-                                ],
-                            #expressions=[
-                                #("diff_rho", "rho-true_rho"),
-                                #("diff_e", "e-true_e"),
-                                #("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
-
-                                #("p", "0.4*(e- 0.5*(rho_u*u))"),
-                                #],
-                            time=t, step=step
-                            )
-                    visf.close()
-
-                fields = stepper(fields, t, dt, rhs)
-                # fields = limiter(fields)
-                # fields = mode_filter(fields)
-
-                assert not numpy.isnan(numpy.sum(fields[0]))
-        finally:
-            vis.close()
-            logmgr.close()
-            discr.close()
-
-        # not solution, just to check against when making code changes
-        true_fields = sod_field.volume_interpolant(t, discr)
-        print(discr.norm(fields-true_fields))
-
-if __name__ == "__main__":
-    main()

examples/gas_dynamics/euler/vortex-adaptive-grid.py

-# Copyright (C) 2008 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-import numpy.linalg as la
-
-
-
-
-def main(write_output=True):
-    from pytools import add_python_path_relative_to_script
-    add_python_path_relative_to_script("..")
-
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context()
-
-    from grudge.tools import EOCRecorder
-    eoc_rec = EOCRecorder()
-
-
-    if rcon.is_head_rank:
-        from grudge.mesh.generator import \
-                make_rect_mesh, \
-                make_centered_regular_rect_mesh
-
-        refine = 4
-        mesh = make_centered_regular_rect_mesh((0,-5), (10,5), n=(9,9),
-                post_refine_factor=refine)
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    # a second mesh to regrid to
-    if rcon.is_head_rank:
-        from grudge.mesh.generator import \
-                make_rect_mesh, \
-                make_centered_regular_rect_mesh
-
-        refine = 4
-        mesh2 = make_centered_regular_rect_mesh((0,-5), (10,5), n=(8,8),
-                post_refine_factor=refine)
-        mesh_data2 = rcon.distribute_mesh(mesh2)
-    else:
-        mesh_data2 = rcon.receive_mesh()
-
-
-
-    for order in [3,4]:
-        discr = rcon.make_discretization(mesh_data, order=order,
-                        default_scalar_type=numpy.float64,
-                        quad_min_degrees={
-                            "gasdyn_vol": 3*order,
-                            "gasdyn_face": 3*order,
-                            })
-
-        discr2 = rcon.make_discretization(mesh_data2, order=order,
-                        default_scalar_type=numpy.float64,
-                        quad_min_degrees={
-                            "gasdyn_vol": 3*order,
-                            "gasdyn_face": 3*order,
-                            })
-
-
-        from grudge.visualization import SiloVisualizer, VtkVisualizer
-        vis = VtkVisualizer(discr, rcon, "vortex-%d" % order)
-        #vis = SiloVisualizer(discr, rcon)
-
-        from gas_dynamics_initials import Vortex
-        vortex = Vortex()
-        fields = vortex.volume_interpolant(0, discr)
-
-        from grudge.models.gas_dynamics import GasDynamicsOperator
-        from grudge.mesh import BTAG_ALL
-
-        op = GasDynamicsOperator(dimensions=2, gamma=vortex.gamma, mu=vortex.mu,
-                prandtl=vortex.prandtl, spec_gas_const=vortex.spec_gas_const,
-                bc_inflow=vortex, bc_outflow=vortex, bc_noslip=vortex,
-                inflow_tag=BTAG_ALL, source=None)
-
-        euler_ex = op.bind(discr)
-
-        max_eigval = [0]
-        def rhs(t, q):
-            ode_rhs, speed = euler_ex(t, q)
-            max_eigval[0] = speed
-            return ode_rhs
-        rhs(0, fields)
-
-
-        if rcon.is_head_rank:
-            print("---------------------------------------------")
-            print("order %d" % order)
-            print("---------------------------------------------")
-            print("#elements for mesh 1 =", len(mesh.elements))
-            print("#elements for mesh 2 =", len(mesh2.elements))
-
-
-        # limiter ------------------------------------------------------------
-        from grudge.models.gas_dynamics import SlopeLimiter1NEuler
-        limiter = SlopeLimiter1NEuler(discr, vortex.gamma, 2, op)
-
-        from grudge.timestep import SSPRK3TimeStepper
-        #stepper = SSPRK3TimeStepper(limiter=limiter)
-        stepper = SSPRK3TimeStepper()
-
-        #from grudge.timestep import RK4TimeStepper
-        #stepper = RK4TimeStepper()
-
-        # diagnostics setup ---------------------------------------------------
-        from pytools.log import LogManager, add_general_quantities, \
-                add_simulation_quantities, add_run_info
-
-        if write_output:
-            log_file_name = "euler-%d.dat" % order
-        else:
-            log_file_name = None
-
-        logmgr = LogManager(log_file_name, "w", rcon.communicator)
-        add_run_info(logmgr)
-        add_general_quantities(logmgr)
-        add_simulation_quantities(logmgr)
-        discr.add_instrumentation(logmgr)
-        stepper.add_instrumentation(logmgr)
-
-        logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
-
-        # timestep loop -------------------------------------------------------
-        try:
-            final_time = 0.2
-            from grudge.timestep import times_and_steps
-            step_it = times_and_steps(
-                    final_time=final_time, logmgr=logmgr,
-                    max_dt_getter=lambda t: op.estimate_timestep(discr,
-                        stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
-
-            for step, t, dt in step_it:
-                if step % 10 == 0 and write_output:
-                #if False:
-                    visf = vis.make_file("vortex-%d-%04d" % (order, step))
-
-                    #true_fields = vortex.volume_interpolant(t, discr)
-
-                    from pyvisfile.silo import DB_VARTYPE_VECTOR
-                    vis.add_data(visf,
-                            [
-                                ("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
-                                ("e", discr.convert_volume(op.e(fields), kind="numpy")),
-                                ("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
-                                ("u", discr.convert_volume(op.u(fields), kind="numpy")),
-
-                                #("true_rho", discr.convert_volume(op.rho(true_fields), kind="numpy")),
-                                #("true_e", discr.convert_volume(op.e(true_fields), kind="numpy")),
-                                #("true_rho_u", discr.convert_volume(op.rho_u(true_fields), kind="numpy")),
-                                #("true_u", discr.convert_volume(op.u(true_fields), kind="numpy")),
-
-                                #("rhs_rho", discr.convert_volume(op.rho(rhs_fields), kind="numpy")),
-                                #("rhs_e", discr.convert_volume(op.e(rhs_fields), kind="numpy")),
-                                #("rhs_rho_u", discr.convert_volume(op.rho_u(rhs_fields), kind="numpy")),
-                                ],
-                            #expressions=[
-                                #("diff_rho", "rho-true_rho"),
-                                #("diff_e", "e-true_e"),
-                                #("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
-
-                                #("p", "0.4*(e- 0.5*(rho_u*u))"),
-                                #],
-                            time=t, step=step
-                            )
-                    visf.close()
-
-                fields = stepper(fields, t, dt, rhs)
-                #fields = limiter(fields)
-
-                #regrid to discr2 at some arbitrary time
-                if step == 21:
-
-                    #get interpolated fields
-                    fields = discr.get_regrid_values(fields, discr2, dtype=None, use_btree=True, thresh=1e-8)
-                    #get new stepper (old one has reference to discr
-                    stepper = SSPRK3TimeStepper()
-                    #new bind
-                    euler_ex = op.bind(discr2)
-                    #new rhs
-                    max_eigval = [0]
-                    def rhs(t, q):
-                        ode_rhs, speed = euler_ex(t, q)
-                        max_eigval[0] = speed
-                        return ode_rhs
-                    rhs(t+dt, fields)
-                    #add logmanager
-                    #discr2.add_instrumentation(logmgr)
-                    #new step_it
-                    step_it = times_and_steps(
-                        final_time=final_time, logmgr=logmgr,
-                        max_dt_getter=lambda t: op.estimate_timestep(discr2,
-                            stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
-
-                    #new visualization
-                    vis.close()
-                    vis = VtkVisualizer(discr2, rcon, "vortexNewGrid-%d" % order)
-                    discr=discr2
-
-
-
-                assert not numpy.isnan(numpy.sum(fields[0]))
-
-            true_fields = vortex.volume_interpolant(final_time, discr)
-            l2_error = discr.norm(fields-true_fields)
-            l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields))
-            l2_error_e = discr.norm(op.e(fields)-op.e(true_fields))
-            l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields))
-            l2_error_u = discr.norm(op.u(fields)-op.u(true_fields))
-
-            eoc_rec.add_data_point(order, l2_error)
-            print()
-            print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
-
-            logmgr.set_constant("l2_error", l2_error)
-            logmgr.set_constant("l2_error_rho", l2_error_rho)
-            logmgr.set_constant("l2_error_e", l2_error_e)
-            logmgr.set_constant("l2_error_rhou", l2_error_rhou)
-            logmgr.set_constant("l2_error_u", l2_error_u)
-            logmgr.set_constant("refinement", refine)
-
-        finally:
-            if write_output:
-                vis.close()
-
-            logmgr.close()
-            discr.close()
-
-
-
-    # after order loop
-    # assert eoc_rec.estimate_order_of_convergence()[0,1] > 6
-
-
-
-
-if __name__ == "__main__":
-    main()