diff --git a/.ci-support/production-testing-env.sh b/.ci-support/production-testing-env.sh index 130a3c53e..e07e9b4af 100755 --- a/.ci-support/production-testing-env.sh +++ b/.ci-support/production-testing-env.sh @@ -8,7 +8,7 @@ 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 @@ -16,7 +16,7 @@ set -x # "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" diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index a4088fa2a..2dd54d6fb 100644 --- a/mirgecom/boundary.py +++ b/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) diff --git a/mirgecom/fluid.py b/mirgecom/fluid.py index 9399df7db..62dd4f51b 100644 --- a/mirgecom/fluid.py +++ b/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) ) diff --git a/mirgecom/initializers.py b/mirgecom/initializers.py index cc3f7141f..29e6e09be 100644 --- a/mirgecom/initializers.py +++ b/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) diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index c04eaf578..eef2f9c2c 100644 --- a/mirgecom/inviscid.py +++ b/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): diff --git a/mirgecom/viscous.py b/mirgecom/viscous.py index e1383ed99..6da04ae29 100644 --- a/mirgecom/viscous.py +++ b/mirgecom/viscous.py @@ -126,7 +126,6 @@ def diffusive_flux(discr, eos, cv, grad_cv): The species diffusive flux vector, $\mathbf{J}_{\alpha}$ """ transport = eos.transport_model() - grad_y = species_mass_fraction_gradient(discr, cv, grad_cv) d = transport.species_diffusivity(eos, cv) return -cv.mass*d.reshape(-1, 1)*grad_y @@ -253,8 +252,12 @@ def viscous_flux(discr, eos, cv, grad_cv, grad_t): dim = cv.dim viscous_mass_flux = 0 * cv.momentum - j = diffusive_flux(discr, eos, cv, grad_cv) - heat_flux_diffusive = diffusive_heat_flux(discr, eos, cv, j) + heat_flux_diffusive = 0 + species_mass_flux = None + if cv.has_multispecies: + j = diffusive_flux(discr, eos, cv, grad_cv) + heat_flux_diffusive = diffusive_heat_flux(discr, eos, cv, j) + species_mass_flux = -j tau = viscous_stress_tensor(discr, eos, cv, grad_cv) viscous_energy_flux = ( @@ -264,9 +267,9 @@ def viscous_flux(discr, eos, cv, grad_cv, grad_t): ) return make_conserved(dim, - mass=viscous_mass_flux, - energy=viscous_energy_flux, - momentum=tau, species_mass=-j) + mass=viscous_mass_flux, + energy=viscous_energy_flux, + momentum=tau, species_mass=species_mass_flux) def viscous_facial_flux(discr, eos, cv_tpair, grad_cv_tpair, grad_t_tpair, diff --git a/test/test_euler.py b/test/test_euler.py index f93cbe70c..60ee7f99a 100644 --- a/test/test_euler.py +++ b/test/test_euler.py @@ -100,12 +100,14 @@ def test_uniform_rhs(actx_factory, nspecies, dim, order): mom_input = make_obj_array( [discr.zeros(actx) for i in range(discr.dim)] ) + num_equations = dim + 2 + nspecies - mass_frac_input = flat_obj_array( - [ones / ((i + 1) * 10) for i in range(nspecies)] - ) - species_mass_input = mass_input * mass_frac_input - num_equations = dim + 2 + len(species_mass_input) + species_mass_input = None + if nspecies > 0: + mass_frac_input = flat_obj_array( + [ones / ((i + 1) * 10) for i in range(nspecies)] + ) + species_mass_input = mass_input * mass_frac_input cv = make_conserved( dim, mass=mass_input, energy=energy_input, momentum=mom_input, @@ -124,18 +126,15 @@ def test_uniform_rhs(actx_factory, nspecies, dim, order): rho_resid = rhs_resid.mass rhoe_resid = rhs_resid.energy mom_resid = rhs_resid.momentum - rhoy_resid = rhs_resid.species_mass rho_rhs = inviscid_rhs.mass rhoe_rhs = inviscid_rhs.energy rhov_rhs = inviscid_rhs.momentum - rhoy_rhs = inviscid_rhs.species_mass logger.info( f"rho_rhs = {rho_rhs}\n" f"rhoe_rhs = {rhoe_rhs}\n" f"rhov_rhs = {rhov_rhs}\n" - f"rhoy_rhs = {rhoy_rhs}\n" ) def inf_norm(x): @@ -145,8 +144,13 @@ def inf_norm(x): assert inf_norm(rhoe_resid) < tolerance for i in range(dim): assert inf_norm(mom_resid[i]) < tolerance - for i in range(nspecies): - assert inf_norm(rhoy_resid[i]) < tolerance + + if nspecies > 0: + rhoy_resid = rhs_resid.species_mass + rhoy_rhs = inviscid_rhs.species_mass + logger.info(f"rhoy_rhs = {rhoy_rhs}") + for i in range(nspecies): + assert inf_norm(rhoy_resid[i]) < tolerance err_max = inf_norm(rho_resid) eoc_rec0.add_data_point(1.0 / nel_1d, err_max) @@ -167,15 +171,17 @@ def inf_norm(x): rho_resid = rhs_resid.mass rhoe_resid = rhs_resid.energy mom_resid = rhs_resid.momentum - rhoy_resid = rhs_resid.species_mass assert inf_norm(rho_resid) < tolerance assert inf_norm(rhoe_resid) < tolerance for i in range(dim): assert inf_norm(mom_resid[i]) < tolerance - for i in range(nspecies): - assert inf_norm(rhoy_resid[i]) < tolerance + + if nspecies > 0: + rhoy_resid = rhs_resid.species_mass + for i in range(nspecies): + assert inf_norm(rhoy_resid[i]) < tolerance err_max = inf_norm(rho_resid) eoc_rec1.add_data_point(1.0 / nel_1d, err_max) diff --git a/test/test_init.py b/test/test_init.py index c65c445db..f63058113 100644 --- a/test/test_init.py +++ b/test/test_init.py @@ -79,7 +79,8 @@ def test_uniform_init(ctx_factory, dim, nspecies): from mirgecom.initializers import Uniform mass_fracs = np.array([float(ispec+1) for ispec in range(nspecies)]) - initializer = Uniform(dim=dim, mass_fracs=mass_fracs, velocity=velocity) + initializer = Uniform(dim=dim, mass_fracs=mass_fracs, velocity=velocity, + nspecies=nspecies) cv = initializer(nodes) def inf_norm(data): @@ -98,13 +99,14 @@ def inf_norm(data): exp_energy = 2.5 + .5 * dim eerrmax = inf_norm(cv.energy - exp_energy) - exp_species_mass = exp_mass * mass_fracs - mferrmax = inf_norm(cv.species_mass - exp_species_mass) - assert perrmax < 1e-15 assert merrmax < 1e-15 assert eerrmax < 1e-15 - assert mferrmax < 1e-15 + + if nspecies > 0: + exp_species_mass = exp_mass * mass_fracs + mferrmax = inf_norm(cv.species_mass - exp_species_mass) + assert mferrmax < 1e-15 def test_lump_init(ctx_factory): diff --git a/test/test_inviscid.py b/test/test_inviscid.py index c97b37501..76e47ffa6 100644 --- a/test/test_inviscid.py +++ b/test/test_inviscid.py @@ -316,7 +316,8 @@ def inf_norm(data): assert inf_norm(interior_face_flux.mass) < tolerance assert inf_norm(interior_face_flux.energy) < tolerance - assert inf_norm(interior_face_flux.species_mass) < tolerance + if cv.has_multispecies: + assert inf_norm(interior_face_flux.species_mass) < tolerance # The expected pressure is 1.0 (by design). And the flux diagonal is # [rhov_x*v_x + p] (etc) since we have zero velocities it's just p. @@ -354,7 +355,8 @@ def inf_norm(data): assert inf_norm(boundary_flux.mass) < tolerance assert inf_norm(boundary_flux.energy) < tolerance - assert inf_norm(boundary_flux.species_mass) < tolerance + if cv.has_multispecies: + assert inf_norm(boundary_flux.species_mass) < tolerance nhat = thaw(actx, discr.normal(BTAG_ALL)) mom_flux_exact = discr.project(BTAG_ALL, "all_faces", p0*nhat)