Y1

.github/workflows/ci.yaml

@@ -40,14 +40,13 @@ jobs:
           run: |
             . .ci-support/install.sh
 
-        - name: Pylint
+        - name: Run Pylint
           run: |
             MINIFORGE_INSTALL_DIR=.miniforge3
             . "$MINIFORGE_INSTALL_DIR/bin/activate" testing
             python -m pip install wheel
             python -m pip install pylint pyyaml
-            curl -L -O https://raw.githubusercontent.com/inducer/ci-support/master/run-pylint.py
-            python run-pylint.py --yaml-rcfile=.pylintrc.yml "$(basename $GITHUB_REPOSITORY)" test examples && echo "Pylint found no errors."
+            ./run-pylint.sh
 
     pydocstyle:
         runs-on: ubuntu-latest
@@ -57,7 +56,7 @@ jobs:
           uses: actions/setup-python@v1
           with:
             python-version: '3.x'
-        - name: Lint
+        - name: Run Pydocstyle
           run: |
             python3 -m venv myenv
             source myenv/bin/activate
@@ -87,7 +86,7 @@ jobs:
               export PYOPENCL_TEST=port:pthread
             fi
 
-            python -m pytest --durations=10 --tb=native --junitxml=pytest.xml --doctest-modules -rxsw .
+            python -m pytest --durations=10 --tb=native --junitxml=pytest.xml --doctest-modules -rxsw . ../doc/*.rst ../doc/*/*.rst
 
     examples:
         name: Examples ${{ matrix.os }}

.gitignore

@@ -30,6 +30,7 @@ run-debug-*
 *.vtk
 *.silo
 *.dat
+*.pkl
 
 .cache
 .pytest_cache
@@ -43,3 +44,6 @@ test/nodal-dg
 # Emacs backups
 .#*
 \#*
+
+# Linting temp files
+.run-pylint.py

doc/conf.py

@@ -78,7 +78,7 @@
     "https://numpy.org/doc/stable/": None,
     "https://documen.tician.de/pyopencl/": None,
     "https://documen.tician.de/modepy/": None,
-    "https://documen.tician.de/arraycontext/": None,
+    "arraycontext": ("https://documen.tician.de/arraycontext/", None),
     "https://documen.tician.de/meshmode/": None,
     "https://documen.tician.de/grudge/": None,
     "https://documen.tician.de/pytato/": None,
@@ -98,7 +98,7 @@
 
 nitpicky = True
 
-mathjax_config = {
+mathjax3_config = {
     "tex2jax": {
         "inlineMath": [["\\(", "\\)"]],
         "displayMath": [["\\[", "\\]"]],

doc/misc.rst

@@ -71,3 +71,7 @@ References
     `(DOI) <https://doi.org/10.1007/978-3-642-59721-3_14>`__
 .. [Ihme_2014] Yu Lv and Matthias Ihme (2014) Journal of Computationsl Physics 270 105 \
     `(DOI) <http://dx.doi.org/10.1016/j.jcp.2014.03.029>`__
+.. [Persson_2012] P. Persson and J. Peraire, AIAA 44 \
+    `(DOI) <https://doi.org/10.2514/6.2006-112>`__
+.. [Woodward_1984] Woodward and Colella, Journal of Computational Physics, 54 \
+    `(DOI) <https://doi.org/10.1006/jcph.1996.5572>`__

doc/operators/artificial_viscosity.rst

@@ -0,0 +1,4 @@
+Artificial Viscosity
+====================
+
+.. automodule:: mirgecom.artificial_viscosity

doc/operators/operators.rst

@@ -6,4 +6,5 @@ Operators
 
    wave-eq
    diffusion
+   artificial_viscosity
    gas-dynamics

doc/support/arrays.rst

@@ -0,0 +1,102 @@
+Arrays and Array Containers in :mod:`mirgecom`
+==============================================
+
+:mod:`mirgecom` is quite flexible in terms of what arrays it can use; in fact,
+no aspect of the code prescribes the use of a specific type of array.
+Instead, we use a :class:`arraycontext.ArrayContext` to determine the
+types of arrays on which the computation should take place. At its simplest,
+this can be a :class:`arraycontext.PyOpenCLArrayContext`, which uses
+:class:`pyopencl.array.Array`\ s that are eagerly evaluated. "Eager evaluation"
+means that when you type ``a + b`` with two arrays ``a`` and ``b``, then the
+result is computed as soon as the expression is evaluated.
+
+Lazy/deferred evaluation
+------------------------
+For efficiency, it is sometimes better to defer the evaluation of the
+expression, to have an opportunity to apply some optimizations. We call this
+*lazy* or *deferred evaluation*; it is a common trick in efficiency-minded
+software (e.g. Google's Tensorflow, or expression templates in C++).
+This is realized in :mod:`mirgecom` via an :class:`~arraycontext.ArrayContext`
+based on :mod:`pytato` (not merged to ``main``) at the time of this writing.
+
+Frozen and thawed arrays
+------------------------
+
+All array contexts have a notion of *thawed* and *frozen* arrays that is important
+to understand, see :ref:`arraycontext:freeze-thaw` for more details.
+
+You might be wondering whether a given function or method will return frozen or
+thawed data. In general, the documentation should state that, however since
+a fair bit of functionality predates these concepts, we are still catching up
+in terms of updating the documentation. (Help welcome!)
+
+In the meantime, these rules of thumb should cover most cases:
+
+* If you did not provide an array context to the function
+  (explicitly or implicitly via an input array), you will receive frozen data.
+* Any data that is cached/long-lived/"at rest" is going to be frozen.
+* If the data is the result of a `memoized` function, then it will likely
+  be frozen.
+
+To demonstrate the effect of this, first we need some setup:
+
+.. doctest::
+
+   >>> import pyopencl as cl
+   >>> from arraycontext import PyOpenCLArrayContext, thaw
+   >>> ctx = cl.create_some_context()
+   >>> queue = cl.CommandQueue(ctx)
+   >>> actx = PyOpenCLArrayContext(queue)
+   >>> from meshmode.mesh.generation import generate_regular_rect_mesh
+   >>> mesh = generate_regular_rect_mesh(a=(0, 0), b=(1, 1), nelements_per_axis=(10, 10))
+   >>> from grudge import DiscretizationCollection
+   >>> dcoll = DiscretizationCollection(actx, mesh, order=5)
+
+Most quantities that are maintained by the discretization will be frozen. For example,
+if one wanted to grab the nodes of the mesh or normals of a named surface in the mesh:
+
+.. doctest::
+
+   >>> from grudge.dof_desc import DOFDesc
+   >>> nodes = thaw(dcoll.nodes(), actx)
+   >>> dd = DOFDesc("all_faces")
+   >>> nhat = thaw(dcoll.normal(dd), actx)
+
+What can go wrong?  Attempts to operate on frozen data will yield errors similar to
+the following:
+
+.. doctest::
+
+   >>> dcoll.nodes() * 5
+   Traceback (most recent call last):
+    ...
+   ValueError: array containers with frozen arrays cannot be operated upon
+
+Fortunately, recovering from this is straightforward:
+
+.. doctest::
+
+   >>> nodes = thaw(dcoll.nodes(), actx)
+   >>> result = nodes * 5
+
+Array Containers
+----------------
+
+Arrays in :mod:`mirgecom` live in (somewhat) deeply nested data structures
+that are :class:`~arraycontext.ArrayContainer`\ s. Array containers typically
+support arithmetic and can be passed to most methods in
+:class:`~arraycontext.ArrayContext` that take arrays: they simply get applied to
+all arrays in the container. The same goes for many of the
+discretization-focused functions in :mod:`grudge` and :mod:`meshmode`.
+
+For example, think of the solver state for Euler's equations of gas dynamics:
+
+* At the outermost level, there is :class:`mirgecom.fluid.ConservedVars`,
+  which contains...
+* :attr:`~mirgecom.fluid.ConservedVars.momentum`, which is a :class:`numpy.ndarray`
+  of :class:`~numpy.dtype` "object" (an "object array" for short), which contains...
+* :class:`meshmode.dof_array.DOFArray`\ s, i.e. arrays representing a scalar
+  solution field on a :class:`meshmode.discretization.Discretization`. These
+  contain...
+* the actual arrays managed by the array context, typically
+  two-dimensional arrays of shape ``(num_elements, num_dofs_per_element)``.

doc/support/support.rst

@@ -4,6 +4,7 @@ Simulation Support
 
 .. toctree::
 
+   arrays
    time-integrators
    integration-loops
    filter

doc/support/time-integrators.rst

@@ -1,4 +1,8 @@
 Time Integrators
 ================
 
+In addition to the time integrators in this module, :mod:`mirgecom` can also
+make use of the time integrators in :mod:`leap` by way of
+:func:`mirgecom.steppers.advance_state`.
+
 .. automodule:: mirgecom.integrators

examples/README.md

@@ -0,0 +1,16 @@
+# *MIRGE-Com* examples
+
+This directory has a collection of examples to demonstrate and test *MIRGE-Com*
+capabilities.  Examples using the "-mpi" naming convention are MPI parallel and
+are able to run on mulitple GPUs or CPUs in a suitable MPI environemnt. All of
+the example exercise some unique feature of *MIRGE-Com*.  The examples and the
+unique features they exercise are as follows:
+
+- `autoignition-mpi.py`: Chemistry with Pyrometheus
+- `heat-source-mpi.py`: Diffusion operator
+- `lump-mpi.py`: Lump advection, advection verification case
+- `mixture-mpi.py`: Mixture EOS with Pyrometheus
+- `scalar-lump-mpi.py`: Scalar component lump advection verification case
+- `pulse-mpi.py`: Acoustic pulse in a box, wall boundary test case
+- `sod-mpi.py`: Sod's shock case: Fluid test case with strong shock
+- `vortex-mpi.py`: Isentropic vortex advection: outflow boundaries, verification

examples/autoignition-mpi.py

@@ -52,14 +52,16 @@
 from mirgecom.boundary import AdiabaticSlipBoundary
 from mirgecom.initializers import MixtureInitializer
 from mirgecom.eos import PyrometheusMixture
+from mirgecom.inviscid import get_inviscid_cfl
 import cantera
 import pyrometheus as pyro
 
 logger = logging.getLogger(__name__)
 
 
 @mpi_entry_point
-def main(ctx_factory=cl.create_some_context, use_leap=False):
+def main(ctx_factory=cl.create_some_context, casename="autoignition", use_leap=False,
+         restart_step=None, restart_name=None):
     """Drive example."""
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
@@ -73,14 +75,15 @@ def main(ctx_factory=cl.create_some_context, use_leap=False):
     # This example runs only 3 steps by default (to keep CI ~short)
     # With the mixture defined below, equilibrium is achieved at ~40ms
     # To run to equlibrium, set t_final >= 40ms.
-    t_final = 3e-9
+    t_final = 1e-8
     current_cfl = 1.0
     velocity = np.zeros(shape=(dim,))
     current_dt = 1e-9
     current_t = 0
     constant_cfl = False
     nstatus = 1
     nviz = 5
+    nrestart = 5
     rank = 0
     checkpoint_t = current_t
     current_step = 0
@@ -97,12 +100,31 @@ def main(ctx_factory=cl.create_some_context, use_leap=False):
     from mpi4py import MPI
     comm = MPI.COMM_WORLD
     rank = comm.Get_rank()
-
-    from meshmode.mesh.generation import generate_regular_rect_mesh
-    generate_mesh = partial(generate_regular_rect_mesh, a=(box_ll,) * dim,
-                            b=(box_ur,) * dim, nelements_per_axis=(nel_1d,) * dim)
-    local_mesh, global_nelements = generate_and_distribute_mesh(comm, generate_mesh)
-    local_nelements = local_mesh.nelements
+    nproc = comm.Get_size()
+
+    restart_file_pattern = "{casename}-{step:04d}-{rank:04d}.pkl"
+    restart_path = "restart_data/"
+    if restart_step:
+        if not restart_name:
+            restart_name = casename
+        rst_filename = (
+            restart_path
+            + restart_file_pattern.format(casename=restart_name,
+                                          step=restart_step, rank=rank)
+        )
+        from mirgecom.restart import read_restart_data
+        restart_data = read_restart_data(actx, rst_filename)
+        local_mesh = restart_data["local_mesh"]
+        local_nelements = local_mesh.nelements
+        global_nelements = restart_data["global_nelements"]
+        assert restart_data["nparts"] == nproc
+    else:
+        from meshmode.mesh.generation import generate_regular_rect_mesh
+        generate_mesh = partial(generate_regular_rect_mesh, a=(box_ll,)*dim,
+                                b=(box_ur,) * dim, nelements_per_axis=(nel_1d,)*dim)
+        local_mesh, global_nelements = generate_and_distribute_mesh(comm,
+                                                                    generate_mesh)
+        local_nelements = local_mesh.nelements
 
     discr = EagerDGDiscretization(
         actx, local_mesh, order=order, mpi_communicator=comm
@@ -160,7 +182,6 @@ def main(ctx_factory=cl.create_some_context, use_leap=False):
     # Create a Pyrometheus EOS with the Cantera soln. Pyrometheus uses Cantera and
     # generates a set of methods to calculate chemothermomechanical properties and
     # states for this particular mechanism.
-    casename = "autoignition"
     pyrometheus_mechanism = pyro.get_thermochem_class(cantera_soln)(actx.np)
     eos = PyrometheusMixture(pyrometheus_mechanism,
                              temperature_guess=init_temperature)
@@ -178,7 +199,14 @@ def main(ctx_factory=cl.create_some_context, use_leap=False):
 
     my_boundary = AdiabaticSlipBoundary()
     boundaries = {BTAG_ALL: my_boundary}
-    current_state = initializer(eos=eos, x_vec=nodes, time=0)
+
+    if restart_step:
+        current_t = restart_data["t"]
+        current_step = restart_step
+        current_state = restart_data["state"]
+    else:
+        # Set the current state from time 0
+        current_state = initializer(eos=eos, x_vec=nodes, time=0)
 
     # Inspection at physics debugging time
     if debug:
@@ -223,10 +251,34 @@ def my_rhs(t, state):
                 + eos.get_species_source_terms(state))
 
     def my_checkpoint(step, t, dt, state):
+        if check_step(step, nrestart) and step != restart_step:
+            rst_filename = (
+                restart_path
+                + restart_file_pattern.format(casename=casename, step=step,
+                                              rank=rank)
+            )
+            rst_data = {
+                "local_mesh": local_mesh,
+                "state": current_state,
+                "t": t,
+                "step": step,
+                "global_nelements": global_nelements,
+                "num_parts": nproc
+            }
+            from mirgecom.restart import write_restart_file
+            write_restart_file(actx, rst_data, rst_filename, comm)
+
+        # awful - computes potentially expensive viz quantities
+        #         regardless of whether it is time to viz
         reaction_rates = eos.get_production_rates(state)
+        local_cfl = get_inviscid_cfl(discr, eos=eos, dt=current_dt, cv=state)
+        viz_fields = [("reaction_rates", reaction_rates),
+                      ("cfl", local_cfl)]
+        from grudge.op import nodal_max
+        max_cfl = nodal_max(discr, "vol", local_cfl)
         viz_fields = [("reaction_rates", reaction_rates)]
         return sim_checkpoint(discr, visualizer, eos, cv=state,
-                              vizname=casename, step=step,
+                              vizname=casename, step=step, cfl=max_cfl,
                               t=t, dt=dt, nstatus=nstatus, nviz=nviz,
                               constant_cfl=constant_cfl, comm=comm,
                               viz_fields=viz_fields)