Stab initially at making THG ESDG work.

examples/vortex-mpi.py

@@ -32,7 +32,7 @@
 from grudge.shortcuts import make_visualizer
 
 from mirgecom.discretization import create_discretization_collection
-from mirgecom.euler import euler_operator
+from mirgecom.euler import euler_operator,  entropy_stable_euler_operator
 from mirgecom.simutil import (
     get_sim_timestep,
     generate_and_distribute_mesh,
@@ -68,9 +68,12 @@ class MyRuntimeError(RuntimeError):
 
 
 @mpi_entry_point
-def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True,
-         use_leap=False, use_profiling=False, casename=None, lazy=False,
-         rst_filename=None):
+def main(actx_class, use_logmgr=True, ctx_factory=cl.create_some_context, 
+          use_leap=False, use_profiling=False, casename=None, lazy=False,
+          use_overintegration=False, use_esdg=False,
+          rst_filename=None):
+
+    """Drive the example."""
     """Drive the example."""
     cl_ctx = ctx_factory()
 
@@ -97,6 +100,20 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True,
     from mirgecom.simutil import get_reasonable_memory_pool
     alloc = get_reasonable_memory_pool(cl_ctx, queue)
 
+    if lazy:
+        actx = actx_class(comm, queue, mpi_base_tag=12000, allocator=alloc)
+    else:
+        actx = actx_class(comm, queue, allocator=alloc, force_device_scalars=True)
+
+    if use_esdg and not actx.supports_nonscalar_broadcasting:
+        raise RuntimeError(
+            f"{actx} is not a suitable array context for using flux-differencing. "
+            "The underlying array context must be capable of performing basic "
+            "array broadcasting operations. Use PytatoPyOpenCLArrayContext instead."
+        )
+    from mirgecom.simutil import get_reasonable_memory_pool
+    alloc = get_reasonable_memory_pool(cl_ctx, queue)
+
     if lazy:
         actx = actx_class(comm, queue, mpi_base_tag=12000, allocator=alloc)
     else:
@@ -109,7 +126,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True,
         timestepper = RK4MethodBuilder("state")
     else:
         timestepper = rk4_step
-    t_final = 0.01
+    t_final = 1.0
     current_cfl = 1.0
     current_dt = .001
     current_t = 0
@@ -147,13 +164,19 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True,
         local_mesh, global_nelements = generate_and_distribute_mesh(comm,
                                                                     generate_mesh)
         local_nelements = local_mesh.nelements
+        from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD
+        from meshmode.discretization.poly_element import \
+        default_simplex_group_factory, QuadratureSimplexGroupFactory
 
     order = 3
     dcoll = create_discretization_collection(actx, local_mesh, order=order)
     nodes = actx.thaw(dcoll.nodes())
 
     vis_timer = None
-
+    if use_overintegration:
+        quadrature_tag = DISCR_TAG_QUAD
+    else:
+        quadrature_tag = None
     if logmgr:
         logmgr_add_cl_device_info(logmgr, queue)
         logmgr_add_device_memory_usage(logmgr, queue)
@@ -184,7 +207,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True,
     eos = IdealSingleGas()
     vel = np.zeros(shape=(dim,))
     orig = np.zeros(shape=(dim,))
-    vel[:dim] = 1.0
+    vel[0] = 1.0
     initializer = Vortex2D(center=orig, velocity=vel)
     gas_model = GasModel(eos=eos)
 
@@ -198,7 +221,10 @@ def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs):
     boundaries = {
         BTAG_ALL: PrescribedFluidBoundary(boundary_state_func=boundary_solution)
     }
-
+    if use_esdg:
+        operator_rhs = entropy_stable_euler_operator
+    else:
+        operator_rhs = euler_operator
     if rst_filename:
         current_t = restart_data["t"]
         current_step = restart_data["step"]
@@ -276,11 +302,11 @@ def my_health_check(pressure, component_errors):
         health_error = False
         from mirgecom.simutil import check_naninf_local, check_range_local
         if check_naninf_local(dcoll, "vol", pressure) \
-           or check_range_local(dcoll, "vol", pressure, .2, 1.02):
+           or check_range_local(dcoll, "vol", pressure, .2, 1.2):
             health_error = True
             logger.info(f"{rank=}: Invalid pressure data found.")
 
-        exittol = .1
+        exittol = 1.0
         if max(component_errors) > exittol:
             health_error = True
             if rank == 0:
@@ -356,7 +382,7 @@ def my_post_step(step, t, dt, state):
 
     def my_rhs(t, state):
         fluid_state = make_fluid_state(state, gas_model)
-        return euler_operator(dcoll, state=fluid_state, time=t,
+        return operator_rhs(dcoll, state=fluid_state, time=t,
                               boundaries=boundaries, gas_model=gas_model)
 
     current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt,
@@ -393,6 +419,10 @@ def my_rhs(t, state):
     import argparse
     casename = "vortex"
     parser = argparse.ArgumentParser(description=f"MIRGE-Com Example: {casename}")
+    parser.add_argument("--overintegration", action="store_true",
+        help="use overintegration in the RHS computations"),
+    parser.add_argument("--esdg", action="store_true",
+        help="use flux-differencing/entropy stable DG for inviscid computations.")
     parser.add_argument("--lazy", action="store_true",
         help="switch to a lazy computation mode")
     parser.add_argument("--profiling", action="store_true",
@@ -419,7 +449,13 @@ def my_rhs(t, state):
     if args.restart_file:
         rst_filename = args.restart_file
 
+#    main(use_logmgr=args.log, use_overintegration=args.overintegration,
+#         use_esdg=args.esdg, use_leap=args.leap, use_profiling=args.profiling,
+#         casename=casename, rst_filename=rst_filename, actx_class=actx_class)
+
     main(actx_class, use_logmgr=args.log, use_leap=args.leap, lazy=lazy,
-         use_profiling=args.profiling, casename=casename, rst_filename=rst_filename)
+         use_profiling=args.profiling,
+         use_overintegration=args.overintegration, 
+         use_esdg=args.esdg, casename=casename, rst_filename=rst_filename)
 
 # vim: foldmethod=marker

mirgecom/euler.py

@@ -53,25 +53,41 @@
 """
 
 import numpy as np  # noqa
-
+from arraycontext import map_array_container
 from meshmode.discretization.connection import FACE_RESTR_ALL
 from grudge.dof_desc import (
     DD_VOLUME_ALL,
     VolumeDomainTag,
     DISCR_TAG_BASE,
 )
-
-from mirgecom.gas_model import make_operator_fluid_states
+from grudge.dof_desc import DOFDesc, as_dofdesc
+from mirgecom.gas_model import (
+    project_fluid_state,
+    make_operator_fluid_states,
+    make_fluid_state_trace_pairs,
+    make_entropy_projected_fluid_state,
+    conservative_to_entropy_vars,
+    entropy_to_conservative_vars
+)
 from mirgecom.inviscid import (
     inviscid_flux,
     inviscid_facial_flux_rusanov,
-    inviscid_flux_on_element_boundary
+    inviscid_flux_on_element_boundary,
+    entropy_conserving_flux_chandrashekar,
+    entropy_stable_inviscid_flux_rusanov
 )
+from meshmode.dof_array import DOFArray
 
+from functools import partial
+import grudge.op as op
 from mirgecom.operators import div_operator
 from mirgecom.utils import normalize_boundaries
-
-
+from grudge.projection import volume_quadrature_project
+from grudge.flux_differencing import volume_flux_differencing
+from grudge.trace_pair import (
+    TracePair,
+    interior_trace_pairs
+)
 def euler_operator(dcoll, state, gas_model, boundaries, time=0.0,
                    inviscid_numerical_flux_func=inviscid_facial_flux_rusanov,
                    quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL,
@@ -157,7 +173,162 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0,
     return -div_operator(dcoll, dd_vol_quad, dd_allfaces_quad,
                          inviscid_flux_vol, inviscid_flux_bnd)
 
-
+def entropy_stable_euler_operator(
+        discr, state, gas_model, boundaries, time=0.0,
+        inviscid_numerical_flux_func=entropy_stable_inviscid_flux_rusanov,
+        quadrature_tag=None):
+    """Compute RHS of the Euler flow equations using flux-differencing.
+    Parameters
+    ----------
+    state: :class:`~mirgecom.gas_model.FluidState`
+        Fluid state object with the conserved state, and dependent
+        quantities.
+    boundaries
+        Dictionary of boundary functions, one for each valid btag
+    time
+        Time
+    gas_model: :class:`~mirgecom.gas_model.GasModel`
+        Physical gas model including equation of state, transport,
+        and kinetic properties as required by fluid state
+    quadrature_tag
+        An optional identifier denoting a particular quadrature
+        discretization to use during operator evaluations.
+        The default value is *None*.
+    Returns
+    -------
+    :class:`mirgecom.fluid.ConservedVars`
+        Agglomerated object array of DOF arrays representing the RHS of the Euler
+        flow equations.
+    """
+    dd_base = as_dofdesc("vol")
+    dd_vol = DOFDesc("vol", quadrature_tag)
+    dd_faces = DOFDesc("all_faces", quadrature_tag)
+    # NOTE: For single-gas this is just a fixed scalar.
+    # However, for mixtures, gamma is a DOFArray. For now,
+    # we are re-using gamma from here and *not* recomputing
+    # after applying entropy projections. It is unclear at this
+    # time whether it's strictly necessary or if this is good enough
+    gamma = gas_model.eos.gamma(state.cv, state.temperature)
+    state_quad = project_fluid_state(discr, "vol", dd_vol, state, gas_model)
+
+    # Compute the projected (nodal) entropy variables
+    entropy_vars = volume_quadrature_project(
+        discr, dd_vol,
+        # Map to entropy variables
+        conservative_to_entropy_vars(gamma, state_quad))
+
+    modified_conserved_fluid_state = \
+        make_entropy_projected_fluid_state(discr, dd_vol, dd_faces,
+                                           state, entropy_vars, gamma, gas_model)
+
+    def _reshape(shape, ary):
+        if not isinstance(ary, DOFArray):
+            return map_array_container(partial(_reshape, shape), ary)
+
+        return DOFArray(ary.array_context, data=tuple(
+            subary.reshape(grp.nelements, *shape)
+            # Just need group for determining the number of elements
+            for grp, subary in zip(discr.discr_from_dd("vol").groups, ary)))
+
+    flux_matrices = entropy_conserving_flux_chandrashekar(
+        gas_model,
+        _reshape((1, -1), modified_conserved_fluid_state),
+        _reshape((-1, 1), modified_conserved_fluid_state))
+
+    # Compute volume derivatives using flux differencing
+    inviscid_flux_vol = \
+        -volume_flux_differencing(discr, dd_vol, dd_faces, flux_matrices)
+
+    def interp_to_surf_quad(utpair):
+        local_dd = utpair.dd
+        local_dd_quad = local_dd.with_discr_tag(quadrature_tag)
+        return TracePair(
+            local_dd_quad,
+            interior=op.project(discr, local_dd, local_dd_quad, utpair.int),
+            exterior=op.project(discr, local_dd, local_dd_quad, utpair.ext)
+        )
+
+    tseed_interior_pairs = None
+    if state.is_mixture:
+        # If this is a mixture, we need to exchange the temperature field because
+        # mixture pressure (used in the inviscid flux calculations) depends on
+        # temperature and we need to seed the temperature calculation for the
+        # (+) part of the partition boundary with the remote temperature data.
+        tseed_interior_pairs = [
+            # Get the interior trace pairs onto the surface quadrature
+            # discretization (if any)
+            interp_to_surf_quad(tpair)
+            for tpair in interior_trace_pairs(discr, state.temperature)
+        ]
+
+    def interp_to_surf_modified_conservedvars(gamma, utpair):
+        """Takes a trace pair containing the projected entropy variables
+        and converts them into conserved variables on the quadrature grid.
+        """
+        local_dd = utpair.dd
+        local_dd_quad = local_dd.with_discr_tag(quadrature_tag)
+        # Interpolate entropy variables to the surface quadrature grid
+        vtilde_tpair = op.project(discr, local_dd, local_dd_quad, utpair)
+        if isinstance(gamma, DOFArray):
+            gamma = op.project(discr, dd_base, local_dd_quad, gamma)
+        return TracePair(
+            local_dd_quad,
+            # Convert interior and exterior states to conserved variables
+            interior=entropy_to_conservative_vars(gamma, vtilde_tpair.int),
+            exterior=entropy_to_conservative_vars(gamma, vtilde_tpair.ext)
+        )
+
+    cv_interior_pairs = [
+        # Compute interior trace pairs using modified conservative
+        # variables on the quadrature grid
+        # (obtaining state from projected entropy variables)
+        interp_to_surf_modified_conservedvars(gamma, tpair)
+        for tpair in interior_trace_pairs(discr, entropy_vars)
+    ]
+
+    boundary_states = {
+        # TODO: Use modified conserved vars as the input state?
+        # Would need to make an "entropy-projection" variant
+        # of *project_fluid_state*
+        btag: project_fluid_state(
+            discr, dd_base,
+            # Make sure we get the state on the quadrature grid
+            # restricted to the tag *btag*
+            as_dofdesc(btag).with_discr_tag(quadrature_tag),
+            state, gas_model) for btag in boundaries
+    }
+
+    # Interior interface state pairs consisting of modified conservative
+    # variables and the corresponding temperature seeds
+    interior_states = make_fluid_state_trace_pairs(cv_interior_pairs,
+                                                   gas_model,
+                                                   tseed_interior_pairs)
+
+    # Surface contributions
+    inviscid_flux_bnd = (
+
+        # Domain boundaries
+        sum(boundaries[btag].inviscid_divergence_flux(
+            discr,
+            # Make sure we get the state on the quadrature grid
+            # restricted to the tag *btag*
+            as_dofdesc(btag).with_discr_tag(quadrature_tag),
+            gas_model,
+            state_minus=boundary_states[btag],
+            time=time,
+            numerical_flux_func=inviscid_numerical_flux_func)
+            for btag in boundaries)
+
+        # Interior boundaries (using entropy stable numerical flux)
+        + sum(inviscid_facial_flux(discr, gas_model=gas_model, state_pair=state_pair,
+                                   numerical_flux_func=inviscid_numerical_flux_func)
+              for state_pair in interior_states)
+    )
+
+    return op.inverse_mass(
+        discr,
+        inviscid_flux_vol - op.face_mass(discr, dd_faces, inviscid_flux_bnd)
+    )
 # By default, run unitless
 NAME_TO_UNITS = {
     "mass": "",