From 09a1a4cb08a282cb927719e16eba989a5623c1ae Mon Sep 17 00:00:00 2001
From: anderson2981 <41346365+anderson2981@users.noreply.github.com>
Date: Thu, 9 Nov 2023 16:11:29 -0600
Subject: [PATCH 1/2] Add av diff (#978)

* modify transport

* add support for av=3, diffusivity

* busy work

* stuff

* removed density from artificial diffusivity

---------

Co-authored-by: Mike Campbell <mtcampbe@illinois.edu>
---
 examples/ablation-workshop.py |   4 +-
 mirgecom/eos.py               |  10 +++
 mirgecom/gas_model.py         |  44 +++++++++++-
 mirgecom/transport.py         | 132 ++++++++++++++++++++++++++++++++++
 4 files changed, 188 insertions(+), 2 deletions(-)

diff --git a/examples/ablation-workshop.py b/examples/ablation-workshop.py
index bce9a1361..6668745a3 100644
--- a/examples/ablation-workshop.py
+++ b/examples/ablation-workshop.py
@@ -488,7 +488,8 @@ def get_species_source_terms(self, cv: ConservedVars, temperature: DOFArray):
         raise NotImplementedError
 
     def dependent_vars(self, cv: ConservedVars, temperature_seed=None,
-            smoothness_mu=None, smoothness_kappa=None, smoothness_beta=None):
+            smoothness_mu=None, smoothness_kappa=None,
+            smoothness_d=None, smoothness_beta=None):
         raise NotImplementedError
 
 
@@ -755,6 +756,7 @@ def make_state(cv, temperature_seed, material_densities):
             smoothness_mu=zeros,
             smoothness_kappa=zeros,
             smoothness_beta=zeros,
+            smoothness_d=zeros,
             species_enthalpies=cv.species_mass,  # empty array
         )
 
diff --git a/mirgecom/eos.py b/mirgecom/eos.py
index 149e19687..d7de23cee 100644
--- a/mirgecom/eos.py
+++ b/mirgecom/eos.py
@@ -77,6 +77,7 @@ class GasDependentVars:
     .. attribute:: smoothness_mu
     .. attribute:: smoothness_kappa
     .. attribute:: smoothness_beta
+    .. attribute:: smoothness_d
     """
 
     temperature: DOFArray
@@ -84,6 +85,7 @@ class GasDependentVars:
     speed_of_sound: DOFArray
     smoothness_mu: DOFArray
     smoothness_kappa: DOFArray
+    smoothness_d: DOFArray
     smoothness_beta: DOFArray
 
 
@@ -179,6 +181,7 @@ def dependent_vars(
             temperature_seed: Optional[DOFArray] = None,
             smoothness_mu: Optional[DOFArray] = None,
             smoothness_kappa: Optional[DOFArray] = None,
+            smoothness_d: Optional[DOFArray] = None,
             smoothness_beta: Optional[DOFArray] = None) -> GasDependentVars:
         """Get an agglomerated array of the dependent variables.
 
@@ -194,6 +197,8 @@ def dependent_vars(
             smoothness_mu = zeros
         if smoothness_kappa is None:
             smoothness_kappa = zeros
+        if smoothness_d is None:
+            smoothness_d = zeros
         if smoothness_beta is None:
             smoothness_beta = zeros
 
@@ -203,6 +208,7 @@ def dependent_vars(
             speed_of_sound=self.sound_speed(cv, temperature),
             smoothness_mu=smoothness_mu,
             smoothness_kappa=smoothness_kappa,
+            smoothness_d=smoothness_d,
             smoothness_beta=smoothness_beta
         )
 
@@ -257,6 +263,7 @@ def dependent_vars(
             temperature_seed: Optional[DOFArray] = None,
             smoothness_mu: Optional[DOFArray] = None,
             smoothness_kappa: Optional[DOFArray] = None,
+            smoothness_d: Optional[DOFArray] = None,
             smoothness_beta: Optional[DOFArray] = None) -> MixtureDependentVars:
         """Get an agglomerated array of the dependent variables.
 
@@ -272,6 +279,8 @@ def dependent_vars(
             smoothness_mu = zeros
         if smoothness_kappa is None:
             smoothness_kappa = zeros
+        if smoothness_d is None:
+            smoothness_d = zeros
         if smoothness_beta is None:
             smoothness_beta = zeros
 
@@ -282,6 +291,7 @@ def dependent_vars(
             species_enthalpies=self.species_enthalpies(cv, temperature),
             smoothness_mu=smoothness_mu,
             smoothness_kappa=smoothness_kappa,
+            smoothness_d=smoothness_d,
             smoothness_beta=smoothness_beta
         )
 
diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py
index 4791270a7..46563d2bc 100644
--- a/mirgecom/gas_model.py
+++ b/mirgecom/gas_model.py
@@ -118,6 +118,7 @@ class FluidState:
     .. autoattribute:: temperature
     .. autoattribute:: smoothness_mu
     .. autoattribute:: smoothness_kappa
+    .. autoattribute:: smoothness_d
     .. autoattribute:: smoothness_beta
     .. autoattribute:: velocity
     .. autoattribute:: speed
@@ -169,6 +170,11 @@ def smoothness_kappa(self):
         """Return the smoothness_kappa field."""
         return self.dv.smoothness_kappa
 
+    @property
+    def smoothness_d(self):
+        """Return the smoothness_d field."""
+        return self.dv.smoothness_d
+
     @property
     def smoothness_beta(self):
         """Return the smoothness_beta field."""
@@ -297,6 +303,7 @@ def make_fluid_state(cv, gas_model,
                      temperature_seed=None,
                      smoothness_mu=None,
                      smoothness_kappa=None,
+                     smoothness_d=None,
                      smoothness_beta=None,
                      material_densities=None,
                      limiter_func=None, limiter_dd=None):
@@ -328,6 +335,11 @@ def make_fluid_state(cv, gas_model,
         Optional array containing the smoothness parameter for extra thermal
         conductivity in the artificial viscosity.
 
+    smoothness_d: :class:`~meshmode.dof_array.DOFArray`
+
+        Optional array containing the smoothness parameter for extra species
+        diffusivity in the artificial viscosity.
+
     smoothness_beta: :class:`~meshmode.dof_array.DOFArray`
 
         Optional array containing the smoothness parameter for extra bulk
@@ -357,6 +369,8 @@ def make_fluid_state(cv, gas_model,
                         is None else smoothness_kappa)
     smoothness_beta = (actx.np.zeros_like(cv.mass) if smoothness_beta
                        is None else smoothness_beta)
+    smoothness_d = (actx.np.zeros_like(cv.mass) if smoothness_d
+                        is None else smoothness_d)
 
     if isinstance(gas_model, GasModel):
         temperature = gas_model.eos.temperature(cv=cv,
@@ -373,6 +387,7 @@ def make_fluid_state(cv, gas_model,
             speed_of_sound=gas_model.eos.sound_speed(cv, temperature),
             smoothness_mu=smoothness_mu,
             smoothness_kappa=smoothness_kappa,
+            smoothness_d=smoothness_d,
             smoothness_beta=smoothness_beta
         )
 
@@ -384,6 +399,7 @@ def make_fluid_state(cv, gas_model,
                 speed_of_sound=dv.speed_of_sound,
                 smoothness_mu=dv.smoothness_mu,
                 smoothness_kappa=dv.smoothness_kappa,
+                smoothness_d=dv.smoothness_d,
                 smoothness_beta=dv.smoothness_beta,
                 species_enthalpies=gas_model.eos.species_enthalpies(cv, temperature)
             )
@@ -428,6 +444,7 @@ def make_fluid_state(cv, gas_model,
             speed_of_sound=gas_model.eos.sound_speed(cv, temperature),
             smoothness_mu=smoothness_mu,
             smoothness_kappa=smoothness_kappa,
+            smoothness_d=smoothness_d,
             smoothness_beta=smoothness_beta,
             species_enthalpies=gas_model.eos.species_enthalpies(cv, temperature),
         )
@@ -507,6 +524,10 @@ def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None,
     if state.dv.smoothness_kappa is not None:
         smoothness_kappa = op.project(dcoll, src, tgt, state.dv.smoothness_kappa)
 
+    smoothness_d = None
+    if state.dv.smoothness_d is not None:
+        smoothness_d = op.project(dcoll, src, tgt, state.dv.smoothness_d)
+
     smoothness_beta = None
     if state.dv.smoothness_beta is not None:
         smoothness_beta = op.project(dcoll, src, tgt, state.dv.smoothness_beta)
@@ -519,6 +540,7 @@ def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None,
                             temperature_seed=temperature_seed,
                             smoothness_mu=smoothness_mu,
                             smoothness_kappa=smoothness_kappa,
+                            smoothness_d=smoothness_d,
                             smoothness_beta=smoothness_beta,
                             material_densities=material_densities,
                             limiter_func=limiter_func, limiter_dd=tgt)
@@ -535,6 +557,7 @@ def make_fluid_state_trace_pairs(cv_pairs, gas_model,
                                  temperature_seed_pairs=None,
                                  smoothness_mu_pairs=None,
                                  smoothness_kappa_pairs=None,
+                                 smoothness_d_pairs=None,
                                  smoothness_beta_pairs=None,
                                  material_densities_pairs=None,
                                  limiter_func=None):
@@ -579,6 +602,8 @@ def make_fluid_state_trace_pairs(cv_pairs, gas_model,
         smoothness_mu_pairs = [None] * len(cv_pairs)
     if smoothness_kappa_pairs is None:
         smoothness_kappa_pairs = [None] * len(cv_pairs)
+    if smoothness_d_pairs is None:
+        smoothness_d_pairs = [None] * len(cv_pairs)
     if smoothness_beta_pairs is None:
         smoothness_beta_pairs = [None] * len(cv_pairs)
     if material_densities_pairs is None:
@@ -590,6 +615,7 @@ def make_fluid_state_trace_pairs(cv_pairs, gas_model,
             temperature_seed=_getattr_ish(tseed_pair, "int"),
             smoothness_mu=_getattr_ish(smoothness_mu_pair, "int"),
             smoothness_kappa=_getattr_ish(smoothness_kappa_pair, "int"),
+            smoothness_d=_getattr_ish(smoothness_d_pair, "int"),
             smoothness_beta=_getattr_ish(smoothness_beta_pair, "int"),
             material_densities=_getattr_ish(material_densities_pair, "int"),
             limiter_func=limiter_func, limiter_dd=cv_pair.dd),
@@ -598,6 +624,7 @@ def make_fluid_state_trace_pairs(cv_pairs, gas_model,
             temperature_seed=_getattr_ish(tseed_pair, "ext"),
             smoothness_mu=_getattr_ish(smoothness_mu_pair, "ext"),
             smoothness_kappa=_getattr_ish(smoothness_kappa_pair, "ext"),
+            smoothness_d=_getattr_ish(smoothness_d_pair, "ext"),
             smoothness_beta=_getattr_ish(smoothness_beta_pair, "ext"),
             material_densities=_getattr_ish(material_densities_pair, "ext"),
             limiter_func=limiter_func, limiter_dd=cv_pair.dd))
@@ -605,10 +632,11 @@ def make_fluid_state_trace_pairs(cv_pairs, gas_model,
             tseed_pair,
             smoothness_mu_pair,
             smoothness_kappa_pair,
+            smoothness_d_pair,
             smoothness_beta_pair,
             material_densities_pair in zip(
                 cv_pairs, temperature_seed_pairs,
-                smoothness_mu_pairs, smoothness_kappa_pairs,
+                smoothness_mu_pairs, smoothness_kappa_pairs, smoothness_d_pairs,
                 smoothness_beta_pairs, material_densities_pairs)]
 
 
@@ -628,6 +656,10 @@ class _FluidSmoothnessKappaTag:
     pass
 
 
+class _FluidSmoothnessDiffTag:
+    pass
+
+
 class _FluidSmoothnessBetaTag:
     pass
 
@@ -758,6 +790,14 @@ def make_operator_fluid_states(
                 dcoll, volume_state.smoothness_kappa, volume_dd=dd_vol,
                 tag=(_FluidSmoothnessKappaTag, comm_tag))]
 
+    smoothness_d_interior_pairs = None
+    if volume_state.smoothness_d is not None:
+        smoothness_d_interior_pairs = [
+            interp_to_surf_quad(tpair=tpair)
+            for tpair in interior_trace_pairs(
+                dcoll, volume_state.smoothness_d, volume_dd=dd_vol,
+                tag=(_FluidSmoothnessDiffTag, comm_tag))]
+
     smoothness_beta_interior_pairs = None
     if volume_state.smoothness_beta is not None:
         smoothness_beta_interior_pairs = [
@@ -780,6 +820,7 @@ def make_operator_fluid_states(
         temperature_seed_pairs=tseed_interior_pairs,
         smoothness_mu_pairs=smoothness_mu_interior_pairs,
         smoothness_kappa_pairs=smoothness_kappa_interior_pairs,
+        smoothness_d_pairs=smoothness_d_interior_pairs,
         smoothness_beta_pairs=smoothness_beta_interior_pairs,
         material_densities_pairs=material_densities_interior_pairs,
         limiter_func=limiter_func)
@@ -874,6 +915,7 @@ def replace_fluid_state(
         temperature_seed=new_tseed,
         smoothness_mu=state.smoothness_mu,
         smoothness_kappa=state.smoothness_kappa,
+        smoothness_d=state.smoothness_d,
         smoothness_beta=state.smoothness_beta,
         material_densities=material_densities,
         limiter_func=limiter_func,
diff --git a/mirgecom/transport.py b/mirgecom/transport.py
index 4f6e5bf6f..9260488fb 100644
--- a/mirgecom/transport.py
+++ b/mirgecom/transport.py
@@ -16,6 +16,7 @@
 .. autoclass:: MixtureAveragedTransport
 .. autoclass:: ArtificialViscosityTransportDiv
 .. autoclass:: ArtificialViscosityTransportDiv2
+.. autoclass:: ArtificialViscosityTransportDiv3
 
 Exceptions
 ^^^^^^^^^^
@@ -698,3 +699,134 @@ def species_diffusivity(self, cv: ConservedVars,
                             eos: Optional[GasEOS] = None) -> DOFArray:
         r"""Get the vector of species diffusivities, ${d}_{\alpha}$."""
         return self._physical_transport.species_diffusivity(cv, dv, eos)
+
+
+class ArtificialViscosityTransportDiv3(TransportModel):
+    r"""Transport model for add artificial viscosity.
+
+    Inherits from (and implements) :class:`TransportModel`.
+
+    Takes a physical transport model and adds the artificial viscosity
+    contribution to it. Defaults to simple transport with inviscid settings.
+    This is equivalent to inviscid flow with artifical viscosity enabled.
+
+    .. automethod:: __init__
+    .. automethod:: bulk_viscosity
+    .. automethod:: viscosity
+    .. automethod:: volume_viscosity
+    .. automethod:: thermal_conductivity
+    .. automethod:: species_diffusivity
+    """
+
+    def __init__(self,
+                 av_mu, av_kappa, av_beta, av_d, av_prandtl,
+                 physical_transport=None,
+                 av_species_diffusivity=None):
+        """Initialize uniform, constant transport properties."""
+        if physical_transport is None:
+            self._physical_transport = SimpleTransport()
+        else:
+            self._physical_transport = physical_transport
+
+        if av_species_diffusivity is None:
+            av_species_diffusivity = np.empty((0,), dtype=object)
+
+        self._av_mu = av_mu
+        self._av_beta = av_beta
+        self._av_kappa = av_kappa
+        self._av_d = av_d
+        self._av_prandtl = av_prandtl
+
+    def av_mu(self, cv, dv, eos):
+        r"""Get the shear artificial viscosity for the gas."""
+        actx = cv.array_context
+        return (self._av_mu * cv.mass
+                * actx.np.sqrt(np.dot(cv.velocity, cv.velocity)
+                               + dv.speed_of_sound**2))
+
+    def av_beta(self, cv, dv, eos):
+        r"""Get the shear artificial viscosity for the gas."""
+        actx = cv.array_context
+        return (self._av_beta * cv.mass
+                * actx.np.sqrt(np.dot(cv.velocity, cv.velocity)
+                               + dv.speed_of_sound**2))
+
+    def av_kappa(self, cv, dv, eos):
+        r"""Get the shear artificial viscosity for the gas."""
+        actx = cv.array_context
+        return (self._av_kappa * cv.mass
+                * actx.np.sqrt(np.dot(cv.velocity, cv.velocity)
+                               + dv.speed_of_sound**2))
+
+    def av_d(self, cv, dv, eos):
+        r"""Get the shear artificial viscosity for the gas."""
+        actx = cv.array_context
+        return (self._av_d * actx.np.sqrt(np.dot(cv.velocity, cv.velocity)
+                               + dv.speed_of_sound**2))
+
+    def bulk_viscosity(self, cv: ConservedVars,  # type: ignore[override]
+                       dv: GasDependentVars,
+                       eos: GasEOS) -> DOFArray:
+        r"""Get the bulk viscosity for the gas, $\mu_{B}$."""
+        actx = cv.array_context
+        smoothness_beta = actx.np.where(
+            actx.np.greater(dv.smoothness_beta, 0.), dv.smoothness_beta, 0.)
+        return (smoothness_beta*self.av_beta(cv, dv, eos)
+                + self._physical_transport.bulk_viscosity(cv, dv))
+
+    def viscosity(self, cv: ConservedVars,  # type: ignore[override]
+                  dv: GasDependentVars,
+                  eos: GasEOS) -> DOFArray:
+        r"""Get the gas dynamic viscosity, $\mu$."""
+        actx = cv.array_context
+        smoothness_mu = actx.np.where(
+            actx.np.greater(dv.smoothness_mu, 0.), dv.smoothness_mu, 0.)
+        return (smoothness_mu*self.av_mu(cv, dv, eos)
+                + self._physical_transport.viscosity(cv, dv))
+
+    def volume_viscosity(self, cv: ConservedVars,  # type: ignore[override]
+                         dv: GasDependentVars,
+                         eos: GasEOS) -> DOFArray:
+        r"""Get the 2nd viscosity coefficent, $\lambda$.
+
+        In this transport model, the second coefficient of viscosity is:
+
+        $\lambda = \left(\mu_{B} - \frac{2\mu}{3}\right)$
+        """
+        actx = cv.array_context
+        smoothness_mu = actx.np.where(
+            actx.np.greater(dv.smoothness_mu, 0.), dv.smoothness_mu, 0.)
+        smoothness_beta = actx.np.where(
+            actx.np.greater(dv.smoothness_beta, 0.), dv.smoothness_beta, 0.)
+
+        return (smoothness_beta*self.av_beta(cv, dv, eos)
+                - 2*smoothness_mu*self.av_mu(cv, dv, eos)/3
+                + self._physical_transport.volume_viscosity(cv, dv))
+
+    def thermal_conductivity(self, cv: ConservedVars,  # type: ignore[override]
+                             dv: GasDependentVars,
+                             eos: GasEOS) -> DOFArray:
+        r"""Get the gas thermal_conductivity, $\kappa$."""
+        cp = eos.heat_capacity_cp(cv, dv.temperature)
+        actx = cv.array_context
+        smoothness_beta = actx.np.where(
+            actx.np.greater(dv.smoothness_beta, 0.), dv.smoothness_beta, 0.)
+        smoothness_kappa = actx.np.where(
+            actx.np.greater(dv.smoothness_kappa, 0.), dv.smoothness_kappa, 0.)
+        av_kappa = (
+            cp*(smoothness_beta*self.av_beta(cv, dv, eos)/self._av_prandtl
+                + smoothness_kappa*self.av_kappa(cv, dv, eos))
+        )
+
+        return (av_kappa
+                + self._physical_transport.thermal_conductivity(cv, dv, eos))
+
+    def species_diffusivity(self, cv: ConservedVars,  # type: ignore[override]
+                            dv: GasDependentVars,
+                            eos: GasEOS) -> DOFArray:
+        r"""Get the vector of species diffusivities, ${d}_{\alpha}$."""
+        actx = cv.array_context
+        smoothness_d = actx.np.where(
+            actx.np.greater(dv.smoothness_d, 0.), dv.smoothness_d, 0.)
+        return (smoothness_d*self.av_d(cv, dv, eos)
+                + self._physical_transport.species_diffusivity(cv, dv, eos))

From d70ac5e23184ad770b18832e881e4a153031ae6b Mon Sep 17 00:00:00 2001
From: Matthias Diener <mdiener@illinois.edu>
Date: Mon, 13 Nov 2023 15:05:22 -0600
Subject: [PATCH 2/2] Clean up temp files on porter before CI runs (#980)

---
 .github/workflows/ci.yaml | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index c2d570390..aa0fcc7bf 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -119,12 +119,16 @@ jobs:
 
         - name: Run examples
           run: |
+            set -x
             MINIFORGE_INSTALL_DIR=.miniforge3
             . "$MINIFORGE_INSTALL_DIR/bin/activate" testing
             export XDG_CACHE_HOME=/tmp
             mamba install vtk  # needed for the accuracy comparison
             [[ $(hostname) == "porter" ]] && export PYOPENCL_TEST="port:nv" && unset XDG_CACHE_HOME
-            # && export POCL_DEBUG=cuda
+
+             # This is only possible because actions run sequentially on porter
+            [[ $(hostname) == "porter" ]] && rm -rf /tmp/githubrunner/pocl-scratch && rm -rf /tmp/githubrunner/xdg-scratch
+
             scripts/run-integrated-tests.sh --examples
 
     doc:
@@ -195,8 +199,13 @@ jobs:
             fetch-depth: '0'
         - name: Prepare production environment
           run: |
+            set -x
             [[ $(uname) == Linux ]] && [[ $(hostname) != "porter" ]] && sudo apt-get update && sudo apt-get install -y openmpi-bin libopenmpi-dev
             [[ $(uname) == Darwin ]] && brew upgrade && brew install mpich
+
+            # This is only possible because actions run sequentially on porter
+            [[ $(hostname) == "porter" ]] && rm -rf /tmp/githubrunner/pocl-scratch && rm -rf /tmp/githubrunner/xdg-scratch
+
             MIRGEDIR=$(pwd)
             cat scripts/production-testing-env.sh
             . scripts/production-testing-env.sh