Add entropy-stable dg operators

doc/misc.rst

@@ -87,4 +87,9 @@ References
 .. [Giles_1988] Michael Giles (1988), Non-Reflecting Boundary Conditions for the Euler \
     Equations, CFDL-TR-88-1
 .. [Lachaud_2014] Jean Lachaud and Nagi Mansour (2014), Porous-Material Analysis Toolbox Based
-    on OpenFOAM and Applications, Journal of Thermophysics and Heat Transfer 28 2 
+    on OpenFOAM and Applications, Journal of Thermophysics and Heat Transfer 28 2
+.. [Chandrashekar_2013] Praveen Chandrashekar, Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes \
+   for Compressible Euler and Navier-Stokes Equations, Communications in Computational Physics 14, 5 \
+   `(DOI) <https://doi.org/10.4208/cicp.170712.010313a>`__
+.. [Renac_2021] Florent Renac, Entropy stable, robust and high-order DGSEM for the compressible multicomponent \
+   Euler equations, Journal of Computational Physics, 445 `(DOI) <https://doi.org/10.1016/j.jcp.2021.110584>`__

examples/autoignition-mpi.py

@@ -38,7 +38,8 @@
 from mirgecom.simutil import (
     get_sim_timestep,
     generate_and_distribute_mesh,
-    write_visfile
+    write_visfile,
+    ApplicationOptionsError
 )
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
@@ -74,10 +75,9 @@ class MyRuntimeError(RuntimeError):
 
 
 @mpi_entry_point
-def main(actx_class, use_logmgr=True,
-         use_leap=False, use_overintegration=False,
+def main(actx_class, use_leap=False, use_overintegration=False,
          casename=None, rst_filename=None, log_dependent=True,
-         viscous_terms_on=False):
+         viscous_terms_on=False, use_esdg=False):
     """Drive example."""
     if casename is None:
         casename = "mirgecom"
@@ -90,7 +90,7 @@ def main(actx_class, use_logmgr=True,
     from mirgecom.simutil import global_reduce as _global_reduce
     global_reduce = partial(_global_reduce, comm=comm)
 
-    logmgr = initialize_logmgr(use_logmgr,
+    logmgr = initialize_logmgr(True,
         filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm)
 
     from mirgecom.array_context import initialize_actx, actx_class_is_profiling
@@ -584,10 +584,16 @@ def my_post_step(step, t, dt, state):
 
         return make_obj_array([cv, fluid_state.temperature]), dt
 
-    from mirgecom.inviscid import inviscid_facial_flux_rusanov as inv_num_flux_func
+    from mirgecom.inviscid import (
+        inviscid_facial_flux_rusanov,
+        entropy_stable_inviscid_facial_flux_rusanov
+    )
     from mirgecom.gas_model import make_operator_fluid_states
     from mirgecom.navierstokes import ns_operator
 
+    inv_num_flux_func = entropy_stable_inviscid_facial_flux_rusanov if use_esdg \
+        else inviscid_facial_flux_rusanov
+
     fluid_operator = euler_operator
     if viscous_terms_on:
         fluid_operator = ns_operator
@@ -606,8 +612,9 @@ def my_rhs(t, state):
         fluid_rhs = fluid_operator(
             dcoll, state=fluid_state, gas_model=gas_model, time=t,
             boundaries=boundaries, operator_states_quad=fluid_operator_states,
-            quadrature_tag=quadrature_tag,
+            quadrature_tag=quadrature_tag, use_esdg=use_esdg,
             inviscid_numerical_flux_func=inv_num_flux_func)
+
         chem_rhs = eos.get_species_source_terms(cv, fluid_state.temperature)
         tseed_rhs = fluid_state.temperature - tseed
         cv_rhs = fluid_rhs + chem_rhs
@@ -662,8 +669,8 @@ def my_rhs(t, state):
                         help="turns on compressible Navier-Stokes RHS")
     parser.add_argument("--profiling", action="store_true",
         help="turn on detailed performance profiling")
-    parser.add_argument("--log", action="store_true", default=True,
-        help="turn on logging")
+    parser.add_argument("--esdg", action="store_true",
+        help="use flux-differencing/entropy stable DG for inviscid computations.")
     parser.add_argument("--leap", action="store_true",
         help="use leap timestepper")
     parser.add_argument("--numpy", action="store_true",
@@ -673,6 +680,13 @@ def my_rhs(t, state):
     args = parser.parse_args()
     from warnings import warn
     warn("Automatically turning off DV logging. MIRGE-Com Issue(578)")
+
+    if args.esdg:
+        if not args.lazy and not args.numpy:
+            raise ApplicationOptionsError("ESDG requires lazy or numpy context.")
+        if not args.overintegration:
+            warn("ESDG requires overintegration, enabling --overintegration.")
+
     log_dependent = False
     viscous_terms_on = args.navierstokes
 
@@ -687,9 +701,9 @@ def my_rhs(t, state):
     if args.restart_file:
         rst_filename = args.restart_file
 
-    main(actx_class, use_logmgr=args.log, use_leap=args.leap,
-         use_overintegration=args.overintegration,
-         casename=casename, rst_filename=rst_filename,
+    main(actx_class, use_leap=args.leap,
+         use_overintegration=args.overintegration or args.esdg,
+         casename=casename, rst_filename=rst_filename, use_esdg=args.esdg,
          log_dependent=log_dependent, viscous_terms_on=args.navierstokes)
 
 # vim: foldmethod=marker

examples/combozzle-mpi.py

@@ -29,8 +29,6 @@
 import numpy as np
 from functools import partial
 
-from meshmode.array_context import PyOpenCLArrayContext
-
 from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
 from grudge.shortcuts import make_visualizer
 from grudge.dof_desc import BoundaryDomainTag, DISCR_TAG_QUAD
@@ -41,12 +39,14 @@
 from mirgecom.euler import extract_vars_for_logging, units_for_logging
 from mirgecom.euler import euler_operator
 from mirgecom.navierstokes import ns_operator
+
 from mirgecom.simutil import (
     get_sim_timestep,
     generate_and_distribute_mesh,
     write_visfile,
     force_evaluation,
-    get_box_mesh
+    get_box_mesh,
+    ApplicationOptionsError
 )
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
@@ -146,11 +146,10 @@ def __call__(self, x_vec, *, time=0.0):
 
 
 @mpi_entry_point
-def main(use_logmgr=True,
+def main(actx_class, rst_filename=None,
          use_overintegration=False, casename=None,
-         rst_filename=None, actx_class=PyOpenCLArrayContext,
          log_dependent=False, input_file=None,
-         force_eval=True):
+         force_eval=True, use_esdg=False):
     """Drive example."""
     if casename is None:
         casename = "mirgecom"
@@ -171,7 +170,7 @@ def main(use_logmgr=True,
 
     # {{{ Some discretization parameters
 
-    dim = 3
+    dim = 2
     order = 3
 
     # - scales the size of the domain
@@ -675,7 +674,7 @@ def vol_max(x):
 
     casename = f"{casename}-d{dim}p{order}e{global_nelements}n{nparts}"
 
-    logmgr = initialize_logmgr(use_logmgr,
+    logmgr = initialize_logmgr(True,
         filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm)
 
     if logmgr:
@@ -1077,7 +1076,12 @@ def my_post_step(step, t, dt, state):
 
         return state, dt
 
-    from mirgecom.inviscid import inviscid_facial_flux_rusanov
+    from mirgecom.inviscid import (
+        inviscid_facial_flux_rusanov,
+        entropy_stable_inviscid_facial_flux_rusanov
+    )
+    inv_num_flux_func = entropy_stable_inviscid_facial_flux_rusanov if use_esdg \
+        else inviscid_facial_flux_rusanov
 
     def dummy_pre_step(step, t, dt, state):
         if logmgr:
@@ -1106,18 +1110,18 @@ def cfd_rhs(t, state):
         from mirgecom.gas_model import make_fluid_state
         fluid_state = make_fluid_state(cv=cv, gas_model=gas_model,
                                        temperature_seed=tseed)
-        fluid_operator_states = make_operator_fluid_states(dcoll, fluid_state,
-                                                           gas_model, boundaries,
-                                                           quadrature_tag)
+        fluid_operator_states = make_operator_fluid_states(
+            dcoll, fluid_state, gas_model, boundaries, quadrature_tag=quadrature_tag)
 
         if inviscid_only:
             fluid_rhs = \
                 euler_operator(
                     dcoll, state=fluid_state, time=t,
                     boundaries=boundaries, gas_model=gas_model,
-                    inviscid_numerical_flux_func=inviscid_facial_flux_rusanov,
+                    inviscid_numerical_flux_func=inv_num_flux_func,
                     quadrature_tag=quadrature_tag,
-                    operator_states_quad=fluid_operator_states)
+                    operator_states_quad=fluid_operator_states,
+                    use_esdg=use_esdg)
         else:
             grad_cv = grad_cv_operator(dcoll, gas_model, boundaries, fluid_state,
                                        time=t,
@@ -1128,33 +1132,29 @@ def cfd_rhs(t, state):
                 ns_operator(
                     dcoll, state=fluid_state, time=t, boundaries=boundaries,
                     gas_model=gas_model, quadrature_tag=quadrature_tag,
-                    inviscid_numerical_flux_func=inviscid_facial_flux_rusanov)
+                    inviscid_numerical_flux_func=inv_num_flux_func,
+                    operator_states_quad=fluid_operator_states,
+                    use_esdg=use_esdg)
 
-        if inert_only:
-            chem_rhs = 0*fluid_rhs
-        else:
-            chem_rhs = eos.get_species_source_terms(cv, fluid_state.temperature)
+        if not inert_only:
+            fluid_rhs = fluid_rhs + eos.get_species_source_terms(
+                cv, fluid_state.temperature)
 
         if av_on:
             alpha_f = compute_av_alpha_field(fluid_state)
             indicator = smoothness_indicator(dcoll, fluid_state.mass_density,
                                              kappa=kappa_sc, s0=s0_sc)
-            av_rhs = av_laplacian_operator(
+            fluid_rhs = fluid_rhs + av_laplacian_operator(
                 dcoll, fluid_state=fluid_state, boundaries=boundaries, time=t,
                 gas_model=gas_model, grad_cv=grad_cv,
                 operator_states_quad=fluid_operator_states,
                 alpha=alpha_f, s0=s0_sc, kappa=kappa_sc,
                 indicator=indicator)
-        else:
-            av_rhs = 0*fluid_rhs
 
         if sponge_on:
-            sponge_rhs = _sponge(fluid_state.cv)
-        else:
-            sponge_rhs = 0*fluid_rhs
+            fluid_rhs = fluid_rhs + _sponge(fluid_state.cv)
 
-        fluid_rhs = fluid_rhs + chem_rhs + av_rhs + sponge_rhs
-        tseed_rhs = fluid_state.temperature - tseed
+        tseed_rhs = actx.zeros_like(fluid_state.temperature)
 
         return make_obj_array([fluid_rhs, tseed_rhs])
 
@@ -1245,17 +1245,25 @@ def dummy_rhs(t, state):
         help="Turn off force lazy eval between timesteps")
     parser.add_argument("--profiling", action="store_true",
         help="turn on detailed performance profiling")
-    parser.add_argument("--log", action="store_true", default=True,
-        help="turn on logging")
+    parser.add_argument("--esdg", action="store_true",
+        help="use entropy-stable for inviscid terms")
     parser.add_argument("--leap", action="store_true",
         help="use leap timestepper")
     parser.add_argument("--numpy", action="store_true",
         help="use numpy-based eager actx.")
     parser.add_argument("--restart_file", help="root name of restart file")
     parser.add_argument("--casename", help="casename to use for i/o")
     args = parser.parse_args()
+
     from warnings import warn
     warn("Automatically turning off DV logging. MIRGE-Com Issue(578)")
+
+    if args.esdg:
+        if not args.lazy and not args.numpy:
+            raise ApplicationOptionsError("ESDG requires lazy or numpy context.")
+        if not args.overintegration:
+            warn("ESDG requires overintegration, enabling --overintegration.")
+
     log_dependent = False
     force_eval = not args.no_force
 
@@ -1279,9 +1287,9 @@ def dummy_rhs(t, state):
 
     print(f"Calling main: {time.ctime(time.time())}")
 
-    main(use_logmgr=args.log, input_file=input_file,
-         use_overintegration=args.overintegration,
-         casename=casename, rst_filename=rst_filename, actx_class=actx_class,
-         log_dependent=log_dependent, force_eval=force_eval)
+    main(actx_class, input_file=input_file,
+         use_overintegration=args.overintegration or args.esdg,
+         casename=casename, rst_filename=rst_filename,
+         log_dependent=log_dependent, force_eval=force_eval, use_esdg=args.esdg)
 
 # vim: foldmethod=marker

examples/doublemach-mpi.py

@@ -47,7 +47,7 @@
 from mirgecom.initializers import DoubleMachReflection
 from mirgecom.eos import IdealSingleGas
 from mirgecom.transport import SimpleTransport
-from mirgecom.simutil import get_sim_timestep
+from mirgecom.simutil import get_sim_timestep, ApplicationOptionsError
 from logpyle import set_dt
 from mirgecom.euler import extract_vars_for_logging, units_for_logging
 from mirgecom.logging_quantities import (
@@ -115,13 +115,10 @@ def get_doublemach_mesh():
 
 
 @mpi_entry_point
-def main(use_logmgr=True,
+def main(actx_class, use_esdg=False,
          use_leap=False, use_overintegration=False,
-         casename=None, rst_filename=None, actx_class=None):
+         casename=None, rst_filename=None):
     """Drive the example."""
-    if actx_class is None:
-        raise RuntimeError("Array context class missing.")
-
     if casename is None:
         casename = "mirgecom"
 
@@ -130,7 +127,7 @@ def main(use_logmgr=True,
     rank = comm.Get_rank()
     nparts = comm.Get_size()
 
-    logmgr = initialize_logmgr(use_logmgr,
+    logmgr = initialize_logmgr(True,
         filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm)
 
     from mirgecom.array_context import initialize_actx, actx_class_is_profiling
@@ -391,7 +388,8 @@ def my_rhs(t, state):
         return (
             euler_operator(dcoll, state=fluid_state, time=t,
                            boundaries=boundaries,
-                           gas_model=gas_model, quadrature_tag=quadrature_tag)
+                           gas_model=gas_model, quadrature_tag=quadrature_tag,
+                           use_esdg=use_esdg)
             + av_laplacian_operator(dcoll, fluid_state=fluid_state,
                                     boundaries=boundaries,
                                     time=t, gas_model=gas_model,
@@ -437,8 +435,8 @@ def my_rhs(t, state):
         help="switch to a lazy computation mode")
     parser.add_argument("--profiling", action="store_true",
         help="turn on detailed performance profiling")
-    parser.add_argument("--log", action="store_true", default=True,
-        help="turn on logging")
+    parser.add_argument("--esdg", action="store_true",
+        help="use flux-differencing/entropy stable DG for inviscid computations.")
     parser.add_argument("--leap", action="store_true",
         help="use leap timestepper")
     parser.add_argument("--numpy", action="store_true",
@@ -447,6 +445,13 @@ def my_rhs(t, state):
     parser.add_argument("--casename", help="casename to use for i/o")
     args = parser.parse_args()
 
+    from warnings import warn
+    if args.esdg:
+        if not args.lazy and not args.numpy:
+            raise ApplicationOptionsError("ESDG requires lazy or numpy context.")
+        if not args.overintegration:
+            warn("ESDG requires overintegration, enabling --overintegration.")
+
     from mirgecom.array_context import get_reasonable_array_context_class
     actx_class = get_reasonable_array_context_class(lazy=args.lazy, distributed=True,
                                                     profiling=args.profiling,
@@ -459,8 +464,8 @@ def my_rhs(t, state):
     if args.restart_file:
         rst_filename = args.restart_file
 
-    main(use_logmgr=args.log, use_leap=args.leap,
-         use_overintegration=args.overintegration,
-         casename=casename, rst_filename=rst_filename, actx_class=actx_class)
+    main(actx_class, use_leap=args.leap, use_esdg=args.esdg,
+         use_overintegration=args.overintegration or args.esdg,
+         casename=casename, rst_filename=rst_filename)
 
 # vim: foldmethod=marker