Split CV

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

@@ -8,15 +8,15 @@ 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-state-handling"   # 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
 # list of driver locations, where each driver location is of the form:
 # "fork/repo@branch". The defaults are provided below as an example. Provide custom
 # production drivers in this variable:
 #
-# export PRODUCTION_DRIVERS=""
+export PRODUCTION_DRIVERS="illinois-ceesd/drivers_y1-nozzle@parallel-lazy-state-handling:illinois-ceesd/drivers_flame1d@state-handling:illinois-ceesd/drivers_y2-isolator@state-handling"
 #
 # Example:
 # PRODUCTION_DRIVERS="illinois-ceesd/drivers_y1-nozzle@main:w-hagen/isolator@NS"

mirgecom/boundary.py

@@ -219,8 +219,10 @@ def adiabatic_slip_pair(self, discr, cv, btag, **kwargs):
         mom_normcomp = np.dot(int_cv.momentum, nhat)  # wall-normal component
         wnorm_mom = nhat * mom_normcomp  # wall-normal mom vec
         ext_mom = int_cv.momentum - 2.0 * wnorm_mom  # prescribed ext momentum
+        species_mass = int_cv.species_mass if int_cv.has_multispecies else None
 
         # Form the external boundary solution with the new momentum
         ext_cv = make_conserved(dim=dim, mass=int_cv.mass, energy=int_cv.energy,
-                                momentum=ext_mom, species_mass=int_cv.species_mass)
+                                momentum=ext_mom, species_mass=species_mass)
+
         return TracePair(btag, interior=int_cv, exterior=ext_cv)

mirgecom/fluid.py

@@ -4,8 +4,7 @@
 ^^^^^^^^^^^^^^^^^^^^^
 
 .. autoclass:: ConservedVars
-.. autofunction:: split_conserved
-.. autofunction:: join_conserved
+.. autoclass:: MixtureConservedVars
 .. autofunction:: make_conserved
 
 Helper Functions
@@ -61,15 +60,9 @@ class ConservedVars:
     for the canonical conserved quantities (mass, energy, momentum,
     and species masses) per unit volume: $(\rho,\rho{E},\rho\vec{V},
     \rho{Y_s})$ from an agglomerated object array.  This data structure is intimately
-    related to helper functions :func:`join_conserved` and :func:`split_conserved`,
-    which pack and unpack (respectively) quantities to-and-from flat object
-    array representations of the data.
-
-    .. note::
-
-        The use of the terms pack and unpack here is misleading. In reality, there
-        is no data movement when using this dataclass to view your data. This
-        dataclass is entirely about the representation of the data.
+    related to the helper function :func:`make_conserved` which helps form a
+    conservation-law-specific representation of the simulation state data for use by
+    specific operators (fluids in this case).
 
     .. attribute:: dim
 
@@ -94,12 +87,9 @@ class ConservedVars:
         ``(ndim, ndim)`` respectively for scalar or vector quantities corresponding
         to the ndim equations of momentum conservation.
 
-    .. attribute:: species_mass
-
-        Object array (:class:`~numpy.ndarray`) with shape ``(nspecies,)``
-        of :class:`~meshmode.dof_array.DOFArray`, or an object array with shape
-        ``(nspecies, ndim)`` respectively for scalar or vector quantities
-        corresponding to the `nspecies` species mass conservation equations.
+    .. autoattribute:: has_multispecies
+    .. autoattribute:: array_context
+    .. autoattribute:: velocity
 
     :example::
 
@@ -124,11 +114,11 @@ class ConservedVars:
         $N_{\text{eq}}$ :class:`~meshmode.dof_array.DOFArray`.
 
         To use this dataclass for a fluid CV-specific view on the content of
-        $\mathbf{Q}$, one can call :func:`split_conserved` to get a `ConservedVars`
+        $\mathbf{Q}$, one can call :func:`make_conserved` to get a `ConservedVars`
         dataclass object that resolves the fluid CV associated with each conservation
         equation::
 
-            fluid_cv = split_conserved(ndim, Q),
+            fluid_cv = make_conserved(dim=ndim, q=Q),
 
         after which::
 
@@ -143,37 +133,6 @@ class ConservedVars:
         - :mod:`~mirgecom.initializers`
         - :mod:`~mirgecom.simutil`
 
-    :example::
-
-        Use `join_conserved` to create an agglomerated $\mathbf{Q}$ array from the
-        fluid conserved quantities (CV).
-
-        See the first example for the definition of CV, $\mathbf{Q}$, `ndim`,
-        `nspecies`, and $N_{\text{eq}}$.
-
-        Often, a user starts with the fluid conserved quantities like mass and
-        energy densities, and it is desired to glom those quantities together into
-        a *MIRGE*-compatible $\mathbf{Q}$ data structure.
-
-        For example, a solution initialization routine may set the fluid
-        quantities::
-
-            rho = ... # rho is a DOFArray with fluid density
-            v = ... # v is an ndim-vector of DOFArray with components of velocity
-            e = ... # e is a DOFArray with fluid energy
-
-        An agglomerated array of fluid independent variables can then be
-        created with::
-
-            q = join_conserved(ndim, mass=rho, energy=rho*e, momentum=rho*v)
-
-        after which *q* will be an obj array of $N_{\text{eq}}$ DOFArrays containing
-        the fluid conserved state data.
-
-        Examples of this sort of use for `join_conserved` can be found in:
-
-        - :mod:`~mirgecom.initializers`
-
     :example::
 
         Use `ConservedVars` to access a vector quantity for each fluid equation.
@@ -203,11 +162,11 @@ class ConservedVars:
 
         Presuming that `grad_q` is the agglomerated *MIRGE* data structure with the
         gradient data, this dataclass can be used to get a fluid CV-specific view on
-        the content of $\nabla\mathbf{Q}$. One can call :func:`split_conserved` to
+        the content of $\nabla\mathbf{Q}$. One can call :func:`make_conserved` to
         get a `ConservedVars` dataclass object that resolves the vector quantity
         associated with each conservation equation::
 
-            grad_cv = split_conserved(ndim, grad_q),
+            grad_cv = make_conserved(dim=ndim, q=grad_q),
 
         after which::
 
@@ -226,7 +185,11 @@ class ConservedVars:
     mass: DOFArray
     energy: DOFArray
     momentum: np.ndarray
-    species_mass: np.ndarray = np.empty((0,), dtype=object)  # empty = immutable
+
+    @property
+    def has_multispecies(self):
+        """Return whether this state has multiple species."""
+        return False
 
     @property
     def array_context(self):
@@ -243,6 +206,52 @@ def velocity(self):
         """Return the fluid velocity = momentum / mass."""
         return self.momentum / self.mass
 
+    def join(self):
+        """Return a flat object array of the CV components."""
+        return _join_conserved(
+            dim=self.dim,
+            mass=self.mass,
+            energy=self.energy,
+            momentum=self.momentum)
+
+    def __reduce__(self):
+        """Return a tuple reproduction of self for pickling."""
+        return (type(self), tuple(getattr(self, f.name)
+                                  for f in fields(type(self))))
+
+    def replace(self, **kwargs):
+        """Return a copy of *self* with the attributes in *kwargs* replaced."""
+        from dataclasses import replace
+        return replace(self, **kwargs)
+
+
+@with_container_arithmetic(bcast_obj_array=False,
+                           bcast_container_types=(DOFArray, np.ndarray),
+                           matmul=True,
+                           rel_comparison=True)
+@dataclass_array_container
+@dataclass(frozen=True)
+class MixtureConservedVars(ConservedVars):
+    r"""Conserved mixture components for multi-species states.
+
+    .. attribute:: species_mass
+
+        Object array (:class:`~numpy.ndarray`) with shape ``(nspecies,)``
+        of :class:`~meshmode.dof_array.DOFArray`, or an object array with shape
+        ``(nspecies, ndim)`` respectively for scalar or vector quantities
+        corresponding to the `nspecies` species mass conservation equations.
+
+    .. automethod:: join
+    .. automethod:: replace
+    """
+
+    species_mass: np.ndarray = np.empty((0,), dtype=object)  # empty = immutable
+
+    @property
+    def has_multispecies(self):
+        """Return whether this CV has multiple species."""
+        return True
+
     @property
     def nspecies(self):
         """Return the number of mixture species."""
@@ -254,24 +263,14 @@ def species_mass_fractions(self):
         return self.species_mass / self.mass
 
     def join(self):
-        """Call :func:`join_conserved` on *self*."""
-        return join_conserved(
+        """Return a flat object array of CV components."""
+        return _join_conserved(
             dim=self.dim,
             mass=self.mass,
             energy=self.energy,
             momentum=self.momentum,
             species_mass=self.species_mass)
 
-    def __reduce__(self):
-        """Return a tuple reproduction of self for pickling."""
-        return (ConservedVars, tuple(getattr(self, f.name)
-                                    for f in fields(ConservedVars)))
-
-    def replace(self, **kwargs):
-        """Return a copy of *self* with the attributes in *kwargs* replaced."""
-        from dataclasses import replace
-        return replace(self, **kwargs)
-
 
 def _aux_shape(ary, leading_shape):
     """:arg leading_shape: a tuple with which ``ary.shape`` is expected to begin."""
@@ -295,7 +294,7 @@ def get_num_species(dim, q):
     return len(q) - (dim + 2)
 
 
-def split_conserved(dim, q):
+def _split_conserved(dim, q):
     """Get quantities corresponding to fluid conservation equations.
 
     Return a :class:`ConservedVars` with quantities corresponding to the
@@ -306,51 +305,57 @@ def split_conserved(dim, q):
     array.
     """
     nspec = get_num_species(dim, q)
-    return ConservedVars(mass=q[0], energy=q[1], momentum=q[2:2+dim],
-                         species_mass=q[2+dim:2+dim+nspec])
+    if nspec > 0:
+        return MixtureConservedVars(mass=q[0], energy=q[1], momentum=q[2:2+dim],
+                                    species_mass=q[2+dim:2+dim+nspec])
+    return ConservedVars(mass=q[0], energy=q[1], momentum=q[2:2+dim])
 
 
 def _join_conserved(dim, mass, energy, momentum, species_mass=None):
-    if species_mass is None:  # empty: immutable
-        species_mass = np.empty((0,), dtype=object)
-
-    nspec = len(species_mass)
-    aux_shapes = [
-        _aux_shape(mass, ()),
-        _aux_shape(energy, ()),
-        _aux_shape(momentum, (dim,)),
-        _aux_shape(species_mass, (nspec,))]
+    multigas = False
+    if species_mass is not None:
+        nspec = len(species_mass)
+        if nspec > 0:
+            multigas = True
+
+    if multigas:  # empty: immutable
+        aux_shapes = [
+            _aux_shape(mass, ()),
+            _aux_shape(energy, ()),
+            _aux_shape(momentum, (dim,)),
+            _aux_shape(species_mass, (nspec,))]
+    else:
+        nspec = 0
+        aux_shapes = [
+            _aux_shape(mass, ()),
+            _aux_shape(energy, ()),
+            _aux_shape(momentum, (dim,))]
 
     from pytools import single_valued
     aux_shape = single_valued(aux_shapes)
-
     result = np.empty((2+dim+nspec,) + aux_shape, dtype=object)
     result[0] = mass
     result[1] = energy
     result[2:dim+2] = momentum
-    result[dim+2:] = species_mass
-
-    return result
 
+    if multigas:
+        result[dim+2:] = species_mass
 
-def join_conserved(dim, mass, energy, momentum, species_mass=None):
-    """Create agglomerated array from quantities for each conservation eqn."""
-    return _join_conserved(dim, mass=mass, energy=energy,
-                           momentum=momentum, species_mass=species_mass)
+    return result
 
 
 def make_conserved(dim, mass=None, energy=None, momentum=None, species_mass=None,
                    q=None, scalar_quantities=None, vector_quantities=None):
     """Create :class:`ConservedVars` from separated conserved quantities."""
     if scalar_quantities is not None:
-        return split_conserved(dim, q=scalar_quantities)
+        return _split_conserved(dim, q=scalar_quantities)
     if vector_quantities is not None:
-        return split_conserved(dim, q=vector_quantities)
+        return _split_conserved(dim, q=vector_quantities)
     if q is not None:
-        return split_conserved(dim, q=q)
+        return _split_conserved(dim, q=q)
     if mass is None or energy is None or momentum is None:
         raise ValueError("Must have one of *q* or *mass, energy, momentum*.")
-    return split_conserved(
+    return _split_conserved(
         dim, _join_conserved(dim, mass=mass, energy=energy,
                              momentum=momentum, species_mass=species_mass)
     )

mirgecom/initializers.py

@@ -709,11 +709,10 @@ def __init__(
         else:
             self._velocity = np.zeros(shape=(dim,))
 
+        self._nspecies = nspecies
         if mass_fracs is not None:
-            self._nspecies = len(mass_fracs)
             self._mass_fracs = mass_fracs
         else:
-            self._nspecies = nspecies
             self._mass_fracs = np.zeros(shape=(nspecies,))
 
         if self._velocity.shape != (dim,):
@@ -742,7 +741,7 @@ def __call__(self, x_vec, *, eos=None, **kwargs):
         mass = 0.0 * x_vec[0] + self._rho
         mom = self._velocity * mass
         energy = (self._p / (gamma - 1.0)) + np.dot(mom, mom) / (2.0 * mass)
-        species_mass = self._mass_fracs * mass
+        species_mass = self._mass_fracs * mass if self._nspecies > 0 else None
 
         return make_conserved(dim=self._dim, mass=mass, energy=energy,
                               momentum=mom, species_mass=species_mass)
@@ -765,7 +764,9 @@ def exact_rhs(self, discr, cv, time=0.0):
         massrhs = 0.0 * mass
         energyrhs = 0.0 * mass
         momrhs = make_obj_array([0 * mass for i in range(self._dim)])
-        yrhs = make_obj_array([0 * mass for i in range(self._nspecies)])
+        yrhs = None
+        if self._nspecies > 0:
+            yrhs = make_obj_array([0 * mass for i in range(self._nspecies)])
 
         return make_conserved(dim=self._dim, mass=massrhs, energy=energyrhs,
                               momentum=momrhs, species_mass=yrhs)
@@ -968,6 +969,9 @@ def exact_grad(self, x_vec, eos, cv_exact):
         dvy = make_obj_array([0*x, 0*x])
         dv = np.stack((dvx, dvy))
         dmom = mass*dv
-        species_mass = velocity*cv_exact.species_mass.reshape(-1, 1)
+        species_mass = None
+        if cv_exact.has_multispecies:
+            species_mass = velocity*cv_exact.species_mass.reshape(-1, 1)
+
         return make_conserved(2, mass=dmass, energy=denergy,
                               momentum=dmom, species_mass=species_mass)

mirgecom/inviscid.py

@@ -61,14 +61,14 @@ def inviscid_flux(discr, eos, cv):
     """
     dim = cv.dim
     p = eos.pressure(cv)
-
     mom = cv.momentum
+    species_mass = \
+        cv.species_mass.reshape(-1, 1)*mom/cv.mass if cv.has_multispecies else None
 
     return make_conserved(
         dim, mass=mom, energy=mom * (cv.energy + p) / cv.mass,
         momentum=np.outer(mom, mom) / cv.mass + np.eye(dim)*p,
-        species_mass=(  # reshaped: (nspecies, dim)
-            (mom / cv.mass) * cv.species_mass.reshape(-1, 1)))
+        species_mass=species_mass)
 
 
 def inviscid_facial_flux(discr, eos, cv_tpair, local=False):