Enable CFL and DT

examples/autoignition-mpi.py

@@ -44,14 +44,15 @@
     generate_and_distribute_mesh,
     ExactSolutionMismatch
 )
+from mirgecom.inviscid import get_inviscid_cfl
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
-
 from mirgecom.integrators import rk4_step
 from mirgecom.steppers import advance_state
 from mirgecom.boundary import AdiabaticSlipBoundary
 from mirgecom.initializers import MixtureInitializer
 from mirgecom.eos import PyrometheusMixture
+
 import cantera
 import pyrometheus as pyro
 
@@ -270,9 +271,13 @@ def my_checkpoint(step, t, dt, state):
         # awful - computes potentially expensive viz quantities
         #         regardless of whether it is time to viz
         reaction_rates = eos.get_production_rates(state)
-        viz_fields = [("reaction_rates", reaction_rates)]
+        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)
         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)

examples/lump-mpi.py

@@ -51,6 +51,7 @@
 from mirgecom.boundary import PrescribedBoundary
 from mirgecom.initializers import Lump
 from mirgecom.eos import IdealSingleGas
+from mirgecom.inviscid import get_inviscid_cfl
 
 
 logger = logging.getLogger(__name__)
@@ -131,9 +132,16 @@ def my_rhs(t, state):
                               boundaries=boundaries, eos=eos)
 
     def my_checkpoint(step, t, dt, state):
+        from grudge.op import nodal_max
+        local_cfl = get_inviscid_cfl(discr, eos=eos, dt=current_dt, cv=state)
+        max_cfl = nodal_max(discr, "vol", local_cfl)
+        viz_fields = [
+            ("cfl", local_cfl)
+        ]
         return sim_checkpoint(discr, visualizer, eos, cv=state,
-                              exact_soln=initializer, vizname=casename, step=step,
-                              t=t, dt=dt, nstatus=nstatus, nviz=nviz,
+                              viz_fields=viz_fields, exact_soln=initializer,
+                              vizname=casename, step=step,
+                              t=t, dt=dt, nstatus=nstatus, nviz=nviz, cfl=max_cfl,
                               exittol=exittol, constant_cfl=constant_cfl, comm=comm)
 
     try:

examples/mixture-fixed-cfl-mpi.py

@@ -0,0 +1,315 @@
+"""Demonstrate combustive mixture with constant CFL mode."""
+
+__copyright__ = """
+Copyright (C) 2020 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.
+"""
+import logging
+import numpy as np
+import pyopencl as cl
+import pyopencl.tools as cl_tools
+from functools import partial
+
+from meshmode.array_context import PyOpenCLArrayContext
+from meshmode.dof_array import thaw
+from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
+from grudge.eager import EagerDGDiscretization
+from grudge.shortcuts import make_visualizer
+
+
+from mirgecom.euler import euler_operator
+from mirgecom.simutil import (
+    check_step,
+    generate_and_distribute_mesh,
+    ExactSolutionMismatch
+)
+from mirgecom.inviscid import (
+    get_inviscid_timestep,
+    get_inviscid_cfl
+)
+from mirgecom.io import make_init_message
+from mirgecom.mpi import mpi_entry_point
+from mirgecom.integrators import rk4_step
+from mirgecom.steppers import advance_state
+from mirgecom.boundary import AdiabaticSlipBoundary
+from mirgecom.initializers import MixtureInitializer
+from mirgecom.eos import PyrometheusMixture
+
+import cantera
+import pyrometheus as pyro
+
+logger = logging.getLogger(__name__)
+
+
+@mpi_entry_point
+def main(ctx_factory=cl.create_some_context, use_leap=False):
+    """Drive example."""
+    cl_ctx = ctx_factory()
+    queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue,
+            allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)))
+
+    dim = 2
+    nel_1d = 8
+    order = 1
+
+    # 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 = 1e-7
+    current_cfl = 0.01
+    velocity = np.zeros(shape=(dim,))
+    current_dt = 1e-9
+    current_t = 0
+    constant_cfl = True
+    nstatus = 1
+    nviz = 2
+    rank = 0
+    checkpoint_t = current_t
+    current_step = 0
+    if use_leap:
+        from leap.rk import RK4MethodBuilder
+        timestepper = RK4MethodBuilder("state")
+    else:
+        timestepper = rk4_step
+    box_ll = -0.005
+    box_ur = 0.005
+    error_state = False
+    debug = 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
+
+    discr = EagerDGDiscretization(
+        actx, local_mesh, order=order, mpi_communicator=comm
+    )
+    nodes = thaw(actx, discr.nodes())
+
+    # {{{  Set up initial state using Cantera
+
+    # Use Cantera for initialization
+    # -- Pick up a CTI for the thermochemistry config
+    # --- Note: Users may add their own CTI file by dropping it into
+    # ---       mirgecom/mechanisms alongside the other CTI files.
+    from mirgecom.mechanisms import get_mechanism_cti
+    mech_cti = get_mechanism_cti("uiuc")
+
+    cantera_soln = cantera.Solution(phase_id="gas", source=mech_cti)
+    nspecies = cantera_soln.n_species
+
+    # Initial temperature, pressure, and mixutre mole fractions are needed to
+    # set up the initial state in Cantera.
+    init_temperature = 1500.0  # Initial temperature hot enough to burn
+    # Parameters for calculating the amounts of fuel, oxidizer, and inert species
+    equiv_ratio = 1.0
+    ox_di_ratio = 0.21
+    stoich_ratio = 3.0
+    # Grab the array indices for the specific species, ethylene, oxygen, and nitrogen
+    i_fu = cantera_soln.species_index("C2H4")
+    i_ox = cantera_soln.species_index("O2")
+    i_di = cantera_soln.species_index("N2")
+    x = np.zeros(nspecies)
+    # Set the species mole fractions according to our desired fuel/air mixture
+    x[i_fu] = (ox_di_ratio*equiv_ratio)/(stoich_ratio+ox_di_ratio*equiv_ratio)
+    x[i_ox] = stoich_ratio*x[i_fu]/equiv_ratio
+    x[i_di] = (1.0-ox_di_ratio)*x[i_ox]/ox_di_ratio
+    # Uncomment next line to make pylint fail when it can't find cantera.one_atm
+    one_atm = cantera.one_atm  # pylint: disable=no-member
+    # one_atm = 101325.0
+
+    # Let the user know about how Cantera is being initilized
+    print(f"Input state (T,P,X) = ({init_temperature}, {one_atm}, {x}")
+    # Set Cantera internal gas temperature, pressure, and mole fractios
+    cantera_soln.TPX = init_temperature, one_atm, x
+    # Pull temperature, total density, mass fractions, and pressure from Cantera
+    # We need total density, and mass fractions to initialize the fluid/gas state.
+    can_t, can_rho, can_y = cantera_soln.TDY
+    can_p = cantera_soln.P
+    # *can_t*, *can_p* should not differ (significantly) from user's initial data,
+    # but we want to ensure that we use exactly the same starting point as Cantera,
+    # so we use Cantera's version of these data.
+
+    # }}}
+
+    # {{{ Create Pyrometheus thermochemistry object & EOS
+
+    # 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 = "mixture-adaptive"
+    pyrometheus_mechanism = pyro.get_thermochem_class(cantera_soln)(actx.np)
+    eos = PyrometheusMixture(pyrometheus_mechanism,
+                             temperature_guess=init_temperature)
+
+    # }}}
+
+    # {{{ MIRGE-Com state initialization
+
+    # Initialize the fluid/gas state with Cantera-consistent data:
+    # (density, pressure, temperature, mass_fractions)
+    print(f"Cantera state (rho,T,P,Y) = ({can_rho}, {can_t}, {can_p}, {can_y}")
+    initializer = MixtureInitializer(dim=dim, nspecies=nspecies,
+                                     pressure=can_p, temperature=can_t,
+                                     massfractions=can_y, velocity=velocity)
+
+    my_boundary = AdiabaticSlipBoundary()
+    boundaries = {BTAG_ALL: my_boundary}
+    current_state = initializer(eos=eos, x_vec=nodes, t=0)
+
+    # Inspection at physics debugging time
+    if debug:
+        print("Initial MIRGE-Com state:")
+        print(f"{current_state=}")
+        print(f"Initial DV pressure: {eos.pressure(current_state)}")
+        print(f"Initial DV temperature: {eos.temperature(current_state)}")
+
+    # }}}
+
+    visualizer = make_visualizer(discr)
+    initname = initializer.__class__.__name__
+    eosname = eos.__class__.__name__
+    init_message = make_init_message(dim=dim, order=order,
+                                     nelements=local_nelements,
+                                     global_nelements=global_nelements,
+                                     dt=current_dt, t_final=t_final, nstatus=nstatus,
+                                     nviz=nviz, cfl=current_cfl,
+                                     constant_cfl=constant_cfl, initname=initname,
+                                     eosname=eosname, casename=casename)
+
+    # Cantera equilibrate calculates the expected end state @ chemical equilibrium
+    # i.e. the expected state after all reactions
+    cantera_soln.equilibrate("UV")
+    eq_temperature, eq_density, eq_mass_fractions = cantera_soln.TDY
+    eq_pressure = cantera_soln.P
+
+    # Report the expected final state to the user
+    if rank == 0:
+        logger.info(init_message)
+        logger.info(f"Expected equilibrium state:"
+                    f" {eq_pressure=}, {eq_temperature=},"
+                    f" {eq_density=}, {eq_mass_fractions=}")
+
+    def my_rhs(t, state):
+        return (euler_operator(discr, cv=state, t=t,
+                               boundaries=boundaries, eos=eos)
+                + eos.get_species_source_terms(state))
+
+    def mixture_prestep_function(step, t, dt, state):
+        do_viz = check_step(step, nviz)
+        viz_fields = [("cv", state)]
+        current_dt = dt
+
+        if constant_cfl:
+            local_dt = get_inviscid_timestep(discr, eos=eos, cv=state)
+            from grudge.op import nodal_min
+            current_dt = current_cfl * nodal_min(discr, "vol", local_dt)
+        else:  # constant dt mode
+            if do_viz:   # calculate cfl field only if visualizing
+                local_cfl = get_inviscid_cfl(discr, eos=eos, dt=dt, cv=state)
+        if do_viz:  # extend viz field if viz time
+            # Calculate DV only if needed for visualization
+            dv = eos.dependent_vars(state)
+            viz_fields.append(("dv", dv))
+            # only if vizzing, calculate reaction rates
+            reaction_rates = eos.get_production_rates(state)
+            viz_fields.append(("reaction_rates", reaction_rates))
+            if constant_cfl:
+                viz_fields.append(("dt", local_dt))
+            else:  # constant dt mode
+                viz_fields.append(("cfl", local_cfl))
+
+        errors = current_dt < 0 or np.isnan(current_dt) or current_dt == np.inf
+
+        if do_viz or errors:  # write viz at viztime, or if there were errors
+            from mirgecom.io import make_rank_fname, make_par_fname
+            rank_fn = make_rank_fname(basename=casename, rank=rank, step=step, t=t)
+            visualizer.write_parallel_vtk_file(
+                comm, rank_fn, viz_fields, overwrite=True,
+                par_manifest_filename=make_par_fname(
+                    basename=casename, step=step, t=t
+                )
+            )
+
+        if check_step(step, nstatus) or errors:
+            if not do_viz:  # we already have dv on viz steps
+                dv = eos.dependent_vars(state)
+            from grudge.op import nodal_max
+            min_temperature = nodal_min(discr, "vol", dv.temperature)
+            max_temperature = nodal_max(discr, "vol", dv.temperature)
+            min_pressure = nodal_min(discr, "vol", dv.pressure)
+            max_pressure = nodal_max(discr, "vol", dv.pressure)
+            if rank == 0:
+                logger.info(f"\nStep:{step}, Time:{t}, DT:{current_dt},"
+                            f"CFL:{current_cfl}\n"
+                            f"---- P({min_pressure}, {max_pressure})\n"
+                            f"---- T({min_temperature}, {max_temperature})\n")
+
+        # this exit is safe, errors(current_dt) is already collective
+        if errors:
+            logger.info("Fatal error: Invalid simulation DT")
+            import sys
+            sys.exit()
+
+        t_remaining = max(0, t_final - t)
+        return min(t_remaining, current_dt)
+
+    try:
+        (current_step, current_t, current_state) = \
+            advance_state(rhs=my_rhs, timestepper=timestepper,
+                          checkpoint=mixture_prestep_function, state=current_state,
+                          t=current_t, t_final=t_final)
+
+    except ExactSolutionMismatch as ex:
+        error_state = True
+        current_step = ex.step
+        current_t = ex.t
+        current_state = ex.state
+
+    if not check_step(current_step, nviz):  # If final step not an output step
+        if rank == 0:
+            logger.info("Checkpointing final state ...")
+        mixture_prestep_function(current_step, t=current_t,
+                                 dt=(current_t - checkpoint_t),
+                                 state=current_state)
+
+    if current_t - t_final < 0:
+        error_state = True
+
+    if error_state:
+        raise ValueError("Simulation did not complete successfully.")
+
+    if rank == 0:
+        logger.info(f"Simulation finished at time {current_t=}.")
+
+
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO)
+    main(use_leap=False)
+
+# vim: foldmethod=marker