add hot mixing layer case

y3prediction/mixing_layer.py

@@ -4,6 +4,144 @@
 from functools import partial
 
 
+class MixingLayerHot:
+    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=2, nspecies=0,
+            inflow_profile
+    ):
+        r"""Initialize mixture parameters.
+
+        Parameters
+        ----------
+        dim: int
+            specifies the number of dimensions for the solution
+        nspecies: int
+            specifies the number of mixture species
+        inflow_profile: dict
+            inflow profile
+        """
+
+        self._dim = dim
+        self._nspecies = nspecies
+
+        self._inflow_y = inflow_profile["y"]
+        self._inflow_rho = inflow_profile["rho"]
+        self._inflow_e = inflow_profile["e_int"]
+        self._inflow_temp = inflow_profile["temp"]
+        self._inflow_vel = inflow_profile["u"]
+        self._inflow_mass_frac = inflow_profile["mass_frac"]
+
+    def __call__(self, dcoll, x_vec, eos, *, time=0.0):
+        """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}.")
+
+        y_offset = 0.
+        ypos = x_vec[1] + y_offset
+        actx = ypos.array_context
+
+        zeros = actx.np.zeros_like(ypos)
+        rho = zeros
+        vmag = zeros
+        velocity = np.zeros(self._dim, dtype=object)*zeros
+        energy = zeros
+        temperature = zeros
+        mf = np.zeros(self._nspecies, dtype=object)*zeros
+
+        rho_bottom = self._inflow_rho[0]
+        vel_bottom = self._inflow_vel[0]
+        temp_bottom = self._inflow_temp[0]
+        e_bottom = self._inflow_e[0]
+        mf_bottom = self._inflow_mass_frac[:, 0]
+        y_bottom = self._inflow_y[0] + y_offset
+
+        # iterate over every interval in the profile
+        # this is expensive for large input data sets
+        for ind in range(1, self._inflow_y.shape[0]):
+
+            rho_top = self._inflow_rho[ind]
+            vel_top = self._inflow_vel[ind]
+            temp_top = self._inflow_temp[ind]
+            e_top = self._inflow_e[ind]
+            mf_top = self._inflow_mass_frac[:, ind]
+
+            # interpolate our data
+            y_top = self._inflow_y[ind] + y_offset
+
+            dy = (y_top - y_bottom)
+            drho = rho_top - rho_bottom
+            dvel = vel_top - vel_bottom
+            dtemp = temp_top - temp_bottom
+            de = e_top - e_bottom
+            dmf = mf_top - mf_bottom
+
+            local_rho = rho_bottom + (ypos - y_bottom)*drho/dy
+            local_vel = vel_bottom + (ypos - y_bottom)*dvel/dy
+            local_temp = temp_bottom + (ypos - y_bottom)*dtemp/dy
+            local_e = e_bottom + (ypos - y_bottom)*de/dy
+            local_mf = mf_bottom + (ypos - y_bottom)*dmf/dy
+
+            # extend just a a little bit to catch the edges
+            bottom_edge = actx.np.greater(ypos, y_bottom - 1.e-12)
+            top_edge = actx.np.less(ypos, y_top + 1.e-12)
+            inside_block = bottom_edge*top_edge
+
+            rho = actx.np.where(inside_block, local_rho, rho)
+            vmag = actx.np.where(inside_block, local_vel, vmag)
+            temperature = actx.np.where(inside_block, local_temp, temperature)
+            energy = actx.np.where(inside_block, local_e, energy)
+            for i in range(self._nspecies):
+                mf[i] = actx.np.where(inside_block, local_mf[i], mf[i])
+
+            y_bottom = y_top
+            rho_bottom = rho_top
+            vel_bottom = vel_top
+            temp_bottom = temp_top
+            e_bottom = e_top
+            mf_bottom = mf_top
+
+        velocity[0] = vmag
+        mom = velocity*rho
+        #internal_energy = eos.get_internal_energy(temperature, mf)
+        internal_energy = energy
+
+        kinetic_energy = 0.5 * np.dot(velocity, velocity)
+        total_energy = rho*(internal_energy + kinetic_energy)
+
+        return make_conserved(dim=self._dim,
+                              mass=rho,
+                              energy=total_energy,
+                              momentum=mom,
+                              species_mass=rho*mf)
+
+
 class MixingLayerCold:
     r"""Solution initializer for flow with a discontinuity.
 

y3prediction/prediction.py

@@ -1627,6 +1627,10 @@ def boundary_type_sanity(boundary, boundary_type):
             print("\tInitializing flow to mixing_layer")
             print(f"Vorticity thickness {vorticity_thickness}")
             print(f"Ambient pressure {pres_bkrnd}")
+        elif init_case == "mixing_layer_hot":
+            print("\tInitializing flow to mixing_layer_hot")
+            print(f"Vorticity thickness {vorticity_thickness}")
+            print(f"Ambient pressure {pres_bkrnd}")
         elif init_case == "flame1d":
             print("\tInitializing flow to flame1d")
             print(f"Ambient pressure {pres_bkrnd}")
@@ -1657,6 +1661,7 @@ def boundary_type_sanity(boundary, boundary_type):
                 "\t flame1d"
                 "\t wedge"
                 "\t mixing_layer"
+                "\t mixing_layer_hot"
                 "\t species_diffusion"
             )
         print("#### Simluation initialization data: ####")
@@ -2216,6 +2221,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs):
             temp_wall=temp_bkrnd,
             vel_sigma=vel_sigma,
             temp_sigma=temp_sigma)
+
     if init_case == "mixing_layer":
         temperature = 300.
         pressure = 101325.
@@ -2239,7 +2245,38 @@ def _compiled_stepper_wrapper(state, t, dt, rhs):
             vorticity_thickness=vorticity_thickness,
             pressure=pres_bkrnd
         )
-    if init_case == "flame1d":
+    if init_case == "mixing_layer_hot":
+        if rank == 0:
+            print("Initializing hot mixing layer")
+
+        import h5py
+
+        def get_data_from_hdf5(group):
+            data_dict = {}
+            for key in group.keys():
+                if isinstance(group[key], h5py.Group):
+                    # If the key is a group, recursively explore it
+                    subgroup_data = get_data_from_hdf5(group[key])
+                    data_dict.update(subgroup_data)
+                elif isinstance(group[key], h5py.Dataset):
+                    # If it's a dataset, add it to the dictionary
+                    data_dict[key] = group[key][()]
+            return data_dict
+
+        # Usage example
+        inflow_fname = "r_mixing_layer_inflow.h5"
+        with h5py.File(inflow_fname, "r") as hf:
+            inflow_data = get_data_from_hdf5(hf)
+
+        #print(f"{inflow_data=}")
+
+        from y3prediction.mixing_layer import MixingLayerHot
+        bulk_init = MixingLayerHot(
+            dim=dim, nspecies=nspecies,
+            inflow_profile=inflow_data
+        )
+
+    elif init_case == "flame1d":
 
         # init params
         disc_location = np.zeros(shape=(dim,))
@@ -3371,7 +3408,7 @@ def get_mesh_data():
                     mesh = rotate_mesh_around_axis(mesh, theta=theta)
 
                 return mesh, tag_to_elements, volume_to_tags
-        elif init_case == "mixing_layer":
+        elif init_case == "mixing_layer" or init_case == "mixing_layer_hot":
             if rank == 0:
                 print("Generating mesh from scratch")