Add HLL numerical flux

.github/workflows/ci.yaml

@@ -77,8 +77,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,11 @@ 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.1016/0021-9991(84)90142-6>`__
+.. [Toro_2009] Eleuterio F. Toro (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer \
+    `(DOI) <http://doi.org/10.1007/978-3-540-49834-6>`__
+.. [Mengaldo_2014] G. Mengaldo (2009), A Guide to the Implementation of Boundary Conditions in Compact \
+    High-Order Methods for Compressible Aerodynamics `(DOI) <https://doi.org/10.2514/6.2014-2923>`__

examples/autoignition-mpi.py

@@ -553,6 +553,8 @@ def my_post_step(step, t, dt, state):
             logmgr.tick_after()
         return make_obj_array([cv, fluid_state.temperature]), dt
 
+    from mirgecom.inviscid import inviscid_facial_flux_rusanov
+
     def my_rhs(t, state):
         cv, tseed = state
         from mirgecom.gas_model import make_fluid_state
@@ -561,6 +563,7 @@ def my_rhs(t, state):
         return make_obj_array([
             euler_operator(discr, state=fluid_state, time=t, boundaries=boundaries,
                            gas_model=gas_model,
+                           inviscid_numerical_flux_func=inviscid_facial_flux_rusanov,
                            quadrature_tag=quadrature_tag)
             + eos.get_species_source_terms(cv, fluid_state.temperature),
             0*tseed])

mirgecom/boundary.py

@@ -3,7 +3,7 @@
 Boundary Treatment Interface
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-.. autoclass FluidBoundary
+.. autoclass:: FluidBoundary
 
 Inviscid Boundary Conditions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -38,11 +38,11 @@
 """
 
 import numpy as np
-from meshmode.dof_array import thaw
+from arraycontext import thaw
 from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
 from mirgecom.fluid import make_conserved
 from grudge.trace_pair import TracePair
-from mirgecom.inviscid import inviscid_facial_flux
+from mirgecom.inviscid import inviscid_facial_flux_rusanov  # default num flux
 
 from abc import ABCMeta, abstractmethod
 
@@ -57,6 +57,7 @@ class FluidBoundary(metaclass=ABCMeta):
     def inviscid_divergence_flux(self, discr, btag, eos, cv_minus, dv_minus,
                                  **kwargs):
         """Get the inviscid boundary flux for the divergence operator."""
+        pass
 
 
 class PrescribedFluidBoundary(FluidBoundary):
@@ -71,24 +72,19 @@ def __init__(self,
                  inviscid_boundary_flux_func=None,
                  # returns CV+, to be used in num flux func (prescribed soln)
                  boundary_state_func=None,
-                 # Inviscid facial flux func given CV(+/-)
-                 inviscid_facial_flux_func=None,
                  # Flux to be used in grad(Temperature) op
                  temperature_gradient_flux_func=None,
                  # Function returns boundary temperature_plus
                  boundary_temperature_func=None):
         """Initialize the PrescribedFluidBoundary and methods."""
         self._bnd_state_func = boundary_state_func
         self._inviscid_bnd_flux_func = inviscid_boundary_flux_func
-        self._inviscid_div_flux_func = inviscid_facial_flux_func
         self._bnd_temperature_func = boundary_temperature_func
 
         if not self._inviscid_bnd_flux_func and not self._bnd_state_func:
             from warnings import warn
             warn("Using dummy boundary: copies interior solution.", stacklevel=2)
 
-        if not self._inviscid_div_flux_func:
-            self._inviscid_div_flux_func = inviscid_facial_flux
         if not self._bnd_state_func:
             self._bnd_state_func = self._dummy_state_func
         if not self._bnd_temperature_func:
@@ -100,29 +96,29 @@ def _dummy_temperature_func(self, temperature_minus, **kwargs):
     def _dummy_state_func(self, state_minus, **kwargs):
         return state_minus
 
-    def _boundary_quantity(self, discr, btag, quantity, **kwargs):
+    def _boundary_quantity(self, discr, btag, quantity, local=False, **kwargs):
         """Get a boundary quantity on local boundary, or projected to "all_faces"."""
-        if "local" in kwargs:
-            if kwargs["local"]:
-                return quantity
-        return discr.project(btag, btag.with_dtag("all_faces"), quantity)
+        return quantity if local else discr.project(
+            btag, btag.with_dtag("all_faces"), quantity)
 
     def inviscid_divergence_flux(self, discr, btag, gas_model, state_minus,
+                                 numerical_flux_func=inviscid_facial_flux_rusanov,
                                  **kwargs):
         """Get the inviscid boundary flux for the divergence operator."""
         # This one is when the user specified a function that directly
         # prescribes the flux components at the boundary
         if self._inviscid_bnd_flux_func:
-            return self._inviscid_bnd_flux_func(discr, btag, gas_model, state_minus,
-                                                **kwargs)
+            return self._inviscid_bnd_flux_func(
+                discr, btag, gas_model, state_minus,
+                numerical_flux_func=numerical_flux_func, **kwargs)
 
         state_plus = self._bnd_state_func(discr=discr, btag=btag,
                                           gas_model=gas_model,
                                           state_minus=state_minus, **kwargs)
         boundary_state_pair = TracePair(btag, interior=state_minus,
                                         exterior=state_plus)
-
-        return self._inviscid_div_flux_func(discr, state_tpair=boundary_state_pair)
+        normal = thaw(discr.normal(btag), state_minus.array_context)
+        return numerical_flux_func(boundary_state_pair, gas_model, normal)
 
 
 class DummyBoundary(PrescribedFluidBoundary):
@@ -173,7 +169,7 @@ def adiabatic_slip_state(self, discr, btag, gas_model, state_minus, **kwargs):
         actx = state_minus.array_context
 
         # Grab a unit normal to the boundary
-        nhat = thaw(actx, discr.normal(btag))
+        nhat = thaw(discr.normal(btag), actx)
 
         # Subtract out the 2*wall-normal component
         # of velocity from the velocity at the wall to

mirgecom/euler.py

@@ -63,8 +63,10 @@
 
 from mirgecom.inviscid import (
     inviscid_flux,
-    inviscid_facial_flux
+    inviscid_facial_flux_rusanov,
+    inviscid_flux_on_element_boundary
 )
+
 from mirgecom.operators import div_operator
 from mirgecom.gas_model import (
     project_fluid_state,
@@ -81,6 +83,7 @@ class _EulerTseedTag:
 
 
 def euler_operator(discr, state, gas_model, boundaries, time=0.0,
+                   inviscid_numerical_flux_func=inviscid_facial_flux_rusanov,
                    quadrature_tag=None):
     r"""Compute RHS of the Euler flow equations.
 
@@ -176,20 +179,13 @@ def _interp_to_surf_quad(utpair):
 
     # Compute volume contributions
     inviscid_flux_vol = inviscid_flux(vol_state_quad)
+
     # Compute interface contributions
-    inviscid_flux_bnd = (
-
-        # Interior faces
-        sum(inviscid_facial_flux(discr, state_pair)
-            for state_pair in interior_states_quad)
-
-        # Domain boundary faces
-        + sum(
-            boundaries[btag].inviscid_divergence_flux(
-                discr, as_dofdesc(btag).with_discr_tag(quadrature_tag), gas_model,
-                state_minus=boundary_states_quad[btag], time=time)
-            for btag in boundaries)
-    )
+    inviscid_flux_bnd = \
+        inviscid_flux_on_element_boundary(
+            discr, gas_model, boundaries, interior_states_quad,
+            boundary_states_quad, quadrature_tag=quadrature_tag,
+            numerical_flux_func=inviscid_numerical_flux_func, time=time)
 
     return -div_operator(discr, dd_quad_vol, dd_quad_faces,
                          inviscid_flux_vol, inviscid_flux_bnd)