Add disc init to NS

mirgecom/initializers.py

@@ -11,6 +11,7 @@
 .. autoclass:: AcousticPulse
 .. automethod: make_pulse
 .. autoclass:: MixtureInitializer
+.. autoclass:: PlanarDiscontinuity
 .. autoclass:: PlanarPoiseuille
 """
 
@@ -42,6 +43,7 @@
 from pytools.obj_array import make_obj_array
 from meshmode.dof_array import thaw
 from mirgecom.eos import IdealSingleGas
+from numbers import Number
 from mirgecom.fluid import make_conserved
 
 
@@ -853,6 +855,142 @@ def __call__(self, x_vec, eos, **kwargs):
                               momentum=mom, species_mass=specmass)
 
 
+class PlanarDiscontinuity:
+    r"""Solution initializer for flow with a discontinuity.
+
+    This initializer creates a physics-consistent flow solution
+    given an initial thermal state (pressure, temperature) and an EOS.
+
+    The solution varies across a planar interface defined by a tanh function
+    located at disc_location for pressure, temperature, velocity, and mass fraction
+
+    .. automethod:: __init__
+    .. automethod:: __call__
+    """
+
+    def __init__(
+            self, *, dim=3, normal_dir=0, disc_location=0, nspecies=0,
+            temperature_left, temperature_right,
+            pressure_left, pressure_right,
+            velocity_left=None, velocity_right=None,
+            species_mass_left=None, species_mass_right=None,
+            convective_velocity=None, sigma=0.5
+    ):
+        r"""Initialize mixture parameters.
+
+        Parameters
+        ----------
+        dim: int
+            specifies the number of dimensions for the solution
+        normal_dir: int
+            specifies the direction (plane) the discontinuity is applied in
+        disc_location: float or Callable
+            fixed location of discontinuity or optionally a function that
+            returns the time-dependent location.
+        nspecies: int
+            specifies the number of mixture species
+        pressure_left: float
+            pressure to the left of the discontinuity
+        temperature_left: float
+            temperature to the left of the discontinuity
+        velocity_left: numpy.ndarray
+            velocity (vector) to the left of the discontinuity
+        species_mass_left: numpy.ndarray
+            species mass fractions to the left of the discontinuity
+        pressure_right: float
+            pressure to the right of the discontinuity
+        temperature_right: float
+            temperaure to the right of the discontinuity
+        velocity_right: numpy.ndarray
+            velocity (vector) to the right of the discontinuity
+        species_mass_right: numpy.ndarray
+            species mass fractions to the right of the discontinuity
+        sigma: float
+           sharpness parameter
+        """
+        if velocity_left is None:
+            velocity_left = np.zeros(shape=(dim,))
+        if velocity_right is None:
+            velocity_right = np.zeros(shape=(dim,))
+
+        if species_mass_left is None:
+            species_mass_left = np.zeros(shape=(nspecies,))
+        if species_mass_right is None:
+            species_mass_right = np.zeros(shape=(nspecies,))
+
+        self._nspecies = nspecies
+        self._dim = dim
+        self._disc_location = disc_location
+        self._sigma = sigma
+        self._ul = velocity_left
+        self._ur = velocity_right
+        self._uc = convective_velocity
+        self._pl = pressure_left
+        self._pr = pressure_right
+        self._tl = temperature_left
+        self._tr = temperature_right
+        self._yl = species_mass_left
+        self._yr = species_mass_right
+        self._xdir = normal_dir
+        if self._xdir >= self._dim:
+            self._xdir = self._dim - 1
+
+    def __call__(self, x_vec, eos, *, time=0.0, **kwargs):
+        """Create the mixture state at locations *x_vec*.
+
+        Parameters
+        ----------
+        x_vec: numpy.ndarray
+            Coordinates at which solution is desired
+        eos:
+            Mixture-compatible equation-of-state object must provide
+            these functions:
+            `eos.get_density`
+            `eos.get_internal_energy`
+        time: float
+            Time at which solution is desired. The location is (optionally)
+            dependent on time
+        """
+        if x_vec.shape != (self._dim,):
+            raise ValueError(f"Position vector has unexpected dimensionality,"
+                             f" expected {self._dim}.")
+
+        x = x_vec[self._xdir]
+        actx = x.array_context
+        if isinstance(self._disc_location, Number):
+            x0 = self._disc_location
+        else:
+            x0 = self._disc_location(time)
+
+        xtanh = 1.0/self._sigma*(x0 - x)
+        weight = 0.5*(1.0 - actx.np.tanh(xtanh))
+        pressure = self._pl + (self._pr - self._pl)*weight
+        temperature = self._tl + (self._tr - self._tl)*weight
+        velocity = self._ul + (self._ur - self._ul)*weight
+        y = self._yl + (self._yr - self._yl)*weight
+
+        if self._nspecies:
+            mass = eos.get_density(pressure, temperature,
+                                   species_mass_fractions=y)
+        else:
+            mass = pressure/temperature/eos.gas_const()
+
+        specmass = mass * y
+        mom = mass * velocity
+        if self._nspecies:
+            internal_energy = \
+                eos.get_internal_energy(temperature,
+                                        species_mass_fractions=y)
+        else:
+            internal_energy = pressure/mass/(eos.gamma() - 1)
+
+        kinetic_energy = 0.5 * np.dot(velocity, velocity)
+        energy = mass * (internal_energy + kinetic_energy)
+
+        return make_conserved(dim=self._dim, mass=mass, energy=energy,
+                              momentum=mom, species_mass=specmass)
+
+
 class PlanarPoiseuille:
     r"""Initializer for the planar Poiseuille case.