Port THGs limiter to main

.ci-support/production-testing-env.sh

@@ -8,7 +8,7 @@ set -x
 # The production capability may be in a CEESD-local mirgecom branch or in a
 # fork, and is specified through the following settings:
 #
-# export PRODUCTION_BRANCH=""   # The base production branch to be installed by emirge
+export PRODUCTION_BRANCH="production-thg"   # The base production branch to be installed by emirge
 # export PRODUCTION_FORK=""  # The fork/home of production changes (if any)
 #
 # Multiple production drivers are supported. The user should provide a ':'-delimited

.github/workflows/ci.yaml

@@ -79,8 +79,8 @@ jobs:
             MINIFORGE_INSTALL_DIR=.miniforge3
             . "$MINIFORGE_INSTALL_DIR/bin/activate" testing
             cd test
-            python -m pip install pytest-xdist
-            python -m pytest -n auto --cov=mirgecom --durations=0 --tb=native --junitxml=pytest.xml --doctest-modules -rxsw . ../doc/*.rst ../doc/*/*.rst
+            python -m pip install pytest
+            python -m pytest --cov=mirgecom --durations=0 --tb=native --junitxml=pytest.xml --doctest-modules -rxsw . ../doc/*.rst ../doc/*/*.rst
 
     examples:
         name: Examples

doc/misc.rst

@@ -71,3 +71,5 @@ 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>`__
+.. [Zhang_2011] Zhang and Shu, Proceedings of The Royal Society, 467 \
+    `(DOI) <https://doi.org/10.1098/rspa.2011.0153>`__

doc/operators/limiters.rst

@@ -0,0 +1,4 @@
+Limiters
+========
+
+.. automodule:: mirgecom.limiter

doc/operators/operators.rst

@@ -7,3 +7,4 @@ Operators
    wave-eq
    diffusion
    gas-dynamics
+   limiters

examples/autoignition-mpi.py

@@ -57,7 +57,7 @@
 from mirgecom.initializers import MixtureInitializer
 from mirgecom.eos import PyrometheusMixture
 from mirgecom.gas_model import GasModel
-from arraycontext import thaw
+from arraycontext import thaw, freeze
 
 from mirgecom.logging_quantities import (
     initialize_logmgr,
@@ -558,10 +558,20 @@ def my_post_step(step, t, dt, state):
             logmgr.tick_after()
         return make_obj_array([cv, fluid_state.temperature]), dt
 
+    from mirgecom.limiter import limiter_liu_osher
+
+    def my_limiter(cv):
+        actx = cv.array_context
+        return thaw(freeze(
+            cv.replace(
+                species_mass=cv.mass*limiter_liu_osher(discr,
+                                                       cv.species_mass_fractions)
+            ), actx), actx)
+
     def my_rhs(t, state):
         cv, tseed = state
         from mirgecom.gas_model import make_fluid_state
-        fluid_state = make_fluid_state(cv=cv, gas_model=gas_model,
+        fluid_state = make_fluid_state(cv=my_limiter(cv), gas_model=gas_model,
                                        temperature_seed=tseed)
         return make_obj_array([
             euler_operator(discr, state=fluid_state, time=t, boundaries=boundaries,

examples/pulse-mpi.py

@@ -36,7 +36,7 @@
     PytatoPyOpenCLArrayContext
 )
 from mirgecom.profiling import PyOpenCLProfilingArrayContext
-from arraycontext import thaw
+from arraycontext import thaw, freeze
 from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
 from grudge.eager import EagerDGDiscretization
 from grudge.shortcuts import make_visualizer
@@ -60,6 +60,7 @@
     GasModel,
     make_fluid_state
 )
+from mirgecom.limiter import limiter_liu_osher
 from logpyle import IntervalTimer, set_dt
 from mirgecom.euler import extract_vars_for_logging, units_for_logging
 from mirgecom.logging_quantities import (
@@ -130,7 +131,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True,
     nviz = 10
     nhealth = 1
 
-    dim = 3
+    dim = 2
     rst_path = "restart_data/"
     rst_pattern = (
         rst_path + "{cname}-{step:04d}-{rank:04d}.pkl"
@@ -312,8 +313,12 @@ def my_post_step(step, t, dt, state):
             logmgr.tick_after()
         return state, dt
 
+    def my_limiter(cv):
+        return thaw(freeze(cv.replace(mass=limiter_liu_osher(discr, cv.mass)), actx),
+                    actx)
+
     def my_rhs(t, state):
-        fluid_state = make_fluid_state(cv=state, gas_model=gas_model)
+        fluid_state = make_fluid_state(cv=my_limiter(state), gas_model=gas_model)
         return euler_operator(discr, state=fluid_state, time=t,
                               boundaries=boundaries,
                               gas_model=gas_model,
@@ -325,7 +330,8 @@ def my_rhs(t, state):
     current_step, current_t, current_cv = \
         advance_state(rhs=my_rhs, timestepper=timestepper,
                       pre_step_callback=my_pre_step,
-                      post_step_callback=my_post_step, dt=current_dt,
+                      post_step_callback=my_post_step,
+                      dt=current_dt,
                       state=current_cv, t=current_t, t_final=t_final)
 
     # Dump the final data

mirgecom/integrators/__init__.py

@@ -26,6 +26,7 @@
 
 from .explicit_rk import rk4_step                          # noqa: F401
 from .lsrk import euler_step, lsrk54_step, lsrk144_step    # noqa: F401
+from .ssprk import ssprk43_step                            # noqa: F401
 
 __doc__ = """
 .. automodule:: mirgecom.integrators.explicit_rk

mirgecom/integrators/ssprk.py

@@ -0,0 +1,40 @@
+"""Timestepping routines for strong-stability preserving Runge-Kutta methods.
+
+.. autofunction:: ssprk43_step
+"""
+
+__copyright__ = """
+Copyright (C) 2022 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.
+"""
+
+
+def ssprk43_step(state, t, dt, rhs):
+    """Take one step using an explicit 4-stage, 3rd-order, SSPRK method."""
+
+    def rhs_update(t, y):
+        return y + dt*rhs(t, y)
+
+    y1 = 1/2*state + 1/2*rhs_update(t, state)
+    y2 = 1/2*y1 + 1/2*rhs_update(t + dt/2, y1)
+    y3 = 2/3*state + 1/6*y2 + 1/6*rhs_update(t + dt, y2)
+    return 1/2*y3 + 1/2*rhs_update(t + dt/2, y3)

mirgecom/limiter.py

@@ -0,0 +1,95 @@
+""":mod:`mirgecom.limiter` is for limiters and limiter-related constructs.
+
+Field limiter functions
+-----------------------
+
+.. autofunction:: limiter_liu_osher
+
+"""
+
+__copyright__ = """
+Copyright (C) 2022 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.
+"""
+
+from functools import partial
+
+from arraycontext import map_array_container
+from meshmode.dof_array import DOFArray
+from grudge.discretization import DiscretizationCollection
+import grudge.op as op
+
+
+def limiter_liu_osher(dcoll: DiscretizationCollection, field):
+    """Implement the positivity-preserving limiter of Liu and Osher (1996).
+
+    The limiter is summarized in the review paper [Zhang_2011]_, Section 2.3,
+    equation 2.9, which uses a linear scaling factor.
+
+    .. note:
+        This limiter is applied only to mass fields
+        (e.g. mass or species masses for multi-component flows)
+
+    Parameters
+    ----------
+    dcoll: :class:`grudge.discretization.DiscretizationCollection`
+        Grudge discretization with boundaries object
+    field: meshmode.dof_array.DOFArray or numpy.ndarray
+        A field or collection of scalar fields to limit
+
+    Returns
+    -------
+    meshmode.dof_array.DOFArray or numpy.ndarray
+        An array container containing the limited field(s).
+    """
+    from grudge.geometry import area_element
+
+    pick_apart = partial(limiter_liu_osher, dcoll)
+
+    if not isinstance(field, DOFArray):
+        # vecs is not a DOFArray -> treat as array container
+        return map_array_container(pick_apart, field)
+
+    actx = field.array_context
+    # Compute nodal and elementwise max/mins of the field
+    _mmax = op.nodal_max(dcoll, "banana", field)
+    _mmin = op.nodal_min(dcoll, "avocado", field)
+    _mmax_i = op.elementwise_max(dcoll, field)
+    _mmin_i = op.elementwise_min(dcoll, field)
+
+    # Compute cell averages of the state
+    inv_area_elements = 1./area_element(
+        actx, dcoll,
+        _use_geoderiv_connection=actx.supports_nonscalar_broadcasting)
+    cell_avgs = \
+        inv_area_elements * op.elementwise_integral(dcoll, field)
+
+    # Compute minmod factor (Eq. 2.9)
+    theta = actx.np.minimum(
+        1.,
+        actx.np.minimum(
+            abs((_mmax - cell_avgs)/(_mmax_i - cell_avgs)),
+            abs((_mmin - cell_avgs)/(_mmin_i - cell_avgs))
+        )
+    )
+
+    return theta*(field - cell_avgs) + cell_avgs