From a26d50e6d576bef6a4841ce3a34ab16fb9bb2af7 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:19:42 -0500
Subject: [PATCH 01/48] Fix AW example

---
 examples/ablation-workshop.py | 59 ++++++++++++++++++-----------------
 1 file changed, 30 insertions(+), 29 deletions(-)

diff --git a/examples/ablation-workshop.py b/examples/ablation-workshop.py
index 8fb990169..6f44f3602 100644
--- a/examples/ablation-workshop.py
+++ b/examples/ablation-workshop.py
@@ -77,6 +77,7 @@
     PorousWallTransport
 )
 from mirgecom.fluid import ConservedVars
+from mirgecom.materials.tacot import TacotEOS as OriginalTacotEOS
 from logpyle import IntervalTimer, set_dt
 from typing import Optional, Union
 from pytools.obj_array import make_obj_array
@@ -339,6 +340,7 @@ def thermal_conductivity(self, cv: ConservedVars,  # type: ignore[override]
 
     def species_diffusivity(self, cv: ConservedVars,  # type: ignore[override]
             dv: GasDependentVars, eos: GasEOS) -> DOFArray:
+        """Return the (empty) species mass diffusivities."""
         return cv.species_mass  # empty array
 
 
@@ -347,9 +349,8 @@ class TabulatedGasEOS(MixtureEOS):
 
     This section is to be used when species conservation is not employed and
     the output gas is assumed to be in chemical equilibrium.
-    The table was extracted from the suplementar material from the
-    ablation workshop. Some lines were removed to reduce the number of spline
-    interpolation segments.
+    The table was extracted from the ablation workshop suplementary material.
+    Some lines were removed to reduce the number of spline interpolation segments.
     """
 
     def __init__(self):
@@ -448,6 +449,7 @@ def pressure(self, cv: ConservedVars, temperature: DOFArray) -> DOFArray:
     def gas_const(self, cv: Optional[ConservedVars] = None,
                   temperature: Optional[DOFArray] = None,
                   species_mass_fractions: Optional[np.ndarray] = None) -> DOFArray:
+        """Return the specific gas constant."""
         coeffs = self._cs_molar_mass.c
         bnds = self._cs_molar_mass.x
         molar_mass = eval_spline(temperature, bnds, coeffs)
@@ -565,6 +567,16 @@ def pressure_diffusivity(self, cv: ConservedVars, wv: PorousWallVars,
         return cv.mass*wv.permeability/(viscosity*wv.void_fraction)
 
 
+class TacotEOS(OriginalTacotEOS):
+    """Inherits and modified the original TACOT material."""
+
+    def permeability(self, tau: DOFArray) -> DOFArray:
+        r"""Permeability $K$ of the composite material."""
+        virgin = 1.6e-11
+        char = 2.0e-11
+        return virgin*tau + char*(1.0 - tau)
+
+
 def binary_sum(ary):
     """Sum the elements of an array, creating a log-depth DAG instead of linear."""
     n = len(ary)
@@ -592,13 +604,13 @@ def main(actx_class=None, use_logmgr=True, casename=None, restart_file=None):
     viz_path = "viz_data/"
     vizname = viz_path+casename
 
-    t_final = 4.0e-6
+    t_final = 2.0e-7
 
     dim = 1
 
-    order = 2
-    dt = 4.0e-8
-    pressure_scaling_factor = 0.1  # noqa N806
+    order = 3
+    dt = 2.0e-8
+    pressure_scaling_factor = 1.0  # noqa N806
 
     nviz = 200
     ngarbage = 50
@@ -622,7 +634,7 @@ def main(actx_class=None, use_logmgr=True, casename=None, restart_file=None):
 
     else:  # generate the grid from scratch
         from functools import partial
-        nel_1d = 201
+        nel_1d = 121
 
         from meshmode.mesh.generation import generate_regular_rect_mesh
         generate_mesh = partial(generate_regular_rect_mesh,
@@ -653,11 +665,8 @@ def main(actx_class=None, use_logmgr=True, casename=None, restart_file=None):
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     pressure_boundaries = {
-        BoundaryDomainTag("prescribed"):
-            DirichletDiffusionBoundary(101325),
-        BoundaryDomainTag("neumann"):
-            NeumannDiffusionBoundary(0.0)
-    }
+        BoundaryDomainTag("prescribed"): DirichletDiffusionBoundary(101325.0),
+        BoundaryDomainTag("neumann"): NeumannDiffusionBoundary(0.0)}
 
     def my_presc_bdry(u_minus):
         return +u_minus
@@ -665,12 +674,8 @@ def my_presc_bdry(u_minus):
     def my_wall_bdry(u_minus):
         return -u_minus
 
-    velocity_boundaries = {
-        BoundaryDomainTag("prescribed"):
-            my_presc_bdry,
-        BoundaryDomainTag("neumann"):
-            my_wall_bdry
-    }
+    velocity_boundaries = {BoundaryDomainTag("prescribed"): my_presc_bdry,
+                           BoundaryDomainTag("neumann"): my_wall_bdry}
 
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -680,7 +685,7 @@ def my_wall_bdry(u_minus):
 
     my_gas = TabulatedGasEOS()
     bprime_class = BprimeTable()
-    my_material = my_composite.TacotEOS(char_mass=220.0, virgin_mass=280.0)
+    my_material = TacotEOS(char_mass=220.0, virgin_mass=280.0)
     pyrolysis = my_composite.Pyrolysis()
 
     base_transport = GasTabulatedTransport()
@@ -740,7 +745,6 @@ def make_state(cv, temperature_seed, material_densities):
             tau=tau,
             density=gas_model.solid_density(material_densities),
             void_fraction=gas_model.wall_eos.void_fraction(tau=tau),
-            emissivity=gas_model.wall_eos.emissivity(tau=tau),
             permeability=gas_model.wall_eos.permeability(tau=tau),
             tortuosity=gas_model.wall_eos.tortuosity(tau=tau)
         )
@@ -895,10 +899,8 @@ def ablation_workshop_flux(dcoll, state, gas_model, velocity, bprime_class,
         dd_bdry_quad = dd_vol_quad.with_domain_tag(bdtag)
 
         temperature_bc = op.project(dcoll, dd_wall, dd_bdry_quad, dv.temperature)
-        m_dot_g = op.project(dcoll, dd_wall, dd_bdry_quad, cv.mass*velocity)
-        emissivity = op.project(dcoll, dd_wall, dd_bdry_quad, wv.emissivity)
-
-        m_dot_g = np.dot(m_dot_g, normal_vec)
+        momentum_bc = op.project(dcoll, dd_wall, dd_bdry_quad, cv.mass*velocity)
+        m_dot_g = np.dot(momentum_bc, normal_vec)
 
         # time-dependent function
         weight = actx.np.where(actx.np.less(time, 0.1), (time/0.1)+1e-7, 1.0)
@@ -949,10 +951,12 @@ def ablation_workshop_flux(dcoll, state, gas_model, velocity, bprime_class,
 
         flux = -(conv_coeff*(h_e - h_w) - m_dot*h_w + m_dot_g*h_g)
 
+        tau_bc = op.project(dcoll, dd_wall, dd_bdry_quad, wv.tau)
+        emissivity = my_material.emissivity(tau=tau_bc)
         radiation = emissivity*5.67e-8*(temperature_bc**4 - 300**4)
 
         # this is the physical flux normal to the boundary
-        return flux - radiation
+        return flux + radiation
 
     def phenolics_operator(dcoll, fluid_state, boundaries, gas_model, pyrolysis,
                            quadrature_tag, dd_wall=DD_VOLUME_ALL, time=0.0,
@@ -1054,7 +1058,6 @@ def my_write_viz(step, t, state):
                       ("WV_phase_3", wv.material_densities[2]),
                       ("WV_tau", wv.tau),
                       ("WV_void_fraction", wv.void_fraction),
-                      ("WV_emissivity", wv.emissivity),
                       ("WV_permeability", wv.permeability),
                       ("WV_tortuosity", wv.tortuosity),
                       ("WV_density", wv.density),
@@ -1106,8 +1109,6 @@ def my_pre_step(step, t, dt, state):
 
             if do_garbage:
                 with gc_timer:
-                    logger.info(
-                        "Running gc.collect() to work around memory growth issue ")
                     gc.collect()
 
             if do_health:

From 2759befa9969b6e35e9d6d4dcd2baed0ed8caa70 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:20:56 -0500
Subject: [PATCH 02/48] Remove emissivity from PorousVars

---
 mirgecom/gas_model.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py
index 068608de2..da8157ee5 100644
--- a/mirgecom/gas_model.py
+++ b/mirgecom/gas_model.py
@@ -424,8 +424,6 @@ def make_fluid_state(cv, gas_model,
             tau=tau,
             density=gas_model.solid_density(material_densities),
             void_fraction=gas_model.wall_eos.void_fraction(tau=tau),
-            # FIXME emissivity as a function of temperature...
-            emissivity=gas_model.wall_eos.emissivity(tau=tau),
             permeability=gas_model.wall_eos.permeability(tau=tau),
             tortuosity=gas_model.wall_eos.tortuosity(tau=tau)
         )

From c08165a0e46f08e7f739282fa14be04c4bb18aaf Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:21:13 -0500
Subject: [PATCH 03/48] Update materials

---
 mirgecom/materials/carbon_fiber.py | 96 ++++++++++++++++--------------
 mirgecom/materials/tacot.py        |  5 +-
 2 files changed, 52 insertions(+), 49 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 93464f574..7cd503d60 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -28,8 +28,8 @@
 THE SOFTWARE.
 """
 
+from typing import Optional
 from abc import abstractmethod
-import numpy as np
 from meshmode.dof_array import DOFArray
 from mirgecom.wall_model import PorousWallEOS
 from pytools.obj_array import make_obj_array
@@ -38,48 +38,49 @@
 class Oxidation:
     """Abstract interface for wall oxidation model.
 
+    .. automethod:: effective_surface_area_fiber
     .. automethod:: get_source_terms
     """
 
     @abstractmethod
-    def get_source_terms(self, temperature: DOFArray,
-            tau: DOFArray, rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+    def effective_surface_area_fiber(self, progress: DOFArray) -> DOFArray:
+        r"""Evaluate the effective surface of the fibers."""
+        raise NotImplementedError()
+
+    @abstractmethod
+    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
+            rhoY_o2: Optional[DOFArray] = None) -> DOFArray:  # noqa N803
         r"""Source terms of fiber oxidation."""
         raise NotImplementedError()
 
 
-# TODO per MTC review, can we generalize the oxidation model?
-# should we keep this in the driver?
-class Y2_Oxidation_Model(Oxidation):  # noqa N801
+class OxidationModel(Oxidation):
     """Evaluate the source terms for the Y2 model of carbon fiber oxidation.
 
-    .. automethod:: puma_effective_surface_area
+    The user must specify in the driver the functions for oxidation.
+    (Tentatively) Generalizing this class makes it easier for adding new,
+    more complex models and for UQ runs.
+
+    .. automethod:: effective_surface_area_fiber
     .. automethod:: get_source_terms
     """
 
-    def puma_effective_surface_area(self, progress) -> DOFArray:
-        """Polynomial fit based on PUMA data.
-
-        Parameters
-        ----------
-        progress: meshmode.dof_array.DOFArray
-            the rate of decomposition of the fibers
-        """
-        # Original fit function: -1.1012e5*x**2 - 0.0646e5*x + 1.1794e5
-        # Rescale by x==0 value and rearrange
-        return 1.1794e5*(1.0 - 0.0547736137*progress - 0.9336950992*progress**2)
+    def __init__(self, surface_area_func, oxidation_func):
+        self._surface_func = surface_area_func
+        self._oxidation_func = oxidation_func
 
-    def _get_wall_effective_surface_area_fiber(self, progress) -> DOFArray:
-        """Evaluate the effective surface of the fibers.
+    def effective_surface_area_fiber(self, progress: DOFArray) -> DOFArray:
+        r"""Evaluate the effective surface of the fibers.
 
         Parameters
         ----------
         progress: meshmode.dof_array.DOFArray
-            the rate of decomposition of the fibers
+            the rate of decomposition of the fibers, defined as 1.0 - $\tau$.
         """
-        return self.puma_effective_surface_area(progress)
+        return self._surface_func(progress)
 
-    def get_source_terms(self, temperature, tau, rhoY_o2) -> DOFArray:  # noqa N803
+    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
+            rhoY_o2: Optional[DOFArray] = None) -> DOFArray:  # noqa N803
         """Return the effective source terms for the oxidation.
 
         Parameters
@@ -87,23 +88,14 @@ def get_source_terms(self, temperature, tau, rhoY_o2) -> DOFArray:  # noqa N803
         temperature: meshmode.dof_array.DOFArray
         tau: meshmode.dof_array.DOFArray
             the progress ratio of the oxidation
-        ox_mass: meshmode.dof_array.DOFArray
+        rhoY_o2: meshmode.dof_array.DOFArray
             the mass fraction of oxygen
         """
-        actx = temperature.array_context
-
-        mw_o = 15.999
-        mw_o2 = mw_o*2
-        mw_co = 28.010
-        univ_gas_const = 8314.46261815324
-
-        eff_surf_area = self._get_wall_effective_surface_area_fiber(1.0-tau)
-        alpha = (
-            (0.00143+0.01*actx.np.exp(-1450.0/temperature))
-            / (1.0+0.0002*actx.np.exp(13000.0/temperature)))
-        k = alpha*actx.np.sqrt(
-            (univ_gas_const*temperature)/(2.0*np.pi*mw_o2))
-        return (mw_co/mw_o2 + mw_o/mw_o2 - 1)*rhoY_o2*k*eff_surf_area
+        area = self.effective_surface_area_fiber(1.0-tau)
+        if rhoY_o2:
+            return self._oxidation_func(temperature=temperature, fiber_area=area,
+                                        rhoY_o2=rhoY_o2)
+        return self._oxidation_func(temperature=temperature, fiber_area=area)
 
 
 class FiberEOS(PorousWallEOS):
@@ -129,14 +121,27 @@ class FiberEOS(PorousWallEOS):
     """
 
     def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass):
-        """Bulk density considering the porosity and intrinsic density."""
+        """Bulk density considering the porosity and intrinsic density.
+
+        Parameters
+        ----------
+        dim: int
+            geometrical dimension of the problem.
+        virgin_mass: float
+            initial mass of the material.
+        char_mass: float
+            final mass when the decomposition is complete.
+        anisotropic_direction: int
+            For orthotropic materials, this indicates the normal direction
+            where the properties are different than in-plane.
+        """
         self._char_mass = char_mass
         self._virgin_mass = virgin_mass
         self._dim = dim
         self._anisotropic_dir = anisotropic_direction
 
-        if anisotropic_direction > dim:
-            raise ValueError("Anisotropic axis must be less or equal than dim.")
+        if anisotropic_direction >= dim:
+            raise ValueError("Anisotropic axis must be less than dim.")
 
     def void_fraction(self, tau: DOFArray) -> DOFArray:
         r"""Return the volumetric fraction $\epsilon$ filled with gas.
@@ -197,9 +202,10 @@ def volume_fraction(self, tau: DOFArray) -> DOFArray:
     def permeability(self, tau: DOFArray) -> DOFArray:
         r"""Permeability $K$ of the porous material."""
         # FIXME find a relation to make it change as a function of "tau"
+        # TODO: the relation depends on the coupling model. Postpone it for now.
         actx = tau.array_context
-        permeability = np.zeros(self._dim,)
-        permeability[:] = 5.57e-11 + actx.np.zeros_like(tau)
+        permeability = make_obj_array([5.57e-11 + actx.np.zeros_like(tau)
+                                       for _ in range(0, self._dim)])
         permeability[self._anisotropic_dir] = 2.62e-11 + actx.np.zeros_like(tau)
         return permeability
 
@@ -212,9 +218,7 @@ def emissivity(self, temperature=None, tau=None) -> DOFArray:
 
     def tortuosity(self, tau: DOFArray) -> DOFArray:
         r"""Tortuosity $\eta$ affects the species diffusivity."""
-        # FIXME find a relation to make it change as a function of "tau"
-        actx = tau.array_context
-        return 1.1 + actx.np.zeros_like(tau)
+        return 1.1*tau + 1.0*(1.0 - tau)
 
     def decomposition_progress(self, mass: DOFArray) -> DOFArray:
         r"""Evaluate the mass loss progress ratio $\tau$ of the oxidation."""
diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index 8f555a740..f17b95293 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -38,9 +38,8 @@
 from mirgecom.wall_model import PorousWallEOS
 
 
-# TODO generalize? define only abstract class and keep this in the driver?
 class Pyrolysis:
-    r"""Evaluate the source terms for the pyrolysis decomposition.
+    r"""Evaluate the source terms for the pyrolysis decomposition of TACOT.
 
     The source terms follow as Arrhenius-like equation given by
 
@@ -90,8 +89,8 @@ def get_source_terms(self, temperature, chi):
                 0.0, (
                     -(30.*((chi[0] - 0.00)/30.)**3)*12000.
                     * actx.np.exp(-8556.000/temperature))),
-            actx.np.where(actx.np.less(temperature, self._Tcrit[1]),
             # reaction 2
+            actx.np.where(actx.np.less(temperature, self._Tcrit[1]),
                 0.0, (
                     -(90.*((chi[1] - 60.0)/90.)**3)*4.48e9
                     * actx.np.exp(-20444.44/temperature))),

From 52cc1a0e9f215b2b505f4c0fb109473925d5cf32 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:21:23 -0500
Subject: [PATCH 04/48] Update mechanism

---
 mirgecom/mechanisms/uiuc_8sp_phenol.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mirgecom/mechanisms/uiuc_8sp_phenol.yaml b/mirgecom/mechanisms/uiuc_8sp_phenol.yaml
index 4baad89dc..d32aff486 100644
--- a/mirgecom/mechanisms/uiuc_8sp_phenol.yaml
+++ b/mirgecom/mechanisms/uiuc_8sp_phenol.yaml
@@ -152,8 +152,8 @@ species:
     model: NASA7
     temperature-ranges: [300.0, 825.0, 1300.0]
     data:
-    - [18.72619129224749 , -0.18263954836161317 , 0.0008049729732394601 , -1.4792419208580702e-06 , 1.0019289891493458e-09 , -22603.659214472158, 0.683010238]
-    - [-19060.39392965401 , 74.3800489123027 , -0.10706969363457917 , 6.755482446169158e-05 , -1.5785126105941202e-08 , 3822285.2052348647, -3.20502331]
+    - [16.039299967161238 , 0.0 , 0.0 , 0.0 , 0.0 , -26840.35489814738, 0.0]
+    - [16.039299967161238 , 0.0 , 0.0 , 0.0 , 0.0 , -26840.35489814738, 0.0]
   transport:
     model: gas
     geometry: linear

From 246032fdb7ecee242224a6df01f3f91a93d19155 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:21:57 -0500
Subject: [PATCH 05/48] Fix docs, comments and emissivity for wall_model

---
 mirgecom/wall_model.py | 22 +++++++++++++---------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index bb7983ca5..20d5ac4db 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -107,7 +107,7 @@ def array_context(self):
 
 
 @dataclass_array_container
-@dataclass(frozen=True)
+@dataclass(frozen=True, eq=False)
 class SolidWallDependentVars:
     """Wall dependent variables for heat conduction only materials."""
 
@@ -206,7 +206,6 @@ class PorousWallVars:
     .. attribute:: material_densities
     .. attribute:: tau
     .. attribute:: void_fraction
-    .. attribute:: emissivity
     .. attribute:: permeability
     .. attribute:: tortuosity
     .. attribute:: density
@@ -215,7 +214,6 @@ class PorousWallVars:
     material_densities: Union[DOFArray, np.ndarray]
     tau: DOFArray
     void_fraction: DOFArray
-    emissivity: DOFArray
     permeability: DOFArray
     tortuosity: DOFArray
     density: DOFArray
@@ -288,12 +286,13 @@ def decomposition_progress(self, mass: DOFArray) -> DOFArray:
 
 
 # FIXME: Generalize TransportModel interface to accept state variables
-# other than fluid cv
+# other than fluid cv following
+# https://github.com/illinois-ceesd/mirgecom/pull/935#discussion_r1298730910
 class PorousWallTransport:
     r"""Transport model for porous media flow.
 
-    Takes any transport model and modifies it to consider the interaction
-    with the porous materials.
+    Takes any gas-only transport model and modifies it to consider the
+    interaction with the porous materials.
 
     .. automethod:: __init__
     .. automethod:: bulk_viscosity
@@ -383,6 +382,11 @@ class PorousFlowModel:
         Transport class that governs how the gas flows through the porous
         media. This is accounted for in :class:`PorousWallTransport`
 
+    .. attribute:: temperature_iteration
+
+        Number of iterations for temperature evaluation using Newton's method.
+        Defaults to 3 if not specified.
+
     It also include functions that combine the properties of the porous
     material and the gas that permeates, yielding the actual porous flow EOS:
 
@@ -434,7 +438,7 @@ def get_temperature(self, cv: ConservedVars, wv: PorousWallVars,
             T^{n+1} = T^n -
                 \frac
                 {\epsilon_g \rho_g e_g + \rho_s h_s - \rho e}
-                {\epsilon_g \rho_g C_{p_g} + \epsilon_s \rho_s C_{p_s}}
+                {\epsilon_g \rho_g C_{v_g} + \epsilon_s \rho_s C_{p_s}}
 
         """
         if isinstance(tseed, DOFArray) is False:
@@ -469,7 +473,7 @@ def get_pressure(self, cv: ConservedVars, wv: PorousWallVars,
 
     def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
                  temperature: DOFArray) -> DOFArray:
-        r"""Return the enthalpy of the gas+solid material.
+        r"""Return the internal energy of the gas+solid material.
 
         .. math::
             \rho e = \epsilon_s \rho_s e_s + \epsilon_g \rho_g e_g
@@ -483,7 +487,7 @@ def heat_capacity(self, cv: ConservedVars, wv: PorousWallVars,
         r"""Return the heat capacity of the gas+solid material.
 
         .. math::
-            \rho e = \epsilon_s \rho_s {C_p}_s + \epsilon_g \rho_g {C_v}_g
+            \rho C_v = \epsilon_s \rho_s {C_p}_s + \epsilon_g \rho_g {C_v}_g
         """
         return (cv.mass*self.eos.heat_capacity_cv(cv, temperature)
                 + wv.density*self.wall_eos.heat_capacity(temperature, wv.tau))

From 857459f8ae2a385f30016e0470aaaab5c2434d0d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:25:08 -0500
Subject: [PATCH 06/48] docs

---
 mirgecom/materials/carbon_fiber.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 7cd503d60..f6673282c 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -55,7 +55,7 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
 
 
 class OxidationModel(Oxidation):
-    """Evaluate the source terms for the Y2 model of carbon fiber oxidation.
+    """Evaluate the source terms for the carbon fiber oxidation.
 
     The user must specify in the driver the functions for oxidation.
     (Tentatively) Generalizing this class makes it easier for adding new,

From 6f0b880ffe7d775393ddbb14fd2161fe00c8c774 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:55:02 -0500
Subject: [PATCH 07/48] Add test for wall model

---
 test/test_wallmodel.py | 139 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 139 insertions(+)
 create mode 100644 test/test_wallmodel.py

diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
new file mode 100644
index 000000000..1c86c6bef
--- /dev/null
+++ b/test/test_wallmodel.py
@@ -0,0 +1,139 @@
+"""Test wall-model related functions."""
+
+__copyright__ = """Copyright (C) 2023 University of Illinois Board of Trustees"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import pytest
+import cantera
+import numpy as np
+from pytools.obj_array import make_obj_array
+from grudge import op
+from meshmode.dof_array import DOFArray
+from meshmode.array_context import (  # noqa
+    pytest_generate_tests_for_pyopencl_array_context
+    as pytest_generate_tests)
+from mirgecom.simutil import get_box_mesh
+from mirgecom.discretization import create_discretization_collection
+from mirgecom.eos import PyrometheusMixture
+from mirgecom.transport import SimpleTransport
+from mirgecom.gas_model import make_fluid_state
+from mirgecom.wall_model import PorousFlowModel, PorousWallTransport
+from mirgecom.materials.initializer import PorousWallInitializer
+from mirgecom.mechanisms import get_mechanism_input
+from mirgecom.thermochemistry import get_pyrometheus_wrapper_class_from_cantera
+
+
+@pytest.mark.parametrize("order", [1, 4])
+@pytest.mark.parametrize("my_material", ["fiber", "composite"])
+def test_wall_eos(actx_factory, order, my_material):
+    """Check the wall degradation model."""
+    actx = actx_factory()
+
+    dim = 2
+
+    tol = 1e-8
+
+    nelems = 20
+    mesh = get_box_mesh(dim, -0.1, 0.1, nelems)
+    dcoll = create_discretization_collection(
+        actx, mesh, order=order, quadrature_order=2*order+1)
+
+    nodes = actx.thaw(dcoll.nodes())
+    zeros = actx.np.zeros_like(nodes[0])
+
+    # {{{ EOS initialization
+
+    mech_input = get_mechanism_input("uiuc_8sp_phenol")
+    cantera_soln = cantera.Solution(name="gas", yaml=mech_input)
+    pyro_obj = get_pyrometheus_wrapper_class_from_cantera(
+        cantera_soln, temperature_niter=3)(actx.np)
+
+    nspecies = pyro_obj.num_species
+
+    x = make_obj_array([0.0 for i in range(nspecies)])
+    x[cantera_soln.species_index("O2")] = 0.21
+    x[cantera_soln.species_index("N2")] = 0.79
+
+    cantera_soln.TPX = 900.0, 101325.0, x
+    _, _, y = cantera_soln.TDY
+
+    eos = PyrometheusMixture(pyro_obj, temperature_guess=cantera_soln.T)
+
+    # }}}
+
+    # {{{ Initialize wall model
+
+    if my_material == "fiber":
+        import mirgecom.materials.carbon_fiber as material_sample
+        material = material_sample.FiberEOS(
+            dim=dim, anisotropic_direction=0, char_mass=0., virgin_mass=168.)
+        material_densities = 0.12*1400.0
+
+    if my_material == "composite":
+        import mirgecom.materials.tacot as material_sample
+        material = material_sample.TacotEOS(char_mass=220., virgin_mass=280.)
+        material_densities = np.empty((3,), dtype=object)
+        material_densities[0] = 30.0
+        material_densities[1] = 90.0
+        material_densities[2] = 160.
+
+    # }}}
+
+    # {{{ Gas model
+
+    base_transport = SimpleTransport()
+    solid_transport = PorousWallTransport(base_transport=base_transport)
+    gas_model = PorousFlowModel(eos=eos, transport=solid_transport,
+                                wall_eos=material)
+
+    # }}}
+
+    pressure = 101325.0 + zeros
+    temperature = 900.0 + zeros
+
+    sample_init = PorousWallInitializer(
+        pressure=pressure, temperature=temperature, species=y,
+        material_densities=material_densities)
+
+    cv, solid_densities = sample_init(nodes, gas_model)
+
+    solid_state = make_fluid_state(cv=cv, gas_model=gas_model,
+        material_densities=solid_densities, temperature_seed=900.0)
+
+    wv = solid_state.wv
+
+    assert actx.to_numpy(op.norm(dcoll, wv.tau - 1.0, np.inf)) < tol
+
+    assert isinstance(solid_state.wv.density, DOFArray)
+
+    if my_material == "fiber":
+        assert np.max(actx.to_numpy(wv.density - 168.0)) < tol
+        assert np.max(actx.to_numpy(wv.void_fraction)) - 1.00 < tol
+
+    if my_material == "composite":
+        assert np.max(actx.to_numpy(wv.density - 280.0)) < tol
+        assert np.max(actx.to_numpy(wv.void_fraction)) - 1.00 < tol
+
+    assert actx.to_numpy(
+        op.norm(dcoll, solid_state.pressure - 101325.0, np.inf)) < tol
+    assert actx.to_numpy(
+        op.norm(dcoll, solid_state.temperature - 900.0, np.inf)) < tol

From b049a4b0a4256191f07e4c0f8aeeace3b5f5b69b Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:55:17 -0500
Subject: [PATCH 08/48] Update wall_model

---
 mirgecom/wall_model.py | 71 ++++++++++++++++++++++++++----------------
 1 file changed, 45 insertions(+), 26 deletions(-)

diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 20d5ac4db..464bb29d0 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -111,6 +111,8 @@ def array_context(self):
 class SolidWallDependentVars:
     """Wall dependent variables for heat conduction only materials."""
 
+    tau: DOFArray
+    density: Union[DOFArray, np.ndarray]
     thermal_conductivity: DOFArray
     temperature: DOFArray
 
@@ -132,7 +134,7 @@ class SolidWallModel:
     driver-defined functions that prescribe the actual properties and their
     spatial distribution.
 
-    .. automethod:: density
+    .. automethod:: solid_density
     .. automethod:: heat_capacity
     .. automethod:: enthalpy
     .. automethod:: thermal_diffusivity
@@ -141,55 +143,72 @@ class SolidWallModel:
     .. automethod:: dependent_vars
     """
 
-    def __init__(self, density_func, enthalpy_func, heat_capacity_func,
-                 thermal_conductivity_func):
-        self._density_func = density_func
+    def __init__(self, enthalpy_func, heat_capacity_func,
+                 thermal_conductivity_func, decomposition_func=None):
         self._enthalpy_func = enthalpy_func
         self._heat_capacity_func = heat_capacity_func
         self._thermal_conductivity_func = thermal_conductivity_func
+        self._decomposition_func = decomposition_func
 
-    def density(self) -> DOFArray:
-        """Return the wall density for all components."""
-        return self._density_func()
+    def solid_density(self, material_densities) -> DOFArray:
+        r"""Return the solid density $\rho_s$."""
+        if isinstance(material_densities, DOFArray):
+            return material_densities
+        return sum(material_densities)
+
+    def decomposition_progress(self, material_densities) -> DOFArray:
+        r"""Evaluate the progress ratio $\tau$ of the decomposition."""
+        mass = self.solid_density(material_densities)
+        if self._decomposition_func:
+            # for materials that decompose when heated
+            return self._decomposition_func(mass)
+        actx = mass.array_context
+        return actx.np.zeros_like(mass) + 1.0
 
-    def heat_capacity(self, temperature: DOFArray = None) -> DOFArray:
-        """Return the wall heat_capacity for all components."""
-        return self._heat_capacity_func(temperature)
+    def heat_capacity(self, temperature: DOFArray = None, tau=None) -> DOFArray:
+        """Return the wall heat_capacity."""
+        return self._heat_capacity_func(temperature=temperature, tau=tau)
 
-    def enthalpy(self, temperature: DOFArray) -> DOFArray:
-        """Return the wall enthalpy for all components."""
-        return self._enthalpy_func(temperature)
+    def enthalpy(self, temperature: DOFArray, tau=None) -> DOFArray:
+        """Return the wall enthalpy."""
+        return self._enthalpy_func(temperature=temperature, tau=tau)
 
-    def thermal_diffusivity(self, mass: DOFArray, temperature: DOFArray,
-            thermal_conductivity: DOFArray = None) -> DOFArray:
+    def thermal_diffusivity(self, solid_state: SolidWallState) -> DOFArray:
         """Return the wall thermal diffusivity for all components."""
-        if thermal_conductivity is None:
-            thermal_conductivity = self.thermal_conductivity(temperature)
-        return thermal_conductivity/(mass * self.heat_capacity(temperature))
+        tau = self.decomposition_progress(solid_state.cv.mass)
+        mass = self.solid_density(solid_state.cv.mass)
+        kappa = self.thermal_conductivity(solid_state.dv.temperature, tau)
+        cp = self.heat_capacity(solid_state.dv.temperature, tau)
+        return kappa/(mass * cp)
 
-    def thermal_conductivity(self, temperature: DOFArray) -> DOFArray:
-        """Return the wall thermal conductivity for all components."""
-        return self._thermal_conductivity_func(temperature)
+    def thermal_conductivity(self, temperature: DOFArray, tau=None) -> DOFArray:
+        """Return the wall thermal conductivity."""
+        return self._thermal_conductivity_func(temperature=temperature, tau=tau)
 
     def get_temperature(self, wv: SolidWallConservedVars,
             tseed: DOFArray = None, niter: int = 3) -> DOFArray:
         """Evaluate the temperature based on the energy."""
         if tseed is not None:
+            tau = self.decomposition_progress(wv.mass)
             temp = tseed*1.0
             for _ in range(0, niter):
-                h = self.enthalpy(temp)
-                cp = self.heat_capacity(temp)
-                temp = temp - (h - wv.energy/wv.mass)/cp
+                h = self.enthalpy(temperature=temp, tau=tau)
+                cp = self.heat_capacity(temperature=temp, tau=tau)
+                temp = temp - (h - wv.energy/self.solid_density(wv.mass))/cp
             return temp
 
-        return wv.energy/(self.density()*self.heat_capacity())
+        return wv.energy/(self.solid_density(wv.mass)*self.heat_capacity())
 
     def dependent_vars(self, wv: SolidWallConservedVars,
             tseed: DOFArray = None, niter: int = 3) -> SolidWallDependentVars:
         """Return solid wall dependent variables."""
+        tau = self.decomposition_progress(wv.mass)
+        density = self.solid_density(wv.mass)
         temperature = self.get_temperature(wv, tseed, niter)
-        kappa = self.thermal_conductivity(temperature)
+        kappa = self.thermal_conductivity(temperature=temperature, tau=tau)
         return SolidWallDependentVars(
+            tau=tau,
+            density=density,
             thermal_conductivity=kappa,
             temperature=temperature)
 

From 28a6d27f60be1368ca0404c2e104dc6269f519b1 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 13:55:32 -0500
Subject: [PATCH 09/48] Update materials initializer

---
 mirgecom/materials/initializer.py | 45 ++++++++++++++++++++-----------
 1 file changed, 30 insertions(+), 15 deletions(-)

diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index b0977bdce..7419bc8c8 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -36,8 +36,9 @@
 class SolidWallInitializer:
     """Initializer for heat conduction only materials."""
 
-    def __init__(self, temperature):
+    def __init__(self, temperature, material_densities):
         self._temp = temperature
+        self._mass = material_densities
 
     def __call__(self, x_vec, wall_model):
         """Evaluate the wall+gas properties for porous materials.
@@ -50,8 +51,14 @@ def __call__(self, x_vec, wall_model):
             Equation of state class
         """
         actx = x_vec[0].array_context
-        mass = wall_model.density() + actx.np.zeros_like(x_vec[0])
-        energy = mass * wall_model.enthalpy(self._temp)
+
+        mass = self._mass + actx.np.zeros_like(x_vec[0])
+        solid_mass = wall_model.solid_density(mass)
+        tau = wall_model.decomposition_progress(mass)
+
+        temperature = self._temp + actx.np.zeros_like(x_vec[0])
+        energy = solid_mass * wall_model.enthalpy(temperature=temperature,
+                                                  tau=tau)
         return SolidWallConservedVars(mass=mass, energy=energy)
 
 
@@ -67,7 +74,7 @@ def __init__(self, temperature, species, material_densities,
         self._temp = temperature
         self._wall_density = material_densities
 
-    def __call__(self, dim, x_vec, gas_model):
+    def __call__(self, x_vec, gas_model):
         """Evaluate the wall+gas properties for porous materials.
 
         Parameters
@@ -76,27 +83,33 @@ def __call__(self, dim, x_vec, gas_model):
             Nodal coordinates
         gas_model: :class:`mirgecom.wall_model.PorousFlowModel`
             Equation of state class
+
+        Returns
+        -------
+        cv: :class:`mirgecom.fluid.ConservedVars`
+            The conserved variables for porous flows.
+        wall_density: numpy.ndarray or :class:`meshmode.dof_array.DOFArray`
+            The densities of each one of the materials
         """
         actx = x_vec[0].array_context
         zeros = actx.np.zeros_like(x_vec[0])
+        dim = x_vec.shape[0]
 
         temperature = self._temp + zeros
+
         species_mass_frac = self._y + zeros
-        wall_density = self._wall_density + zeros
 
+        wall_density = self._wall_density + zeros
         tau = gas_model.decomposition_progress(wall_density)
+        eps_rho_solid = gas_model.solid_density(wall_density)
 
         eps_gas = gas_model.wall_eos.void_fraction(tau)
         if self._mass is None:
             pressure = self._pres + zeros
-            eps_rho_gas = eps_gas*gas_model.eos.get_density(pressure,
-                temperature, species_mass_frac)
+            eps_rho_gas = eps_gas*gas_model.eos.get_density(
+                pressure, temperature, species_mass_frac)
         else:
-            density = self._mass + zeros
-            eps_rho_gas = eps_gas*density
-
-        # internal energy (kinetic energy is neglected)
-        eps_rho_solid = gas_model.solid_density(wall_density)
+            eps_rho_gas = eps_gas*self._mass
 
         bulk_energy = (
             eps_rho_solid*gas_model.wall_eos.enthalpy(temperature, tau)
@@ -104,9 +117,11 @@ def __call__(self, dim, x_vec, gas_model):
                                                             species_mass_frac)
         )
 
-        momentum = make_obj_array([zeros, zeros])
+        momentum = make_obj_array([zeros for _ in range(dim)])
 
         species_mass = eps_rho_gas*species_mass_frac
 
-        return make_conserved(dim=dim, mass=eps_rho_gas, energy=bulk_energy,
-            momentum=momentum, species_mass=species_mass)
+        cv = make_conserved(dim=dim, mass=eps_rho_gas, energy=bulk_energy,
+                            momentum=momentum, species_mass=species_mass)
+
+        return cv, wall_density

From 8ec893a5e52fc52202b538de23abaa25a1c7bc9d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 14:27:03 -0500
Subject: [PATCH 10/48] Update inviscid flux argument

---
 mirgecom/inviscid.py     | 35 +++++++++++++++++++++++++----------
 mirgecom/navierstokes.py |  3 ++-
 2 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py
index f1c713a6d..da045614f 100644
--- a/mirgecom/inviscid.py
+++ b/mirgecom/inviscid.py
@@ -44,6 +44,8 @@
 """
 
 import numpy as np
+from arraycontext import outer
+from meshmode.dof_array import DOFArray
 from meshmode.discretization.connection import FACE_RESTR_ALL
 from grudge.dof_desc import (
     DD_VOLUME_ALL,
@@ -56,12 +58,10 @@
     ConservedVars
 )
 from mirgecom.utils import normalize_boundaries
-
-from arraycontext import outer
-from meshmode.dof_array import DOFArray
+from mirgecom.gas_model import PorousFlowFluidState
 
 
-def inviscid_flux(state):
+def inviscid_flux(state, gas_model=None):
     r"""Compute the inviscid flux vectors from fluid conserved vars *cv*.
 
     The inviscid fluxes are
@@ -89,11 +89,25 @@ def inviscid_flux(state):
         conservation equation.
     """
     mass_flux = state.momentum_density
-    energy_flux = state.velocity * (state.energy_density + state.pressure)
+
+    # only the gas part goes into the energy flux
+    if isinstance(state, PorousFlowFluidState):
+        gas_energy = state.mass_density*(
+            gas_model.eos.get_internal_energy(
+                temperature=state.temperature,
+                species_mass_fractions=state.cv.species_mass_fractions)
+            + 0.5*np.dot(state.velocity, state.velocity)
+        )
+        energy_flux = state.velocity * (gas_energy
+                                        + state.wv.void_fraction*state.pressure)
+    else:
+        energy_flux = state.velocity * (state.energy_density + state.pressure)
+
     mom_flux = (
         state.mass_density * np.outer(state.velocity, state.velocity)
         + np.eye(state.dim)*state.pressure
     )
+
     species_mass_flux = \
         state.velocity*state.species_mass_density.reshape(-1, 1)
 
@@ -145,16 +159,17 @@ def inviscid_facial_flux_rusanov(state_pair, gas_model, normal):
     actx = state_pair.int.array_context
     lam = actx.np.maximum(state_pair.int.wavespeed, state_pair.ext.wavespeed)
     from mirgecom.flux import num_flux_lfr
-    return num_flux_lfr(f_minus_normal=inviscid_flux(state_pair.int)@normal,
-                        f_plus_normal=inviscid_flux(state_pair.ext)@normal,
-                        q_minus=state_pair.int.cv,
-                        q_plus=state_pair.ext.cv, lam=lam)
+    return num_flux_lfr(
+        f_minus_normal=inviscid_flux(state_pair.int, gas_model=gas_model)@normal,
+        f_plus_normal=inviscid_flux(state_pair.ext, gas_model=gas_model)@normal,
+        q_minus=state_pair.int.cv,
+        q_plus=state_pair.ext.cv, lam=lam)
 
 
 def inviscid_facial_flux_hll(state_pair, gas_model, normal):
     r"""High-level interface for inviscid facial flux using HLL numerical flux.
 
-    The Harten, Lax, van Leer approximate riemann numerical flux is calculated as:
+    The Harten, Lax, van Leer approximate Riemann numerical flux is calculated as:
 
     .. math::
 
diff --git a/mirgecom/navierstokes.py b/mirgecom/navierstokes.py
index 3d7f26ec3..067b63d10 100644
--- a/mirgecom/navierstokes.py
+++ b/mirgecom/navierstokes.py
@@ -549,7 +549,8 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0,
     # To get separate inviscid operator, set inviscid_fluid_operator explicitly
     # or use ESDG.
     if inviscid_terms_on and inviscid_fluid_operator is None:
-        vol_term = vol_term - inviscid_flux(state=vol_state_quad)
+        vol_term = vol_term - inviscid_flux(state=vol_state_quad,
+                                            gas_model=gas_model)
         bnd_term = bnd_term - inviscid_flux_on_element_boundary(
             dcoll, gas_model, boundaries, inter_elem_bnd_states_quad,
             domain_bnd_states_quad, quadrature_tag=quadrature_tag,

From 441b352c6d9e0ca48d9d68a0cb1a63a64740b7eb Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 14:27:39 -0500
Subject: [PATCH 11/48] Update material and properties when coupling with fluid

---
 mirgecom/materials/carbon_fiber.py | 17 ++++++++++++-----
 mirgecom/materials/tacot.py        |  7 ++++++-
 2 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index f6673282c..2dc1c1e4f 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -109,6 +109,7 @@ class FiberEOS(PorousWallEOS):
         \tau = \frac{m}{m_0} = \frac{\rho_i \epsilon}{\rho_i \epsilon_0}
              = \frac{r^2}{r_0^2}
 
+    .. automethod:: __init__
     .. automethod:: void_fraction
     .. automethod:: enthalpy
     .. automethod:: heat_capacity
@@ -120,7 +121,8 @@ class FiberEOS(PorousWallEOS):
     .. automethod:: decomposition_progress
     """
 
-    def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass):
+    def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass,
+                 timescale=1.0):
         """Bulk density considering the porosity and intrinsic density.
 
         Parameters
@@ -129,16 +131,21 @@ def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass):
             geometrical dimension of the problem.
         virgin_mass: float
             initial mass of the material.
+
         char_mass: float
             final mass when the decomposition is complete.
         anisotropic_direction: int
             For orthotropic materials, this indicates the normal direction
             where the properties are different than in-plane.
+        timescale: float
+            Modifies the thermal conductivity and the radiation emission to
+            increase/decrease the wall time-scale. Defaults to 1.0 (no changes).
         """
         self._char_mass = char_mass
         self._virgin_mass = virgin_mass
         self._dim = dim
         self._anisotropic_dir = anisotropic_direction
+        self._timescale = timescale
 
         if anisotropic_direction >= dim:
             raise ValueError("Anisotropic axis must be less than dim.")
@@ -186,13 +193,12 @@ def thermal_conductivity(self, temperature, tau=None) -> DOFArray:
             + 1.93072961e-10*temperature**3 - 3.52595953e-07*temperature**2
             + 4.54935976e-04*temperature**1 + 5.08960039e-02)
 
-        # initialize with the in-plane value
+        # initialize with the in-plane value then modify the normal direction
         kappa = make_obj_array([kappa_ij for _ in range(self._dim)])
-        # modify only the normal direction
         kappa[self._anisotropic_dir] = kappa_k
 
         # account for fiber shrinkage via "tau"
-        return kappa*tau
+        return kappa*tau*self._timescale
 
     # ~~~~~~~~ other properties
     def volume_fraction(self, tau: DOFArray) -> DOFArray:
@@ -207,11 +213,12 @@ def permeability(self, tau: DOFArray) -> DOFArray:
         permeability = make_obj_array([5.57e-11 + actx.np.zeros_like(tau)
                                        for _ in range(0, self._dim)])
         permeability[self._anisotropic_dir] = 2.62e-11 + actx.np.zeros_like(tau)
+
         return permeability
 
     def emissivity(self, temperature=None, tau=None) -> DOFArray:
         """Emissivity for energy radiation."""
-        return (
+        return self._timescale * (
             + 2.26413679e-18*temperature**5 - 2.03008004e-14*temperature**4
             + 7.05300324e-11*temperature**3 - 1.22131715e-07*temperature**2
             + 1.21137817e-04*temperature**1 + 8.66656964e-01)
diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index f17b95293..058870fd2 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -215,7 +215,12 @@ def decomposition_progress(self, mass: DOFArray) -> DOFArray:
         .. math::
             \tau = \frac{\rho_0}{\rho_0 - \rho_c}
                     \left( 1 - \frac{\rho_c}{\rho(t)} \right)
+
+        Note that $\rho(t)$ is the mass of resin only, not including fibers,
+        permeating gas or gas-only regions.
         """
+        actx = mass.array_context
         char_mass = self._char_mass
         virgin_mass = self._virgin_mass
-        return virgin_mass/(virgin_mass - char_mass)*(1.0 - char_mass/mass)
+        resin_mass = actx.np.maximum(char_mass, mass)
+        return virgin_mass/(virgin_mass - char_mass)*(1.0 - char_mass/resin_mass)

From 18177f23caab8cec2339daed39416e5ed9cbf565 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 14:28:32 -0500
Subject: [PATCH 12/48] Update wall model to handle IdealSingleGas

---
 mirgecom/wall_model.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 464bb29d0..757acf796 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -393,8 +393,7 @@ class PorousFlowModel:
 
     .. attribute:: eos
 
-        The thermophysical properties of the gas and its EOS. For now, only
-        mixtures are considered.
+        The thermophysical properties of the gas and its EOS.
 
     .. attribute:: transport
 
@@ -417,7 +416,7 @@ class PorousFlowModel:
     .. automethod:: heat_capacity
     """
 
-    eos: MixtureEOS
+    eos: Union[GasEOS, MixtureEOS]
     wall_eos: PorousWallEOS
     transport: PorousWallTransport
     temperature_iteration: int = 3
@@ -488,7 +487,9 @@ def get_pressure(self, cv: ConservedVars, wv: PorousWallVars,
         where $\epsilon \rho$ is stored in :class:`~mirgecom.fluid.ConservedVars`
         and $M$ is the molar mass of the mixture.
         """
-        return self.eos.pressure(cv, temperature)/wv.void_fraction
+        if isinstance(self.eos, MixtureEOS):
+            return self.eos.pressure(cv, temperature)/wv.void_fraction
+        return cv.mass*self.eos.gas_const()*temperature/wv.void_fraction
 
     def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
                  temperature: DOFArray) -> DOFArray:

From f53ffb1f6f02256dea3b37a733e29cc7dd6f5667 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 14:29:03 -0500
Subject: [PATCH 13/48] Add test to test_wallmodel

---
 test/test_wallmodel.py | 43 +++++++++++++++++++++++++++---------------
 1 file changed, 28 insertions(+), 15 deletions(-)

diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index 1c86c6bef..2fba383f7 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -33,7 +33,7 @@
     as pytest_generate_tests)
 from mirgecom.simutil import get_box_mesh
 from mirgecom.discretization import create_discretization_collection
-from mirgecom.eos import PyrometheusMixture
+from mirgecom.eos import IdealSingleGas, PyrometheusMixture
 from mirgecom.transport import SimpleTransport
 from mirgecom.gas_model import make_fluid_state
 from mirgecom.wall_model import PorousFlowModel, PorousWallTransport
@@ -44,7 +44,9 @@
 
 @pytest.mark.parametrize("order", [1, 4])
 @pytest.mark.parametrize("my_material", ["fiber", "composite"])
-def test_wall_eos(actx_factory, order, my_material):
+@pytest.mark.parametrize("my_eos", ["ideal", "mixture"])
+@pytest.mark.parametrize("use_diffused_interface", [True, False])
+def test_wall_eos(actx_factory, order, my_material, my_eos, use_diffused_interface):
     """Check the wall degradation model."""
     actx = actx_factory()
 
@@ -62,21 +64,26 @@ def test_wall_eos(actx_factory, order, my_material):
 
     # {{{ EOS initialization
 
-    mech_input = get_mechanism_input("uiuc_8sp_phenol")
-    cantera_soln = cantera.Solution(name="gas", yaml=mech_input)
-    pyro_obj = get_pyrometheus_wrapper_class_from_cantera(
-        cantera_soln, temperature_niter=3)(actx.np)
+    if my_eos == "mixture":
+        mech_input = get_mechanism_input("uiuc_8sp_phenol")
+        cantera_soln = cantera.Solution(name="gas", yaml=mech_input)
+        pyro_obj = get_pyrometheus_wrapper_class_from_cantera(
+            cantera_soln, temperature_niter=3)(actx.np)
 
-    nspecies = pyro_obj.num_species
+        nspecies = pyro_obj.num_species
 
-    x = make_obj_array([0.0 for i in range(nspecies)])
-    x[cantera_soln.species_index("O2")] = 0.21
-    x[cantera_soln.species_index("N2")] = 0.79
+        x = make_obj_array([0.0 for i in range(nspecies)])
+        x[cantera_soln.species_index("O2")] = 0.21
+        x[cantera_soln.species_index("N2")] = 0.79
 
-    cantera_soln.TPX = 900.0, 101325.0, x
-    _, _, y = cantera_soln.TDY
+        cantera_soln.TPX = 900.0, 101325.0, x
+        _, _, y = cantera_soln.TDY
 
-    eos = PyrometheusMixture(pyro_obj, temperature_guess=cantera_soln.T)
+        eos = PyrometheusMixture(pyro_obj, temperature_guess=cantera_soln.T)
+
+    if my_eos == "ideal":
+        eos = IdealSingleGas()
+        y = None
 
     # }}}
 
@@ -110,9 +117,13 @@ def test_wall_eos(actx_factory, order, my_material):
     pressure = 101325.0 + zeros
     temperature = 900.0 + zeros
 
+    if use_diffused_interface:
+        porous_region = 0.5*(1.0 - actx.np.tanh(nodes[0]/0.005))
+    else:
+        porous_region = None
     sample_init = PorousWallInitializer(
         pressure=pressure, temperature=temperature, species=y,
-        material_densities=material_densities)
+        material_densities=material_densities, porous_region=porous_region)
 
     cv, solid_densities = sample_init(nodes, gas_model)
 
@@ -121,7 +132,9 @@ def test_wall_eos(actx_factory, order, my_material):
 
     wv = solid_state.wv
 
-    assert actx.to_numpy(op.norm(dcoll, wv.tau - 1.0, np.inf)) < tol
+    assert isinstance(wv.tau, DOFArray)
+    if use_diffused_interface is False:
+        assert actx.to_numpy(op.norm(dcoll, wv.tau - 1.0, np.inf)) < tol
 
     assert isinstance(solid_state.wv.density, DOFArray)
 

From 723689323a1aef7254782dd835fd42c83e4e229d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 14:29:43 -0500
Subject: [PATCH 14/48] Add material property to state creation in boundary
 conditions

---
 mirgecom/boundary.py | 62 +++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 55 insertions(+), 7 deletions(-)

diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py
index bc467219d..4206bbe07 100644
--- a/mirgecom/boundary.py
+++ b/mirgecom/boundary.py
@@ -64,6 +64,7 @@
 from mirgecom.utils import project_from_base
 from mirgecom.viscous import viscous_facial_flux_central, viscous_flux
 from mirgecom.inviscid import inviscid_facial_flux_rusanov, inviscid_flux
+from mirgecom.gas_model import PorousFlowFluidState
 
 
 def _ldg_bnd_flux_for_grad(internal_quantity, external_quantity):
@@ -1089,6 +1090,8 @@ def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         # set the normal momentum to 0
         mom_bc = self._slip.momentum_bc(state_minus.momentum_density, nhat)
 
+        # this doesn't recompute the internal energy but only substracts the
+        # kinetic energy from the total energy
         energy_bc = (
             gas_model.eos.internal_energy(state_minus.cv)
             + 0.5*np.dot(mom_bc, mom_bc)/state_minus.mass_density)
@@ -1202,11 +1205,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
             momentum=free_stream_density*free_stream_velocity,
             species_mass=free_stream_spec_mass)
 
+        if isinstance(state_minus, PorousFlowFluidState):
+            material_densities = state_minus.wv.material_densities
+        else:
+            material_densities = None
+
         return make_fluid_state(cv=cv_infinity, gas_model=gas_model,
                                 temperature_seed=free_stream_temperature,
                                 smoothness_mu=state_minus.smoothness_mu,
                                 smoothness_kappa=state_minus.smoothness_kappa,
-                                smoothness_beta=state_minus.smoothness_beta)
+                                smoothness_beta=state_minus.smoothness_beta,
+                                material_densities=material_densities)
 
     def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Return BC fluid state."""
@@ -1337,11 +1346,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
                                     energy=total_energy,
                                     species_mass=state_minus.cv.species_mass)
 
+        if isinstance(state_minus, PorousFlowFluidState):
+            material_densities = state_minus.wv.material_densities
+        else:
+            material_densities = None
+
         return make_fluid_state(cv=cv_outflow, gas_model=gas_model,
                                 temperature_seed=state_minus.temperature,
                                 smoothness_mu=state_minus.smoothness_mu,
                                 smoothness_kappa=state_minus.smoothness_kappa,
-                                smoothness_beta=state_minus.smoothness_beta)
+                                smoothness_beta=state_minus.smoothness_beta,
+                                material_densities=material_densities)
 
     def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Return state."""
@@ -1380,11 +1395,17 @@ def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
             momentum=state_minus.momentum_density,
             species_mass=state_minus.species_mass_density)
 
+        if isinstance(state_minus, PorousFlowFluidState):
+            material_densities = state_minus.wv.material_densities
+        else:
+            material_densities = None
+
         return make_fluid_state(cv=cv_plus, gas_model=gas_model,
                                 temperature_seed=state_minus.temperature,
                                 smoothness_mu=state_minus.smoothness_mu,
                                 smoothness_kappa=state_minus.smoothness_kappa,
-                                smoothness_beta=state_minus.smoothness_beta)
+                                smoothness_beta=state_minus.smoothness_beta,
+                                material_densities=material_densities)
 
     def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs):
         """Get temperature value used in grad(T)."""
@@ -1483,11 +1504,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
                                      momentum=rho_boundary * velocity_boundary,
                                      species_mass=species_mass_boundary)
 
+        if isinstance(state_minus, PorousFlowFluidState):
+            material_densities = state_minus.wv.material_densities
+        else:
+            material_densities = None
+
         return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
                                 temperature_seed=state_minus.temperature,
                                 smoothness_mu=state_minus.smoothness_mu,
                                 smoothness_kappa=state_minus.smoothness_kappa,
-                                smoothness_beta=state_minus.smoothness_beta)
+                                smoothness_beta=state_minus.smoothness_beta,
+                                material_densities=material_densities)
 
     def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Return BC fluid state."""
@@ -1600,11 +1627,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
                                      momentum=rho_boundary*velocity_boundary,
                                      species_mass=species_mass_boundary)
 
+        if isinstance(state_minus, PorousFlowFluidState):
+            material_densities = state_minus.wv.material_densities
+        else:
+            material_densities = None
+
         return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
                                 temperature_seed=state_minus.temperature,
                                 smoothness_mu=state_minus.smoothness_mu,
                                 smoothness_kappa=state_minus.smoothness_kappa,
-                                smoothness_beta=state_minus.smoothness_beta)
+                                smoothness_beta=state_minus.smoothness_beta,
+                                material_densities=material_densities)
 
     def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Return BC fluid state."""
@@ -1965,11 +1998,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
                                      momentum=mass*velocity,
                                      species_mass=mass*species_mass_fracs)
 
+        if isinstance(state_minus, PorousFlowFluidState):
+            material_densities = state_minus.wv.material_densities
+        else:
+            material_densities = None
+
         return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
                                 temperature_seed=state_minus.temperature,
                                 smoothness_mu=state_minus.smoothness_mu,
                                 smoothness_kappa=state_minus.smoothness_kappa,
-                                smoothness_beta=state_minus.smoothness_beta)
+                                smoothness_beta=state_minus.smoothness_beta,
+                                material_densities=material_densities)
 
     def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Return BC fluid state."""
@@ -2067,8 +2106,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
                                      momentum=mass*velocity,
                                      species_mass=mass*self._spec_mass_fracs)
 
+        if isinstance(state_minus, PorousFlowFluidState):
+            material_densities = state_minus.wv.material_densities
+        else:
+            material_densities = None
+
         return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
-                                temperature_seed=state_minus.temperature)
+                                temperature_seed=state_minus.temperature,
+                                smoothness_mu=state_minus.smoothness_mu,
+                                smoothness_kappa=state_minus.smoothness_kappa,
+                                smoothness_beta=state_minus.smoothness_beta,
+                                material_densities=material_densities)
 
     def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Return BC fluid state."""

From aebbe359da025136a8f3be37344fd459024fbac4 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 14:30:17 -0500
Subject: [PATCH 15/48] Update gas model to handle IdealSingleGas in porous
 fluid

---
 mirgecom/gas_model.py | 37 ++++++++++++++++++++++++-------------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py
index da8157ee5..fd29fa85c 100644
--- a/mirgecom/gas_model.py
+++ b/mirgecom/gas_model.py
@@ -417,7 +417,6 @@ def make_fluid_state(cv, gas_model,
         # FIXME per previous review, think of a way to de-couple wall and fluid.
         # ~~~ we need to squeeze wall_eos in gas_model because this is easily
         # accessible everywhere in the code
-
         tau = gas_model.decomposition_progress(material_densities)
         wv = PorousWallVars(
             material_densities=material_densities,
@@ -437,16 +436,28 @@ def make_fluid_state(cv, gas_model,
             cv = limiter_func(cv=cv, wv=wv, pressure=pressure,
                               temperature=temperature, dd=limiter_dd)
 
-        dv = MixtureDependentVars(
-            temperature=temperature,
-            pressure=pressure,
-            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),
-        )
+        from mirgecom.eos import MixtureEOS
+        if isinstance(gas_model.eos, MixtureEOS):
+            dv = MixtureDependentVars(
+                temperature=temperature,
+                pressure=pressure,
+                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)
+            )
+        else:
+            dv = GasDependentVars(
+                temperature=temperature,
+                pressure=pressure,
+                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
+            )
 
         tv = gas_model.transport.transport_vars(cv=cv, dv=dv, wv=wv,
             eos=gas_model.eos, wall_eos=gas_model.wall_eos)
@@ -504,7 +515,7 @@ def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None,
     cv_sd = op.project(dcoll, src, tgt, state.cv)
 
     temperature_seed = None
-    if state.is_mixture:
+    if state.is_mixture or isinstance(gas_model, PorousFlowModel):
         temperature_seed = op.project(dcoll, src, tgt, state.dv.temperature)
 
     if entropy_stable:
@@ -760,7 +771,7 @@ def make_operator_fluid_states(
     ]
 
     tseed_interior_pairs = None
-    if volume_state.is_mixture:
+    if volume_state.is_mixture or isinstance(gas_model, PorousFlowModel):
         # If this is a mixture, we need to exchange the temperature field because
         # mixture pressure (used in the inviscid flux calculations) depends on
         # temperature and we need to seed the temperature calculation for the

From 6db936f61917d469960f2896bd09993d0ffcbe96 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 11 Mar 2024 14:30:55 -0500
Subject: [PATCH 16/48] Update materials initializer to handle IdealSingleGas
 in porous fluid

---
 mirgecom/materials/initializer.py | 33 +++++++++++++++++++++++++++----
 1 file changed, 29 insertions(+), 4 deletions(-)

diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 7419bc8c8..88eca819d 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -63,16 +63,30 @@ def __call__(self, x_vec, wall_model):
 
 
 class PorousWallInitializer:
-    """Initializer for porous materials."""
+    """State initializer for porous materials."""
 
-    def __init__(self, temperature, species, material_densities,
-                 pressure=None, density=None):
+    def __init__(self, temperature, material_densities, species=None,
+                 pressure=None, density=None, porous_region=None):
+        """Initialize the object for porous materials.
 
+        Parameters
+        ----------
+        porous_region:
+            Field describing the homogeneous fluid (0) and porous material (1)
+            portions of the flow. Only used for unified-domain solver without
+            explicit coupling.
+        """
         self._pres = pressure
         self._mass = density
-        self._y = species
         self._temp = temperature
         self._wall_density = material_densities
+        self._porous_region = porous_region
+
+        if species is not None:
+            self._y = species
+        else:
+            import numpy as np
+            self._y = np.empty((0,), dtype=object)
 
     def __call__(self, x_vec, gas_model):
         """Evaluate the wall+gas properties for porous materials.
@@ -93,10 +107,21 @@ def __call__(self, x_vec, gas_model):
         """
         actx = x_vec[0].array_context
         zeros = actx.np.zeros_like(x_vec[0])
+        ones = zeros + 1.0
         dim = x_vec.shape[0]
 
         temperature = self._temp + zeros
 
+        nspecies = len(self._y)
+        if nspecies > 0:
+            species_mass_frac = make_obj_array([self._y[i] + zeros
+                                                for i in range(nspecies)])
+        else:
+            species_mass_frac = self._y
+
+        porous_region = ones if self._porous_region is None else self._porous_region
+        wall_density = self._wall_density * porous_region
+
         species_mass_frac = self._y + zeros
 
         wall_density = self._wall_density + zeros

From 3aa5461af1829850bfe5d607407a33eaee310ca3 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 12 Mar 2024 16:14:38 -0500
Subject: [PATCH 17/48] docs for the materials

---
 mirgecom/materials/carbon_fiber.py |  1 +
 mirgecom/materials/initializer.py  |  2 +-
 mirgecom/materials/tacot.py        | 10 +++++++++-
 3 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index f6673282c..073407a89 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -109,6 +109,7 @@ class FiberEOS(PorousWallEOS):
         \tau = \frac{m}{m_0} = \frac{\rho_i \epsilon}{\rho_i \epsilon_0}
              = \frac{r^2}{r_0^2}
 
+    .. automethod:: __init__
     .. automethod:: void_fraction
     .. automethod:: enthalpy
     .. automethod:: heat_capacity
diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index b0977bdce..47e95eca8 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -56,7 +56,7 @@ def __call__(self, x_vec, wall_model):
 
 
 class PorousWallInitializer:
-    """Initializer for porous materials."""
+    """State initializer for porous materials."""
 
     def __init__(self, temperature, species, material_densities,
                  pressure=None, density=None):
diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index f17b95293..823fdbd5b 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -105,6 +105,7 @@ class TacotEOS(PorousWallEOS):
     yield the material properties. Polynomials were generated offline to avoid
     interpolation and they are not valid for temperatures above 3200K.
 
+    .. automethod:: __init__
     .. automethod:: void_fraction
     .. automethod:: enthalpy
     .. automethod:: heat_capacity
@@ -119,7 +120,14 @@ class TacotEOS(PorousWallEOS):
     def __init__(self, char_mass, virgin_mass):
         """Bulk density considering the porosity and intrinsic density.
 
-        The fiber and all resin components must be considered.
+        Parameters
+        ----------
+
+        virgin_mass: float
+            initial mass of the material. The fiber and all resin components
+            must be considered.
+        char_mass: float
+            final mass when the resin decomposition is complete.
         """
         self._char_mass = char_mass
         self._virgin_mass = virgin_mass

From 213ed3dff707bb205c5ffecb115512f892b06006 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 12 Mar 2024 16:23:15 -0500
Subject: [PATCH 18/48] docs for inviscid.py

---
 mirgecom/inviscid.py               | 11 +++++++++--
 mirgecom/materials/carbon_fiber.py |  1 -
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py
index da045614f..aa5d8420b 100644
--- a/mirgecom/inviscid.py
+++ b/mirgecom/inviscid.py
@@ -44,8 +44,6 @@
 """
 
 import numpy as np
-from arraycontext import outer
-from meshmode.dof_array import DOFArray
 from meshmode.discretization.connection import FACE_RESTR_ALL
 from grudge.dof_desc import (
     DD_VOLUME_ALL,
@@ -60,6 +58,9 @@
 from mirgecom.utils import normalize_boundaries
 from mirgecom.gas_model import PorousFlowFluidState
 
+from arraycontext import outer
+from meshmode.dof_array import DOFArray
+
 
 def inviscid_flux(state, gas_model=None):
     r"""Compute the inviscid flux vectors from fluid conserved vars *cv*.
@@ -81,6 +82,12 @@ def inviscid_flux(state, gas_model=None):
 
         Full fluid conserved and thermal state.
 
+    gas_model: :class:`~mirgecom.gas_model.GasModel`
+
+        Physical gas model including equation of state, transport,
+        and kinetic properties as required by fluid state.
+        Only necessary for porous media flows using the unified-domain.
+
     Returns
     -------
     :class:`~mirgecom.fluid.ConservedVars`
diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 2dc1c1e4f..4e4677ad7 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -131,7 +131,6 @@ def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass,
             geometrical dimension of the problem.
         virgin_mass: float
             initial mass of the material.
-
         char_mass: float
             final mass when the decomposition is complete.
         anisotropic_direction: int

From 9a041ba619664a039a1a9a0d32c99afaf590fdb0 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 13 Mar 2024 09:15:50 -0500
Subject: [PATCH 19/48] improve docs in initializer

---
 mirgecom/materials/initializer.py | 44 +++++++++++++++++++++----------
 1 file changed, 30 insertions(+), 14 deletions(-)

diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 47e95eca8..7d91847bb 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -34,7 +34,11 @@
 
 
 class SolidWallInitializer:
-    """Initializer for heat conduction only materials."""
+    """State initializer for wall models solving heat-diffusion equation.
+
+    This class computes the initial condition for either solid or porous
+    materials, and/or their combination, subject or not to ablation.
+    """
 
     def __init__(self, temperature):
         self._temp = temperature
@@ -48,6 +52,11 @@ def __call__(self, x_vec, wall_model):
             Nodal coordinates
         wall_model: :class:`mirgecom.wall_model.SolidWallModel`
             Equation of state class
+
+        Returns
+        -------
+        wv: :class:`mirgecom.wall_model.SolidWallConservedVars`
+            The conserved variables for heat-conduction only materials.
         """
         actx = x_vec[0].array_context
         mass = wall_model.density() + actx.np.zeros_like(x_vec[0])
@@ -56,7 +65,7 @@ def __call__(self, x_vec, wall_model):
 
 
 class PorousWallInitializer:
-    """State initializer for porous materials."""
+    """State initializer for porous materials in the unified-domain solver."""
 
     def __init__(self, temperature, species, material_densities,
                  pressure=None, density=None):
@@ -76,37 +85,44 @@ def __call__(self, dim, x_vec, gas_model):
             Nodal coordinates
         gas_model: :class:`mirgecom.wall_model.PorousFlowModel`
             Equation of state class
+
+        Returns
+        -------
+        cv: :class:`mirgecom.fluid.ConservedVars`
+            The conserved variables for porous-media flows. It depends on
+            both gas and porous material properties.
         """
         actx = x_vec[0].array_context
         zeros = actx.np.zeros_like(x_vec[0])
 
-        temperature = self._temp + zeros
-        species_mass_frac = self._y + zeros
+        # wall-only properties
         wall_density = self._wall_density + zeros
-
         tau = gas_model.decomposition_progress(wall_density)
+        eps_rho_solid = gas_model.solid_density(wall_density)
+
+        # coupled properties
+        temperature = self._temp + zeros
+        species_mass_frac = self._y + zeros
 
         eps_gas = gas_model.wall_eos.void_fraction(tau)
         if self._mass is None:
             pressure = self._pres + zeros
-            eps_rho_gas = eps_gas*gas_model.eos.get_density(pressure,
-                temperature, species_mass_frac)
+            eps_rho_gas = eps_gas*gas_model.eos.get_density(
+                pressure, temperature, species_mass_frac)
         else:
-            density = self._mass + zeros
-            eps_rho_gas = eps_gas*density
+            eps_rho_gas = eps_gas*self._mass
 
-        # internal energy (kinetic energy is neglected)
-        eps_rho_solid = gas_model.solid_density(wall_density)
+        # FIXME: for now, let's always start without velocity
+        momentum = make_obj_array([zeros for _ in range(dim)])
 
+        # internal energy (kinetic energy is absent in here)
         bulk_energy = (
             eps_rho_solid*gas_model.wall_eos.enthalpy(temperature, tau)
             + eps_rho_gas*gas_model.eos.get_internal_energy(temperature,
                                                             species_mass_frac)
         )
 
-        momentum = make_obj_array([zeros, zeros])
-
         species_mass = eps_rho_gas*species_mass_frac
 
         return make_conserved(dim=dim, mass=eps_rho_gas, energy=bulk_energy,
-            momentum=momentum, species_mass=species_mass)
+                              momentum=momentum, species_mass=species_mass)

From 3656771a79631e14616fec96b5092fa3a46ad4c1 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Thu, 14 Mar 2024 09:07:42 -0500
Subject: [PATCH 20/48] Remove dim

---
 mirgecom/materials/initializer.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 7d91847bb..c39973180 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -76,7 +76,7 @@ def __init__(self, temperature, species, material_densities,
         self._temp = temperature
         self._wall_density = material_densities
 
-    def __call__(self, dim, x_vec, gas_model):
+    def __call__(self, x_vec, gas_model):
         """Evaluate the wall+gas properties for porous materials.
 
         Parameters
@@ -94,6 +94,7 @@ def __call__(self, dim, x_vec, gas_model):
         """
         actx = x_vec[0].array_context
         zeros = actx.np.zeros_like(x_vec[0])
+        dim = x_vec.shape[0]
 
         # wall-only properties
         wall_density = self._wall_density + zeros
@@ -112,7 +113,7 @@ def __call__(self, dim, x_vec, gas_model):
         else:
             eps_rho_gas = eps_gas*self._mass
 
-        # FIXME: for now, let's always start without velocity
+        # FIXME: for now, let's always start with zero velocity
         momentum = make_obj_array([zeros for _ in range(dim)])
 
         # internal energy (kinetic energy is absent in here)

From 5fc7bb3ce71539ba963f73daf1289042c016ff35 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Sat, 23 Mar 2024 10:23:30 -0500
Subject: [PATCH 21/48] Is it necessary to have a specific dt calculation for
 porous flows?

---
 mirgecom/materials/carbon_fiber.py |  2 ++
 mirgecom/wall_model.py             | 47 ++++++++++++++++++++++++++----
 2 files changed, 43 insertions(+), 6 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 4e4677ad7..952dba28a 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -197,6 +197,7 @@ def thermal_conductivity(self, temperature, tau=None) -> DOFArray:
         kappa[self._anisotropic_dir] = kappa_k
 
         # account for fiber shrinkage via "tau"
+        # XXX check if here is the best place for timescale
         return kappa*tau*self._timescale
 
     # ~~~~~~~~ other properties
@@ -217,6 +218,7 @@ def permeability(self, tau: DOFArray) -> DOFArray:
 
     def emissivity(self, temperature=None, tau=None) -> DOFArray:
         """Emissivity for energy radiation."""
+        # XXX check if here is the best place for timescale
         return self._timescale * (
             + 2.26413679e-18*temperature**5 - 2.03008004e-14*temperature**4
             + 7.05300324e-11*temperature**3 - 1.22131715e-07*temperature**2
diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 757acf796..2e200391c 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -355,12 +355,9 @@ def thermal_conductivity(self, cv: ConservedVars, dv: GasDependentVars,
         .. math::
             \frac{\rho_s \kappa_s + \rho_g \kappa_g}{\rho_s + \rho_g}
         """
-        y_g = cv.mass/(cv.mass + wv.density)
-        y_s = 1.0 - y_g
         kappa_s = wall_eos.thermal_conductivity(dv.temperature, wv.tau)
         kappa_g = self.base_transport.thermal_conductivity(cv, dv, eos)
-
-        return y_s*kappa_s + y_g*kappa_g
+        return (wv.density*kappa_s + cv.mass*kappa_g)/(cv.mass + wv.density)
 
     def species_diffusivity(self, cv: ConservedVars, dv: GasDependentVars,
             wv: PorousWallVars, eos: GasEOS) -> DOFArray:
@@ -492,7 +489,7 @@ def get_pressure(self, cv: ConservedVars, wv: PorousWallVars,
         return cv.mass*self.eos.gas_const()*temperature/wv.void_fraction
 
     def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
-                 temperature: DOFArray) -> DOFArray:
+                        temperature: DOFArray) -> DOFArray:
         r"""Return the internal energy of the gas+solid material.
 
         .. math::
@@ -503,7 +500,7 @@ def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
                 + wv.density*self.wall_eos.enthalpy(temperature, wv.tau))
 
     def heat_capacity(self, cv: ConservedVars, wv: PorousWallVars,
-                 temperature: DOFArray) -> DOFArray:
+                      temperature: DOFArray) -> DOFArray:
         r"""Return the heat capacity of the gas+solid material.
 
         .. math::
@@ -511,3 +508,41 @@ def heat_capacity(self, cv: ConservedVars, wv: PorousWallVars,
         """
         return (cv.mass*self.eos.heat_capacity_cv(cv, temperature)
                 + wv.density*self.wall_eos.heat_capacity(temperature, wv.tau))
+
+#    def thermal_diffusivity(self, fluid_state) -> DOFArray:
+#        """Return the wall thermal diffusivity for all components."""
+#        tau = self.decomposition_progress(solid_state.cv.mass)
+#        mass = self.solid_density(solid_state.cv.mass)
+#        kappa = fluid_state.tv.thermal_conductivity
+#        cp = self.heat_capacity(fluid_state.cv, fluid_state.wv,
+#                                fluid_state.dv.temperature)
+#        return kappa/(mass * cp)
+
+
+def get_porous_flow_timestep(dcoll, gas_model, state, cfl, dd):
+    r"""Routine returns the the node-local maximum stable viscous timestep."""
+
+    from grudge.dt_utils import characteristic_lengthscales
+    from mirgecom.viscous import get_local_max_species_diffusivity
+    import grudge.op as op
+    from functools import reduce
+
+    actx = state.array_context
+
+    length_scales = characteristic_lengthscales(
+        state.array_context, dcoll, dd=dd)
+
+    kappa = reduce(actx.np.maximum, state.thermal_conductivity)
+    mass_cp = gas_model.heat_capacity(state.cv, state.wv, state.dv.temperature)
+
+    nu = state.viscosity / state.mass_density
+    alpha = kappa / ( mass_cp )
+    d_species_max = \
+        get_local_max_species_diffusivity(actx, state.species_diffusivity)
+    viscous_stuff = nu + alpha + d_species_max
+    vdt = op.elementwise_min(
+        dcoll, dd, 
+        length_scales / (state.wavespeed + ((viscous_stuff) / length_scales))
+    )
+
+    return cfl * vdt

From 531e2788528e6a5d75a4625bcf56bf1cd23d0265 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 25 Mar 2024 08:51:01 -0500
Subject: [PATCH 22/48] flake8

---
 mirgecom/wall_model.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 2e200391c..4d97b9386 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -536,12 +536,12 @@ def get_porous_flow_timestep(dcoll, gas_model, state, cfl, dd):
     mass_cp = gas_model.heat_capacity(state.cv, state.wv, state.dv.temperature)
 
     nu = state.viscosity / state.mass_density
-    alpha = kappa / ( mass_cp )
+    alpha = kappa / mass_cp
     d_species_max = \
         get_local_max_species_diffusivity(actx, state.species_diffusivity)
     viscous_stuff = nu + alpha + d_species_max
     vdt = op.elementwise_min(
-        dcoll, dd, 
+        dcoll, dd,
         length_scales / (state.wavespeed + ((viscous_stuff) / length_scales))
     )
 

From c11bc5f1689b9300d38d4857c1f418c07eb2ce69 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 25 Mar 2024 08:53:00 -0500
Subject: [PATCH 23/48] pydocstyle

---
 mirgecom/wall_model.py | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 4d97b9386..34e2ca057 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -86,6 +86,11 @@
 )
 from mirgecom.transport import GasTransportVars, TransportModel
 
+from grudge.dt_utils import characteristic_lengthscales
+from mirgecom.viscous import get_local_max_species_diffusivity
+import grudge.op as op
+from functools import reduce
+
 
 @with_container_arithmetic(bcast_obj_array=False,
                            bcast_container_types=(DOFArray, np.ndarray),
@@ -509,24 +514,18 @@ def heat_capacity(self, cv: ConservedVars, wv: PorousWallVars,
         return (cv.mass*self.eos.heat_capacity_cv(cv, temperature)
                 + wv.density*self.wall_eos.heat_capacity(temperature, wv.tau))
 
-#    def thermal_diffusivity(self, fluid_state) -> DOFArray:
-#        """Return the wall thermal diffusivity for all components."""
-#        tau = self.decomposition_progress(solid_state.cv.mass)
-#        mass = self.solid_density(solid_state.cv.mass)
-#        kappa = fluid_state.tv.thermal_conductivity
-#        cp = self.heat_capacity(fluid_state.cv, fluid_state.wv,
-#                                fluid_state.dv.temperature)
-#        return kappa/(mass * cp)
+    # def thermal_diffusivity(self, fluid_state) -> DOFArray:
+    #     """Return the wall thermal diffusivity for all components."""
+    #     tau = self.decomposition_progress(solid_state.cv.mass)
+    #     mass = self.solid_density(solid_state.cv.mass)
+    #     kappa = fluid_state.tv.thermal_conductivity
+    #     cp = self.heat_capacity(fluid_state.cv, fluid_state.wv,
+    #                             fluid_state.dv.temperature)
+    #     return kappa/(mass * cp)
 
 
 def get_porous_flow_timestep(dcoll, gas_model, state, cfl, dd):
     r"""Routine returns the the node-local maximum stable viscous timestep."""
-
-    from grudge.dt_utils import characteristic_lengthscales
-    from mirgecom.viscous import get_local_max_species_diffusivity
-    import grudge.op as op
-    from functools import reduce
-
     actx = state.array_context
 
     length_scales = characteristic_lengthscales(

From dcc81e8cf3689f49c0f3d68c87b9bb631067966a Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 26 Mar 2024 08:19:00 -0500
Subject: [PATCH 24/48] Keep Y2 oxidation model for now. Postpone it to Darcy
 flow.

---
 mirgecom/materials/carbon_fiber.py | 58 ++++++++++++++++++------------
 1 file changed, 36 insertions(+), 22 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 073407a89..6845256c8 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -54,33 +54,38 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
         raise NotImplementedError()
 
 
-class OxidationModel(Oxidation):
-    """Evaluate the source terms for the carbon fiber oxidation.
+# TODO per MTC review, can we generalize the oxidation model?
+# should we keep this in the driver?
+class Y2_Oxidation_Model(Oxidation):  # noqa N801
+    """Evaluate the source terms for the Y2 model of carbon fiber oxidation.
 
-    The user must specify in the driver the functions for oxidation.
-    (Tentatively) Generalizing this class makes it easier for adding new,
-    more complex models and for UQ runs.
-
-    .. automethod:: effective_surface_area_fiber
+    .. automethod:: puma_effective_surface_area
     .. automethod:: get_source_terms
     """
 
-    def __init__(self, surface_area_func, oxidation_func):
-        self._surface_func = surface_area_func
-        self._oxidation_func = oxidation_func
+    def puma_effective_surface_area(self, progress) -> DOFArray:
+        """Polynomial fit based on PUMA data.
 
-    def effective_surface_area_fiber(self, progress: DOFArray) -> DOFArray:
-        r"""Evaluate the effective surface of the fibers.
+        Parameters
+        ----------
+        progress: meshmode.dof_array.DOFArray
+            the rate of decomposition of the fibers
+        """
+        # Original fit function: -1.1012e5*x**2 - 0.0646e5*x + 1.1794e5
+        # Rescale by x==0 value and rearrange
+        return 1.1794e5*(1.0 - 0.0547736137*progress - 0.9336950992*progress**2)
+
+    def _get_wall_effective_surface_area_fiber(self, progress) -> DOFArray:
+        """Evaluate the effective surface of the fibers.
 
         Parameters
         ----------
         progress: meshmode.dof_array.DOFArray
-            the rate of decomposition of the fibers, defined as 1.0 - $\tau$.
+            the rate of decomposition of the fibers
         """
-        return self._surface_func(progress)
+        return self.puma_effective_surface_area(progress)
 
-    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: Optional[DOFArray] = None) -> DOFArray:  # noqa N803
+    def get_source_terms(self, temperature, tau, rhoY_o2) -> DOFArray:  # noqa N803
         """Return the effective source terms for the oxidation.
 
         Parameters
@@ -88,14 +93,23 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
         temperature: meshmode.dof_array.DOFArray
         tau: meshmode.dof_array.DOFArray
             the progress ratio of the oxidation
-        rhoY_o2: meshmode.dof_array.DOFArray
+        ox_mass: meshmode.dof_array.DOFArray
             the mass fraction of oxygen
         """
-        area = self.effective_surface_area_fiber(1.0-tau)
-        if rhoY_o2:
-            return self._oxidation_func(temperature=temperature, fiber_area=area,
-                                        rhoY_o2=rhoY_o2)
-        return self._oxidation_func(temperature=temperature, fiber_area=area)
+        actx = temperature.array_context
+
+        mw_o = 15.999
+        mw_o2 = mw_o*2
+        mw_co = 28.010
+        univ_gas_const = 8314.46261815324
+
+        eff_surf_area = self._get_wall_effective_surface_area_fiber(1.0-tau)
+        alpha = (
+            (0.00143+0.01*actx.np.exp(-1450.0/temperature))
+            / (1.0+0.0002*actx.np.exp(13000.0/temperature)))
+        k = alpha*actx.np.sqrt(
+            (univ_gas_const*temperature)/(2.0*np.pi*mw_o2))
+        return (mw_co/mw_o2 + mw_o/mw_o2 - 1)*rhoY_o2*k*eff_surf_area
 
 
 class FiberEOS(PorousWallEOS):

From 94873e352eb0014cf9b40116bac6878e5c48f847 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 26 Mar 2024 09:01:09 -0500
Subject: [PATCH 25/48] Improve docs, mypy in carbon_fiber and wall_model

---
 mirgecom/materials/carbon_fiber.py | 38 ++++++++++++++++--------------
 mirgecom/wall_model.py             | 19 +++++++--------
 2 files changed, 29 insertions(+), 28 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 6845256c8..cae6f0d12 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -30,6 +30,7 @@
 
 from typing import Optional
 from abc import abstractmethod
+import numpy as np
 from meshmode.dof_array import DOFArray
 from mirgecom.wall_model import PorousWallEOS
 from pytools.obj_array import make_obj_array
@@ -38,24 +39,19 @@
 class Oxidation:
     """Abstract interface for wall oxidation model.
 
-    .. automethod:: effective_surface_area_fiber
     .. automethod:: get_source_terms
     """
 
-    @abstractmethod
-    def effective_surface_area_fiber(self, progress: DOFArray) -> DOFArray:
-        r"""Evaluate the effective surface of the fibers."""
-        raise NotImplementedError()
-
     @abstractmethod
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: Optional[DOFArray] = None) -> DOFArray:  # noqa N803
+            rhoY_o2: DOFArray) -> DOFArray:
         r"""Source terms of fiber oxidation."""
         raise NotImplementedError()
 
 
 # TODO per MTC review, can we generalize the oxidation model?
 # should we keep this in the driver?
+# https://github.com/illinois-ceesd/mirgecom/pull/875#discussion_r1261414281
 class Y2_Oxidation_Model(Oxidation):  # noqa N801
     """Evaluate the source terms for the Y2 model of carbon fiber oxidation.
 
@@ -75,7 +71,7 @@ def puma_effective_surface_area(self, progress) -> DOFArray:
         # Rescale by x==0 value and rearrange
         return 1.1794e5*(1.0 - 0.0547736137*progress - 0.9336950992*progress**2)
 
-    def _get_wall_effective_surface_area_fiber(self, progress) -> DOFArray:
+    def _get_wall_effective_surface_area_fiber(self, progress: DOFArray) -> DOFArray:
         """Evaluate the effective surface of the fibers.
 
         Parameters
@@ -85,14 +81,15 @@ def _get_wall_effective_surface_area_fiber(self, progress) -> DOFArray:
         """
         return self.puma_effective_surface_area(progress)
 
-    def get_source_terms(self, temperature, tau, rhoY_o2) -> DOFArray:  # noqa N803
+    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
+            rhoY_o2: DOFArray) -> DOFArray:
         """Return the effective source terms for the oxidation.
 
         Parameters
         ----------
         temperature: meshmode.dof_array.DOFArray
         tau: meshmode.dof_array.DOFArray
-            the progress ratio of the oxidation
+            the progress ratio of the oxidation: 1 for virgin, 0 for fully oxidized
         ox_mass: meshmode.dof_array.DOFArray
             the mass fraction of oxygen
         """
@@ -142,13 +139,13 @@ def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass):
         ----------
         dim: int
             geometrical dimension of the problem.
-        virgin_mass: float
-            initial mass of the material.
-        char_mass: float
-            final mass when the decomposition is complete.
         anisotropic_direction: int
             For orthotropic materials, this indicates the normal direction
             where the properties are different than in-plane.
+        char_mass: float
+            final mass when the decomposition is complete.
+        virgin_mass: float
+            initial mass of the material.
         """
         self._char_mass = char_mass
         self._virgin_mass = virgin_mass
@@ -186,7 +183,7 @@ def heat_capacity(self, temperature, tau=None) -> DOFArray:
             - 3.62422269e+02)
 
     # ~~~~~~~~ fiber conductivity
-    def thermal_conductivity(self, temperature, tau=None) -> DOFArray:
+    def thermal_conductivity(self, temperature, tau=None) -> np.ndarray:
         r"""Evaluate the thermal conductivity $\kappa$ of the fibers.
 
         It accounts for anisotropy and oxidation progress.
@@ -214,7 +211,7 @@ def volume_fraction(self, tau: DOFArray) -> DOFArray:
         r"""Fraction $\phi$ occupied by the solid."""
         return 0.12*tau
 
-    def permeability(self, tau: DOFArray) -> DOFArray:
+    def permeability(self, tau: DOFArray) -> np.ndarray:
         r"""Permeability $K$ of the porous material."""
         # FIXME find a relation to make it change as a function of "tau"
         # TODO: the relation depends on the coupling model. Postpone it for now.
@@ -224,7 +221,8 @@ def permeability(self, tau: DOFArray) -> DOFArray:
         permeability[self._anisotropic_dir] = 2.62e-11 + actx.np.zeros_like(tau)
         return permeability
 
-    def emissivity(self, temperature=None, tau=None) -> DOFArray:
+    def emissivity(self, temperature: DOFArray,  # type: ignore[override]
+            tau: Optional[DOFArray] = None) -> DOFArray:
         """Emissivity for energy radiation."""
         return (
             + 2.26413679e-18*temperature**5 - 2.03008004e-14*temperature**4
@@ -232,7 +230,11 @@ def emissivity(self, temperature=None, tau=None) -> DOFArray:
             + 1.21137817e-04*temperature**1 + 8.66656964e-01)
 
     def tortuosity(self, tau: DOFArray) -> DOFArray:
-        r"""Tortuosity $\eta$ affects the species diffusivity."""
+        r"""Tortuosity $\eta$ affects the species diffusivity.
+
+        .. math:
+            D_{eff} = \frac{D_i^(m)}{\eta}
+        """
         return 1.1*tau + 1.0*(1.0 - tau)
 
     def decomposition_progress(self, mass: DOFArray) -> DOFArray:
diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 20d5ac4db..7011ad3ca 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -214,7 +214,7 @@ class PorousWallVars:
     material_densities: Union[DOFArray, np.ndarray]
     tau: DOFArray
     void_fraction: DOFArray
-    permeability: DOFArray
+    permeability: Union[DOFArray, np.ndarray]
     tortuosity: DOFArray
     density: DOFArray
 
@@ -265,18 +265,19 @@ def volume_fraction(self, tau: DOFArray) -> DOFArray:
         raise NotImplementedError()
 
     @abstractmethod
-    def permeability(self, tau: DOFArray) -> DOFArray:
+    def permeability(self, tau: DOFArray) -> Union[np.ndarray, DOFArray]:
         r"""Permeability $K$ of the porous material."""
         raise NotImplementedError()
 
     @abstractmethod
-    def emissivity(self, temperature=None, tau=None) -> DOFArray:
+    def emissivity(self, temperature: Optional[DOFArray] = None,
+            tau: Optional[DOFArray] = None) -> DOFArray:
         """Emissivity for energy radiation."""
         raise NotImplementedError()
 
     @abstractmethod
     def tortuosity(self, tau: DOFArray) -> DOFArray:
-        """Tortuosity of the porous material."""
+        r"""Tortuosity $\eta$ of the porous material."""
         raise NotImplementedError()
 
     @abstractmethod
@@ -303,7 +304,7 @@ class PorousWallTransport:
     """
 
     def __init__(self, base_transport: TransportModel):
-        """Initialize transport model."""
+        """Initialize base transport model for fluid-only."""
         self.base_transport = base_transport
 
     def bulk_viscosity(self, cv: ConservedVars, dv: GasDependentVars,
@@ -324,7 +325,8 @@ def viscosity(self, cv: ConservedVars, dv: GasDependentVars,
             self.base_transport.viscosity(cv, dv, eos))
 
     def thermal_conductivity(self, cv: ConservedVars, dv: GasDependentVars,
-            wv: PorousWallVars, eos: GasEOS, wall_eos: PorousWallEOS) -> DOFArray:
+            wv: PorousWallVars, eos: GasEOS,
+            wall_eos: PorousWallEOS) -> Union[np.ndarray, DOFArray]:
         r"""Return the effective thermal conductivity of the gas+solid.
 
         It is a function of temperature and degradation progress. As the
@@ -336,12 +338,9 @@ def thermal_conductivity(self, cv: ConservedVars, dv: GasDependentVars,
         .. math::
             \frac{\rho_s \kappa_s + \rho_g \kappa_g}{\rho_s + \rho_g}
         """
-        y_g = cv.mass/(cv.mass + wv.density)
-        y_s = 1.0 - y_g
         kappa_s = wall_eos.thermal_conductivity(dv.temperature, wv.tau)
         kappa_g = self.base_transport.thermal_conductivity(cv, dv, eos)
-
-        return y_s*kappa_s + y_g*kappa_g
+        return (wv.density*kappa_s + cv.mass*kappa_g)/(cv.mass + wv.density)
 
     def species_diffusivity(self, cv: ConservedVars, dv: GasDependentVars,
             wv: PorousWallVars, eos: GasEOS) -> DOFArray:

From 6d4163c142f6c1ec8967abdc4b8ac6c778cc6219 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 26 Mar 2024 09:05:59 -0500
Subject: [PATCH 26/48] flake8

---
 mirgecom/materials/carbon_fiber.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index cae6f0d12..fd4b727ee 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -43,8 +43,8 @@ class Oxidation:
     """
 
     @abstractmethod
-    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: DOFArray) -> DOFArray:
+    def get_source_terms(self, temperature: DOFArray,
+            tau: DOFArray, rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
         r"""Source terms of fiber oxidation."""
         raise NotImplementedError()
 
@@ -82,7 +82,7 @@ def _get_wall_effective_surface_area_fiber(self, progress: DOFArray) -> DOFArray
         return self.puma_effective_surface_area(progress)
 
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: DOFArray) -> DOFArray:
+            rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
         """Return the effective source terms for the oxidation.
 
         Parameters

From fdd7fa8532dca6ff72b8e6f3d99396e3858616ae Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 26 Mar 2024 09:53:15 -0500
Subject: [PATCH 27/48] Use tau instead of progress for Y2 model + docs, mypy

---
 mirgecom/materials/carbon_fiber.py            | 34 +++++++------------
 .../materials/prescribed_porous_material.py   | 12 +++----
 2 files changed, 16 insertions(+), 30 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index fd4b727ee..c4e72636b 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -1,6 +1,7 @@
 r""":mod:`mirgecom.materials.carbon_fiber` evaluate carbon fiber data.
 
 .. autoclass:: Oxidation
+.. autoclass:: Y2_Oxidation_Model
 .. autoclass:: FiberEOS
 """
 
@@ -59,38 +60,27 @@ class Y2_Oxidation_Model(Oxidation):  # noqa N801
     .. automethod:: get_source_terms
     """
 
-    def puma_effective_surface_area(self, progress) -> DOFArray:
-        """Polynomial fit based on PUMA data.
-
-        Parameters
-        ----------
-        progress: meshmode.dof_array.DOFArray
-            the rate of decomposition of the fibers
-        """
+    def puma_effective_surface_area(self, tau: DOFArray) -> DOFArray:
+        """Polynomial fit based on PUMA data."""
         # Original fit function: -1.1012e5*x**2 - 0.0646e5*x + 1.1794e5
         # Rescale by x==0 value and rearrange
+        progress = 1.0-tau
         return 1.1794e5*(1.0 - 0.0547736137*progress - 0.9336950992*progress**2)
 
-    def _get_wall_effective_surface_area_fiber(self, progress: DOFArray) -> DOFArray:
-        """Evaluate the effective surface of the fibers.
-
-        Parameters
-        ----------
-        progress: meshmode.dof_array.DOFArray
-            the rate of decomposition of the fibers
-        """
-        return self.puma_effective_surface_area(progress)
+    def _get_wall_effective_surface_area_fiber(self, tau: DOFArray) -> DOFArray:
+        """Evaluate the effective surface of the fibers."""
+        return self.puma_effective_surface_area(tau)
 
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
             rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
-        """Return the effective source terms for the oxidation.
+        """Return the effective source terms for fiber oxidation.
 
         Parameters
         ----------
-        temperature: meshmode.dof_array.DOFArray
-        tau: meshmode.dof_array.DOFArray
+        temperature:
+        tau:
             the progress ratio of the oxidation: 1 for virgin, 0 for fully oxidized
-        ox_mass: meshmode.dof_array.DOFArray
+        rhoY_o2:
             the mass fraction of oxygen
         """
         actx = temperature.array_context
@@ -100,7 +90,7 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
         mw_co = 28.010
         univ_gas_const = 8314.46261815324
 
-        eff_surf_area = self._get_wall_effective_surface_area_fiber(1.0-tau)
+        eff_surf_area = self._get_wall_effective_surface_area_fiber(tau)
         alpha = (
             (0.00143+0.01*actx.np.exp(-1450.0/temperature))
             / (1.0+0.0002*actx.np.exp(13000.0/temperature)))
diff --git a/mirgecom/materials/prescribed_porous_material.py b/mirgecom/materials/prescribed_porous_material.py
index e63dcc358..f69673637 100644
--- a/mirgecom/materials/prescribed_porous_material.py
+++ b/mirgecom/materials/prescribed_porous_material.py
@@ -71,21 +71,17 @@ def void_fraction(self, tau: DOFArray) -> DOFArray:
         return 1.0 - self.volume_fraction(tau)
 
     def enthalpy(self, temperature: DOFArray, tau: Optional[DOFArray]) -> DOFArray:
-        r"""Evaluate the solid enthalpy $h_s$ of the fibers."""
+        r"""Evaluate the solid enthalpy $h_s$."""
         return self._enthalpy_func(temperature)
 
     def heat_capacity(self, temperature: DOFArray,
                       tau: Optional[DOFArray]) -> DOFArray:
-        r"""Evaluate the heat capacity $C_{p_s}$ of the fibers."""
+        r"""Evaluate the heat capacity $C_{p_s}$."""
         return self._heat_capacity_func(temperature)
 
     def thermal_conductivity(self, temperature: DOFArray,
                                    tau: DOFArray) -> DOFArray:
-        r"""Evaluate the thermal conductivity $\kappa$ of the fibers.
-
-        It employs a rescaling of the experimental data based on the fiber
-        shrinkage during the oxidation.
-        """
+        r"""Evaluate the thermal conductivity $\kappa$"""
         return self._thermal_conductivity_func(temperature)
 
     def volume_fraction(self, tau: DOFArray) -> DOFArray:
@@ -109,5 +105,5 @@ def tortuosity(self, tau: DOFArray) -> DOFArray:
         return self._tortuosity_func(tau) + actx.np.zeros_like(tau)
 
     def decomposition_progress(self, mass: DOFArray) -> DOFArray:
-        r"""Evaluate the progress rate $\tau$."""
+        r"""Evaluate the progress ratio $\tau$."""
         return self._decomposition_progress_func(mass)

From a26e5d54cdc81d72420d3561e3976f874c399c78 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 26 Mar 2024 10:05:35 -0500
Subject: [PATCH 28/48] Organize oxidation models; docs

---
 mirgecom/materials/carbon_fiber.py | 124 ++++++++++++++++++++++++++---
 1 file changed, 113 insertions(+), 11 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 793708ed1..906af0b5a 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -2,6 +2,8 @@
 
 .. autoclass:: Oxidation
 .. autoclass:: Y2_Oxidation_Model
+.. autoclass:: Y3_Oxidation_Model
+.. autoclass:: OxidationModel
 .. autoclass:: FiberEOS
 """
 
@@ -50,27 +52,20 @@ def get_source_terms(self, temperature: DOFArray,
         raise NotImplementedError()
 
 
-# TODO per MTC review, can we generalize the oxidation model?
-# should we keep this in the driver?
-# https://github.com/illinois-ceesd/mirgecom/pull/875#discussion_r1261414281
 class Y2_Oxidation_Model(Oxidation):  # noqa N801
-    """Evaluate the source terms for the Y2 model of carbon fiber oxidation.
+    r"""Evaluate the source terms for carbon fiber oxidation in Y2 model.
 
-    .. automethod:: puma_effective_surface_area
     .. automethod:: get_source_terms
     """
 
-    def puma_effective_surface_area(self, tau: DOFArray) -> DOFArray:
-        """Polynomial fit based on PUMA data."""
+    def _get_wall_effective_surface_area_fiber(self, tau: DOFArray) -> DOFArray:
+        """Evaluate the effective surface of the fibers."""
+        # Polynomial fit based on PUMA data:
         # Original fit function: -1.1012e5*x**2 - 0.0646e5*x + 1.1794e5
         # Rescale by x==0 value and rearrange
         progress = 1.0-tau
         return 1.1794e5*(1.0 - 0.0547736137*progress - 0.9336950992*progress**2)
 
-    def _get_wall_effective_surface_area_fiber(self, tau: DOFArray) -> DOFArray:
-        """Evaluate the effective surface of the fibers."""
-        return self.puma_effective_surface_area(tau)
-
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
             rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
         """Return the effective source terms for fiber oxidation.
@@ -99,6 +94,113 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
         return (mw_co/mw_o2 + mw_o/mw_o2 - 1)*rhoY_o2*k*eff_surf_area
 
 
+class Y3_Oxidation_Model(Oxidation):  # noqa N801
+    r"""Evaluate the source terms for carbon fiber oxidation in Y3 model.
+
+    Follows ``A. Martin, AIAA 2013-2636'', using a single reaction given by
+
+    .. math::
+        C_{(s)} + O_2 \to CO_2
+
+    .. automethod:: get_source_terms
+    """
+
+    def __init__(self, wall_material):
+        self._material = wall_material
+
+    def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
+        r"""Evaluate the effective surface of the fibers.
+
+        The fiber radius as a function of mass loss $\tau$ is given by
+
+        .. math::
+            \tau = \frac{m}{m_0} = \frac{\pi r^2/L}{\pi r_0^2/L} = \frac{r^2}{r_0^2}
+        """
+        actx = tau.array_context
+
+        original_fiber_radius = 5e-6  # half the diameter
+        fiber_radius = original_fiber_radius*actx.np.sqrt(tau)
+
+        epsilon_0 = self._material.volume_fraction(tau=1.0)
+
+        return 2.0*epsilon_0/original_fiber_radius**2*fiber_radius
+
+    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
+            rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+        r"""Return the effective source terms for the oxidation.
+
+        Parameters
+        ----------
+        temperature:
+        tau:
+            the progress ratio of the oxidation
+        rhoY_o2:
+            the mass fraction of oxygen
+
+        Returns
+        -------
+            The tuple (\omega_{C}, \omega_{O_2}, \omega_{CO_2})
+        """
+        actx = temperature.array_context
+
+        mw_c = 12.011
+        mw_o = 15.999
+        mw_o2 = mw_o*2
+        mw_co2 = 44.010
+        univ_gas_const = 8.31446261815324  # J/(K-mol)
+
+        eff_surf_area = self._get_wall_effective_surface_area_fiber(tau)
+
+        k_f = 1.0e5*actx.np.exp(-120000.0/(univ_gas_const*temperature))
+
+        m_dot_c = - rhoY_o2/mw_o2 * mw_c * eff_surf_area * k_f
+        m_dot_o2 = - rhoY_o2/mw_o2 * mw_o2 * eff_surf_area * k_f
+        m_dot_co2 = + rhoY_o2/mw_o2 * mw_co2 * eff_surf_area * k_f
+
+        return m_dot_c, m_dot_o2, m_dot_co2
+
+
+class OxidationModel(Oxidation):
+    """Evaluate the source terms for the carbon fiber oxidation.
+
+    The user must specify in the driver the functions for oxidation.
+    (Tentatively) Generalizing this class makes it easier for adding new,
+    more complex models and for UQ runs.
+
+    .. automethod:: get_source_terms
+    """
+
+    def __init__(self, surface_area_func, oxidation_func):
+        """Initialize the general oxidation class.
+
+        Parameters
+        ----------
+        surface_area_func:
+            Function prescribing how the fiber area changes during oxidation.
+        oxidation_func:
+            Reaction rate for the oxidation model.
+        """
+        self._surface_func = surface_area_func
+        self._oxidation_func = oxidation_func
+
+    # TODO we potentially have to include atomic oxygen as well
+    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
+            rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+        """Return the effective source terms for the oxidation.
+
+        Parameters
+        ----------
+        temperature:
+        tau:
+            the progress ratio of the oxidation
+        rhoY_o2:
+            the mass fraction of oxygen
+        """
+        area = self._surface_func(tau)
+        return self._oxidation_func(temperature=temperature, fiber_area=area,
+                                    rhoY_o2=rhoY_o2)
+
+
 class FiberEOS(PorousWallEOS):
     r"""Evaluate the properties of the solid state containing only fibers.
 

From 73b04dfa4b1553b265b41612b44015bb5ae0355d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 26 Mar 2024 15:04:39 -0500
Subject: [PATCH 29/48] Add velocity to porous flow initializer; fix broken
 documentation

---
 doc/misc.rst                       |  2 ++
 mirgecom/materials/carbon_fiber.py | 11 ++++++-----
 mirgecom/materials/initializer.py  | 18 ++++++++++++------
 3 files changed, 20 insertions(+), 11 deletions(-)

diff --git a/doc/misc.rst b/doc/misc.rst
index df4c005a0..33d1e466d 100644
--- a/doc/misc.rst
+++ b/doc/misc.rst
@@ -96,3 +96,5 @@ References
 .. [Ern_2008] Daniele A. Di Pietro, Alexandre Ern, Jean-Luc Guermond, Discontinuous Galerkin Methods for \
    Anisotropic Semidefinite Diffusion with Advection, SIAM Journal on Numerical Analysis 46 2 \
    `(DOI) <https://www.jstor.org/stable/40233233?seq=12>`__
+.. [Martin_2013] Alexandre Martin, Volume averaged modeling of the oxidation of porous carbon fiber material,
+    AIAA Thermophysics conference 2013-2636
diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 906af0b5a..3d592f293 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -31,7 +31,7 @@
 THE SOFTWARE.
 """
 
-from typing import Optional
+from typing import Optional, Union
 from abc import abstractmethod
 import numpy as np
 from meshmode.dof_array import DOFArray
@@ -47,7 +47,7 @@ class Oxidation:
 
     @abstractmethod
     def get_source_terms(self, temperature: DOFArray,
-            tau: DOFArray, rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+            tau: DOFArray, rhoY_o2: DOFArray) -> Union[DOFArray, tuple]:  # noqa N803
         r"""Source terms of fiber oxidation."""
         raise NotImplementedError()
 
@@ -97,7 +97,7 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
 class Y3_Oxidation_Model(Oxidation):  # noqa N801
     r"""Evaluate the source terms for carbon fiber oxidation in Y3 model.
 
-    Follows ``A. Martin, AIAA 2013-2636'', using a single reaction given by
+    Follows [Martin_2013]_ by using a single reaction given by
 
     .. math::
         C_{(s)} + O_2 \to CO_2
@@ -126,7 +126,7 @@ def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
         return 2.0*epsilon_0/original_fiber_radius**2*fiber_radius
 
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+            rhoY_o2: DOFArray):  # noqa N803
         r"""Return the effective source terms for the oxidation.
 
         Parameters
@@ -139,7 +139,8 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
 
         Returns
         -------
-            The tuple (\omega_{C}, \omega_{O_2}, \omega_{CO_2})
+        sources: tuple
+            the tuple ($\omega_{C}$, $\omega_{O_2}$, $\omega_{CO_2}$)
         """
         actx = temperature.array_context
 
diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 502c48107..1cfff75db 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -28,6 +28,7 @@
 THE SOFTWARE.
 """
 
+import numpy as np
 from pytools.obj_array import make_obj_array
 from mirgecom.fluid import make_conserved
 from mirgecom.wall_model import SolidWallConservedVars
@@ -74,8 +75,8 @@ def __call__(self, x_vec, wall_model):
 class PorousWallInitializer:
     """State initializer for porous materials in the unified-domain solver."""
 
-    def __init__(self, temperature, material_densities, species=None,
-                 pressure=None, density=None, porous_region=None):
+    def __init__(self, *, temperature, material_densities, species=None,
+                 velocity=None, pressure=None, density=None, porous_region=None):
         """Initialize the object for porous materials.
 
         Parameters
@@ -85,6 +86,7 @@ def __init__(self, temperature, material_densities, species=None,
             portions of the flow. Only used for unified-domain solver without
             explicit coupling.
         """
+        self._velocity = velocity
         self._pres = pressure
         self._mass = density
         self._temp = temperature
@@ -120,7 +122,10 @@ def __call__(self, x_vec, gas_model):
         ones = zeros + 1.0
         dim = x_vec.shape[0]
 
-        porous_region = ones if self._porous_region is None else self._porous_region
+        if self._porous_region is not None:
+            porous_region = self._porous_region + zeros
+        else:
+            porous_region = ones
 
         # wall-only properties
         wall_density = self._wall_density * porous_region
@@ -145,14 +150,15 @@ def __call__(self, x_vec, gas_model):
         else:
             eps_rho_gas = eps_gas*self._mass
 
-        # FIXME: for now, let's always start with zero velocity
-        momentum = make_obj_array([zeros for _ in range(dim)])
+        # start with zero velocity inside the material
+        #momentum = make_obj_array([zeros for _ in range(dim)])
+        momentum = (1.0-self._porous_region) * (eps_gas*self._velocity)
 
-        # internal energy (kinetic energy is absent in here)
         bulk_energy = (
             eps_rho_solid*gas_model.wall_eos.enthalpy(temperature, tau)
             + eps_rho_gas*gas_model.eos.get_internal_energy(temperature,
                                                             species_mass_frac)
+            + 0.5*np.dot(momentum, momentum)/eps_rho_gas
         )
 
         species_mass = eps_rho_gas*species_mass_frac

From e33e91c9fc6a88bd4278929048d1ab021cd4bf25 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 27 Mar 2024 09:06:48 -0500
Subject: [PATCH 30/48] Optional velocity to porous material

---
 mirgecom/materials/initializer.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 1cfff75db..79c98574b 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -151,8 +151,10 @@ def __call__(self, x_vec, gas_model):
             eps_rho_gas = eps_gas*self._mass
 
         # start with zero velocity inside the material
-        #momentum = make_obj_array([zeros for _ in range(dim)])
-        momentum = (1.0-self._porous_region) * (eps_gas*self._velocity)
+        if self._velocity is not None:
+            momentum = (1.0-self._porous_region) * (eps_gas*self._velocity)
+        else:
+            momentum = make_obj_array([zeros for _ in range(dim)])
 
         bulk_energy = (
             eps_rho_solid*gas_model.wall_eos.enthalpy(temperature, tau)

From 03f47dc8703602e0b81e3efc5f77cb535e96d417 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 27 Mar 2024 09:35:13 -0500
Subject: [PATCH 31/48] Pass arguments to TACOT degradation

---
 mirgecom/materials/tacot.py | 41 +++++++++++++++++++++++++++----------
 1 file changed, 30 insertions(+), 11 deletions(-)

diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index 823fdbd5b..52063ba74 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -57,9 +57,27 @@ class Pyrolysis:
     .. automethod:: get_source_terms
     """
 
-    def __init__(self):
-        """Temperature in which each reaction starts."""
+    def __init__(
+            self, *, virgin_mass=120.0, char_mass=60.0, fiber_mass=160.0,
+            pre_exponential=(12000.0, 4.48e9)):
+        """Initialize TACOT parameters."""
         self._Tcrit = np.array([333.3, 555.6])
+        self._virgin_mass = virgin_mass
+        self._char_mass = char_mass
+        self._fiber_mass = fiber_mass
+
+        if len(pre_exponential) != 2:
+            raise ValueError("TACOT degradation model requires 2 pre-exponentials.")
+        self._pre_exp = np.array(pre_exponential)
+
+    def get_tacot_parameters(self):
+        """Return the parameters of TACOT decomposition for inspection."""
+        return {"virgin mass": self._virgin_mass,
+                "char_mass": self._char_mass,
+                "fiber_mass": self._fiber_mass,
+                "reaction_weights": np.array([self._virgin_mass*0.75,
+                                              self._virgin_mass*0.25]),
+                "pre_exponential": self._pre_exp}
 
     def get_source_terms(self, temperature, chi):
         r"""Return the source terms of pyrolysis decomposition for TACOT.
@@ -81,20 +99,21 @@ def get_source_terms(self, temperature, chi):
         """
         actx = temperature.array_context
 
-        # The density parameters are hard-coded for TACOT. They depend on the
-        # virgin and char volume fraction.
+        # The density parameters are hard-coded for TACOT, depending on
+        # virgin and char volume fractions.
+        w1 = self._virgin_mass*0.25*(chi[0]/self._virgin_mass*0.25)**3
+        w2 = self._virgin_mass*0.75*(chi[1]/self._virgin_mass*0.75 - 2./3.)**3
+
         return make_obj_array([
             # reaction 1
             actx.np.where(actx.np.less(temperature, self._Tcrit[0]),
-                0.0, (
-                    -(30.*((chi[0] - 0.00)/30.)**3)*12000.
-                    * actx.np.exp(-8556.000/temperature))),
+                0.0, (-w1 * self._pre_exp[0] * actx.np.exp(-8556.000/temperature))
+                ),
             # reaction 2
             actx.np.where(actx.np.less(temperature, self._Tcrit[1]),
-                0.0, (
-                    -(90.*((chi[1] - 60.0)/90.)**3)*4.48e9
-                    * actx.np.exp(-20444.44/temperature))),
-            # fiber oxidation: include in the RHS but dont do anything with it.
+                0.0, (-w2 * self._pre_exp[1] * actx.np.exp(-20444.44/temperature))
+                ),
+            # fiber oxidation: include in the RHS but don't do anything with it.
             actx.np.zeros_like(temperature)])
 
 

From 1dda90e0190458b4def1cc92e177ec2f0cee2a7d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 27 Mar 2024 10:34:10 -0500
Subject: [PATCH 32/48] Add test for TACOT decomposition; fix mistakes in the
 model

---
 mirgecom/materials/tacot.py | 31 +++++++++-------
 test/test_wallmodel.py      | 72 +++++++++++++++++++++++++++++++++++++
 2 files changed, 91 insertions(+), 12 deletions(-)
 create mode 100644 test/test_wallmodel.py

diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index 52063ba74..a9c31b068 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -54,30 +54,37 @@ class Pyrolysis:
     a minimum temperature, are considered based on the resin constituents.
     The third reaction is the fiber oxidation, which is not handled here for now.
 
+    .. automethod:: __init__
+    .. automethod:: get_decomposition_parameters
     .. automethod:: get_source_terms
     """
 
-    def __init__(
-            self, *, virgin_mass=120.0, char_mass=60.0, fiber_mass=160.0,
-            pre_exponential=(12000.0, 4.48e9)):
+    def __init__(self, *,
+                 virgin_mass=120.0, char_mass=60.0, fiber_mass=160.0,
+                 pre_exponential=(12000.0, 4.48e9),
+                 decomposition_temperature=(333.3, 555.6)):
         """Initialize TACOT parameters."""
-        self._Tcrit = np.array([333.3, 555.6])
         self._virgin_mass = virgin_mass
         self._char_mass = char_mass
         self._fiber_mass = fiber_mass
 
+        if len(decomposition_temperature) != 2:
+            raise ValueError("TACOT model requires 2 starting temperatures.")
+        self._Tcrit = np.array(decomposition_temperature)
+
         if len(pre_exponential) != 2:
-            raise ValueError("TACOT degradation model requires 2 pre-exponentials.")
+            raise ValueError("TACOT model requires 2 pre-exponentials.")
         self._pre_exp = np.array(pre_exponential)
 
-    def get_tacot_parameters(self):
+    def get_decomposition_parameters(self):
         """Return the parameters of TACOT decomposition for inspection."""
-        return {"virgin mass": self._virgin_mass,
+        return {"virgin_mass": self._virgin_mass,
                 "char_mass": self._char_mass,
                 "fiber_mass": self._fiber_mass,
-                "reaction_weights": np.array([self._virgin_mass*0.75,
-                                              self._virgin_mass*0.25]),
-                "pre_exponential": self._pre_exp}
+                "reaction_weights": np.array([self._virgin_mass*0.25,
+                                              self._virgin_mass*0.75]),
+                "pre_exponential": self._pre_exp,
+                "temperature": self._Tcrit}
 
     def get_source_terms(self, temperature, chi):
         r"""Return the source terms of pyrolysis decomposition for TACOT.
@@ -101,8 +108,8 @@ def get_source_terms(self, temperature, chi):
 
         # The density parameters are hard-coded for TACOT, depending on
         # virgin and char volume fractions.
-        w1 = self._virgin_mass*0.25*(chi[0]/self._virgin_mass*0.25)**3
-        w2 = self._virgin_mass*0.75*(chi[1]/self._virgin_mass*0.75 - 2./3.)**3
+        w1 = self._virgin_mass*0.25*(chi[0]/(self._virgin_mass*0.25))**3
+        w2 = self._virgin_mass*0.75*(chi[1]/(self._virgin_mass*0.75) - 2./3.)**3
 
         return make_obj_array([
             # reaction 1
diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
new file mode 100644
index 000000000..c4dcea530
--- /dev/null
+++ b/test/test_wallmodel.py
@@ -0,0 +1,72 @@
+"""Test wall-model related functions."""
+
+__copyright__ = """Copyright (C) 2023 University of Illinois Board of Trustees"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import numpy as np
+from pytools.obj_array import make_obj_array
+from meshmode.array_context import (  # noqa
+    pytest_generate_tests_for_pyopencl_array_context
+    as pytest_generate_tests)
+
+
+def test_tacot_decomposition():
+    """Check the wall degradation model."""
+    temperature = 900.0
+
+    from mirgecom.materials.tacot import Pyrolysis
+    decomposition = Pyrolysis()
+    chi = np.array([30.0, 90.0, 160.0])
+
+    tol = 1e-8
+
+    tacot_decomp = decomposition.get_decomposition_parameters()
+
+    print(tacot_decomp)
+
+    # virgin_mass = tacot_decomp["virgin_mass"]
+    # char_mass = tacot_decomp["char_mass"]
+    # fiber_mass = tacot_decomp["fiber_mass"]
+    weights = tacot_decomp["reaction_weights"]
+    pre_exp = tacot_decomp["pre_exponential"]
+    Tcrit = tacot_decomp["temperature"]  # noqa N806
+
+    # The density parameters are hard-coded for TACOT, depending on
+    # virgin and char volume fractions.
+    w1 = weights[0]*(chi[0]/(weights[0]))**3
+    w2 = weights[1]*(chi[1]/(weights[1]) - 2./3.)**3
+
+    solid_mass_rhs = make_obj_array([
+        # reaction 1
+        np.where(np.less(temperature, Tcrit[0]),
+            0.0, (-w1 * pre_exp[0] * np.exp(-8556.000/temperature))),
+        # reaction 2
+        np.where(np.less(temperature, Tcrit[1]),
+            0.0, (-w2 * pre_exp[1] * np.exp(-20444.44/temperature))),
+        # fiber oxidation: include in the RHS but don't do anything with it.
+        np.zeros_like(temperature)])
+
+    sample_source_gas = -sum(solid_mass_rhs)
+
+    assert solid_mass_rhs[0] + 26.7676118539965 < tol
+    assert solid_mass_rhs[1] + 2.03565420370596 < tol
+    assert sample_source_gas - 28.8032660577024 < tol

From 53b9e1bb4da956b5482e23c13ae49b7d240f8c19 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 27 Mar 2024 11:53:26 -0500
Subject: [PATCH 33/48] Add some extra missing mypy

---
 mirgecom/materials/carbon_fiber.py               | 12 +++++++-----
 mirgecom/materials/prescribed_porous_material.py |  3 +--
 mirgecom/materials/tacot.py                      |  7 ++-----
 3 files changed, 10 insertions(+), 12 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index c4e72636b..d1e984cde 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -72,7 +72,7 @@ def _get_wall_effective_surface_area_fiber(self, tau: DOFArray) -> DOFArray:
         return self.puma_effective_surface_area(tau)
 
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+                         rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
         """Return the effective source terms for fiber oxidation.
 
         Parameters
@@ -154,14 +154,16 @@ def void_fraction(self, tau: DOFArray) -> DOFArray:
         """
         return 1.0 - self.volume_fraction(tau)
 
-    def enthalpy(self, temperature, tau=None) -> DOFArray:
+    def enthalpy(self, temperature: DOFArray,
+                 tau: Optional[DOFArray] = None) -> DOFArray:
         r"""Evaluate the solid enthalpy $h_s$ of the fibers."""
         return (
             - 3.37112113e-11*temperature**5 + 3.13156695e-07*temperature**4
             - 1.17026962e-03*temperature**3 + 2.29194901e+00*temperature**2
             - 3.62422269e+02*temperature**1 - 5.96993843e+04)
 
-    def heat_capacity(self, temperature, tau=None) -> DOFArray:
+    def heat_capacity(self, temperature: DOFArray,
+                      tau: Optional[DOFArray] = None) -> DOFArray:
         r"""Evaluate the heat capacity $C_{p_s}$ of the fibers.
 
         The coefficients are obtained with the analytical derivative of the
@@ -173,7 +175,7 @@ def heat_capacity(self, temperature, tau=None) -> DOFArray:
             - 3.62422269e+02)
 
     # ~~~~~~~~ fiber conductivity
-    def thermal_conductivity(self, temperature, tau=None) -> np.ndarray:
+    def thermal_conductivity(self, temperature, tau) -> np.ndarray:
         r"""Evaluate the thermal conductivity $\kappa$ of the fibers.
 
         It accounts for anisotropy and oxidation progress.
@@ -212,7 +214,7 @@ def permeability(self, tau: DOFArray) -> np.ndarray:
         return permeability
 
     def emissivity(self, temperature: DOFArray,  # type: ignore[override]
-            tau: Optional[DOFArray] = None) -> DOFArray:
+                   tau: Optional[DOFArray] = None) -> DOFArray:
         """Emissivity for energy radiation."""
         return (
             + 2.26413679e-18*temperature**5 - 2.03008004e-14*temperature**4
diff --git a/mirgecom/materials/prescribed_porous_material.py b/mirgecom/materials/prescribed_porous_material.py
index f69673637..9b7ced291 100644
--- a/mirgecom/materials/prescribed_porous_material.py
+++ b/mirgecom/materials/prescribed_porous_material.py
@@ -79,8 +79,7 @@ def heat_capacity(self, temperature: DOFArray,
         r"""Evaluate the heat capacity $C_{p_s}$."""
         return self._heat_capacity_func(temperature)
 
-    def thermal_conductivity(self, temperature: DOFArray,
-                                   tau: DOFArray) -> DOFArray:
+    def thermal_conductivity(self, temperature: DOFArray, tau: DOFArray) -> DOFArray:
         r"""Evaluate the thermal conductivity $\kappa$"""
         return self._thermal_conductivity_func(temperature)
 
diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index a9c31b068..b3baa6422 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -148,7 +148,6 @@ def __init__(self, char_mass, virgin_mass):
 
         Parameters
         ----------
-
         virgin_mass: float
             initial mass of the material. The fiber and all resin components
             must be considered.
@@ -181,8 +180,7 @@ def enthalpy(self, temperature: DOFArray, tau: DOFArray) -> DOFArray:
 
         return virgin*tau + char*(1.0 - tau)
 
-    def heat_capacity(self, temperature: DOFArray,
-                      tau: DOFArray) -> DOFArray:
+    def heat_capacity(self, temperature: DOFArray, tau: DOFArray) -> DOFArray:
         r"""Solid heat capacity $C_{p_s}$ as a function of pyrolysis progress."""
         actx = temperature.array_context
 
@@ -199,8 +197,7 @@ def heat_capacity(self, temperature: DOFArray,
 
         return virgin*tau + char*(1.0 - tau)
 
-    def thermal_conductivity(self, temperature: DOFArray,
-                             tau: DOFArray) -> DOFArray:
+    def thermal_conductivity(self, temperature: DOFArray, tau: DOFArray) -> DOFArray:
         """Solid thermal conductivity as a function of pyrolysis progress."""
         virgin = (
             + 2.31290019732353e-17*temperature**5 - 2.167785032562e-13*temperature**4

From 963cc96b0b1a278e5d176cbd542daf4cc848f72f Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Thu, 28 Mar 2024 10:01:10 -0500
Subject: [PATCH 34/48] Some fixes/updates to work with porous slab

---
 mirgecom/boundary.py              | 22 ++++++++++++++++++++--
 mirgecom/materials/initializer.py | 12 ++++++------
 mirgecom/materials/tacot.py       |  8 ++++++--
 mirgecom/wall_model.py            |  1 +
 4 files changed, 33 insertions(+), 10 deletions(-)

diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py
index 4206bbe07..26ea591e4 100644
--- a/mirgecom/boundary.py
+++ b/mirgecom/boundary.py
@@ -1297,6 +1297,8 @@ class PressureOutflowBoundary(MengaldoBoundaryCondition):
     def __init__(self, boundary_pressure=101325):
         """Initialize the boundary condition object."""
         self._pressure = boundary_pressure
+        self._slip = _SlipBoundaryComponent()
+        self._adiabatic = _AdiabaticBoundaryComponent()
 
     def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Get the exterior solution on the boundary.
@@ -1414,11 +1416,27 @@ def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs):
     def grad_cv_bc(self, dcoll, dd_bdry, gas_model, state_minus, grad_cv_minus,
                    normal, **kwargs):
         """Return grad(CV) to be used in the boundary calculation of viscous flux."""
-        return grad_cv_minus
+#        return grad_cv_minus
+
+        dd_bdry = as_dofdesc(dd_bdry)
+        normal = state_minus.array_context.thaw(dcoll.normal(dd_bdry))
+        state_bc = self.state_bc(
+            dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model,
+            state_minus=state_minus, **kwargs)
+
+        grad_v_bc = self._slip.grad_velocity_bc(
+            state_minus, state_bc, grad_cv_minus, normal)
+
+        grad_mom_bc = (
+            state_bc.mass_density * grad_v_bc
+            + np.outer(state_bc.velocity, grad_cv_minus.mass))
+
+        return grad_cv_minus.replace(momentum=grad_mom_bc)
 
     def grad_temperature_bc(self, dcoll, dd_bdry, grad_t_minus, normal, **kwargs):
         """Return grad(temperature) to be used in viscous flux at wall."""
-        return grad_t_minus
+#        return grad_t_minus
+        return self._adiabatic.grad_temperature_bc(grad_t_minus, normal)
 
 
 class RiemannInflowBoundary(MengaldoBoundaryCondition):
diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 79c98574b..ff89f0d83 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -75,8 +75,9 @@ def __call__(self, x_vec, wall_model):
 class PorousWallInitializer:
     """State initializer for porous materials in the unified-domain solver."""
 
-    def __init__(self, *, temperature, material_densities, species=None,
-                 velocity=None, pressure=None, density=None, porous_region=None):
+    def __init__(
+            self, *, temperature, material_densities, species_mass_fractions=None,
+            velocity=None, pressure=None, density=None, porous_region=None):
         """Initialize the object for porous materials.
 
         Parameters
@@ -93,10 +94,9 @@ def __init__(self, *, temperature, material_densities, species=None,
         self._wall_density = material_densities
         self._porous_region = porous_region
 
-        if species is not None:
-            self._y = species
+        if species_mass_fractions is not None:
+            self._y = species_mass_fractions
         else:
-            import numpy as np
             self._y = np.empty((0,), dtype=object)
 
     def __call__(self, x_vec, gas_model):
@@ -152,7 +152,7 @@ def __call__(self, x_vec, gas_model):
 
         # start with zero velocity inside the material
         if self._velocity is not None:
-            momentum = (1.0-self._porous_region) * (eps_gas*self._velocity)
+            momentum = (1.0-self._porous_region) * (eps_rho_gas*self._velocity)
         else:
             momentum = make_obj_array([zeros for _ in range(dim)])
 
diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index 01586a253..25f027c51 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -120,8 +120,12 @@ def get_source_terms(self, temperature, chi):
             actx.np.where(actx.np.less(temperature, self._Tcrit[1]),
                 0.0, (-w2 * self._pre_exp[1] * actx.np.exp(-20444.44/temperature))
                 ),
-            # fiber oxidation: include in the RHS but don't do anything with it.
-            actx.np.zeros_like(temperature)])
+            # TODO: fiber oxidation:
+            # Include a resin-dependent term (probably linear) to slow down
+            # oxidation when there is a lot of resin, but make it goes with
+            # the "right" rate when the resin is consumed.
+            actx.np.zeros_like(temperature)
+            ])
 
 
 class TacotEOS(PorousWallEOS):
diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index bc2e68e6d..ffaddb2b8 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -445,6 +445,7 @@ def decomposition_progress(self, material_densities) -> DOFArray:
         Where $\tau=1$, the material is locally virgin. On the other hand, if
         $\tau=0$, then the fibers were all consumed.
         """
+        # TODO consider the resin+carbon case, where both are degrading.
         mass = self.solid_density(material_densities)
         return self.wall_eos.decomposition_progress(mass)
 

From 43b52a3231bc8d7bdc4d88b030c4c740adb9224d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Thu, 28 Mar 2024 11:53:51 -0500
Subject: [PATCH 35/48] Fix test

---
 test/test_wallmodel.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index 9cda7ceaf..2f924758b 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -122,7 +122,7 @@ def test_wall_eos(actx_factory, order, my_material, my_eos, use_diffused_interfa
     else:
         porous_region = None
     sample_init = PorousWallInitializer(
-        pressure=pressure, temperature=temperature, species=y,
+        pressure=pressure, temperature=temperature, species_mass_fractions=y,
         material_densities=material_densities, porous_region=porous_region)
 
     cv, solid_densities = sample_init(nodes, gas_model)

From 1a083c1c6ca4cc0245925ae5a41b28ab3473be28 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 2 Apr 2024 14:45:20 -0500
Subject: [PATCH 36/48] Force wall to be positive; add option to return whole
 set of WVs

---
 mirgecom/materials/initializer.py | 24 +++++++++++++++++++-----
 mirgecom/wall_model.py            |  7 +++++--
 2 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index ff89f0d83..6f60a9441 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -32,6 +32,7 @@
 from pytools.obj_array import make_obj_array
 from mirgecom.fluid import make_conserved
 from mirgecom.wall_model import SolidWallConservedVars
+from mirgecom.wall_model import PorousWallVars
 
 
 class SolidWallInitializer:
@@ -99,7 +100,7 @@ def __init__(
         else:
             self._y = np.empty((0,), dtype=object)
 
-    def __call__(self, x_vec, gas_model):
+    def __call__(self, x_vec, gas_model, return_wv=False):
         """Evaluate the wall+gas properties for porous materials.
 
         Parameters
@@ -116,6 +117,8 @@ def __call__(self, x_vec, gas_model):
             both gas and porous material properties.
         wall_density: numpy.ndarray or :class:`meshmode.dof_array.DOFArray`
             The densities of each one of the materials
+        wv: :class:`mirgecom.fluid.wall_model.PorousWallVars`
+            Wall dependent variables
         """
         actx = x_vec[0].array_context
         zeros = actx.np.zeros_like(x_vec[0])
@@ -128,9 +131,18 @@ def __call__(self, x_vec, gas_model):
             porous_region = ones
 
         # wall-only properties
-        wall_density = self._wall_density * porous_region
-        tau = gas_model.decomposition_progress(wall_density)
-        eps_rho_solid = gas_model.solid_density(wall_density)
+        material_densities = self._wall_density * porous_region
+        tau = gas_model.decomposition_progress(material_densities)
+        eps_rho_solid = gas_model.solid_density(material_densities)
+
+        wv = PorousWallVars(
+            material_densities=material_densities,
+            tau=tau,
+            density=eps_rho_solid,
+            void_fraction=gas_model.wall_eos.void_fraction(tau=tau),
+            permeability=gas_model.wall_eos.permeability(tau=tau),
+            tortuosity=gas_model.wall_eos.tortuosity(tau=tau)
+        )
 
         # coupled properties
         temperature = self._temp + zeros
@@ -168,4 +180,6 @@ def __call__(self, x_vec, gas_model):
         cv = make_conserved(dim=dim, mass=eps_rho_gas, energy=bulk_energy,
                             momentum=momentum, species_mass=species_mass)
 
-        return cv, wall_density
+        if return_wv:
+            return cv, wv
+        return cv, wv.material_densities
diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index ffaddb2b8..60669808d 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -435,9 +435,12 @@ def solid_density(self, material_densities) -> DOFArray:
         .. math::
             \epsilon_s \rho_s = \sum_i^N \epsilon_i \rho_i
         """
+        actx = material_densities.array_context
         if isinstance(material_densities, DOFArray):
-            return material_densities
-        return sum(material_densities)
+            return actx.np.absolute(material_densities)
+        else:
+            1/0
+        # XXX return sum(material_densities)
 
     def decomposition_progress(self, material_densities) -> DOFArray:
         r"""Evaluate the progress ratio $\tau$ of the decomposition.

From 4b4b80b0d7dd3c5b0ceee7cdf3f90a58b609ba53 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 2 Apr 2024 17:53:16 -0500
Subject: [PATCH 37/48] Actually, dont do anything to the wall density but only
 modify the source term. Still pending add a test to the oxidation model...

---
 mirgecom/materials/carbon_fiber.py | 12 ++++++++----
 mirgecom/wall_model.py             | 13 ++++++++-----
 test/test_wallmodel.py             |  5 +++++
 3 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 9793bf391..aecd7790e 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -105,8 +105,9 @@ class Y3_Oxidation_Model(Oxidation):  # noqa N801
     .. automethod:: get_source_terms
     """
 
-    def __init__(self, wall_material):
+    def __init__(self, wall_material, initial_fiber_radius=5e-6):
         self._material = wall_material
+        self._init_fiber_radius = initial_fiber_radius
 
     def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
         r"""Evaluate the effective surface of the fibers.
@@ -117,13 +118,16 @@ def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
             \tau = \frac{m}{m_0} = \frac{\pi r^2/L}{\pi r_0^2/L} = \frac{r^2}{r_0^2}
         """
         actx = tau.array_context
+        zeros = actx.np.zeros_like(tau)
 
-        original_fiber_radius = 5e-6  # half the diameter
-        fiber_radius = original_fiber_radius*actx.np.sqrt(tau)
+        # half the diameter
+        fiber_radius = \
+            self._init_fiber_radius * actx.np.sqrt(actx.np.minimum(tau, zeros))
 
         epsilon_0 = self._material.volume_fraction(tau=1.0)
 
-        return 2.0*epsilon_0/original_fiber_radius**2*fiber_radius
+        # TODO simplify this...
+        return 2.0 * epsilon_0 * fiber_radius/self._init_fiber_radius**2
 
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
             rhoY_o2: DOFArray):  # noqa N803
diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 60669808d..7f9de1a6c 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -435,12 +435,15 @@ def solid_density(self, material_densities) -> DOFArray:
         .. math::
             \epsilon_s \rho_s = \sum_i^N \epsilon_i \rho_i
         """
-        actx = material_densities.array_context
         if isinstance(material_densities, DOFArray):
-            return actx.np.absolute(material_densities)
-        else:
-            1/0
-        # XXX return sum(material_densities)
+            return material_densities
+        return sum(material_densities)
+        # if isinstance(material_densities, DOFArray):
+        #     actx = material_densities.array_context
+        #     return actx.np.absolute(material_densities)
+        # else:
+        #     actx = material_densities[0].array_context
+        #     return sum(actx.np.absolute(material_densities))
 
     def decomposition_progress(self, material_densities) -> DOFArray:
         r"""Evaluate the progress ratio $\tau$ of the decomposition.
diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index 2f924758b..be12c9ce9 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -193,3 +193,8 @@ def test_tacot_decomposition():
     assert solid_mass_rhs[0] + 26.7676118539965 < tol
     assert solid_mass_rhs[1] + 2.03565420370596 < tol
     assert sample_source_gas - 28.8032660577024 < tol
+
+
+# TODO
+def test_oxidation_model():
+    assert 1 == 1

From 108d8fde6140edf1c082a77e40e957d59fea9b62 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 23 Apr 2024 16:03:51 -0500
Subject: [PATCH 38/48] temp

---
 mirgecom/materials/carbon_fiber.py | 113 +++++++++++++++--------------
 test/test_wallmodel.py             |   2 -
 2 files changed, 57 insertions(+), 58 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index aecd7790e..60e3956fe 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -40,15 +40,28 @@
 
 
 class Oxidation:
-    """Abstract interface for wall oxidation model.
+    """Abstract interface for wall oxidation models.
 
     .. automethod:: get_source_terms
     """
 
     @abstractmethod
-    def get_source_terms(self, temperature: DOFArray,
-            tau: DOFArray, rhoY_o2: DOFArray) -> Union[DOFArray, tuple]:  # noqa N803
-        r"""Source terms of fiber oxidation."""
+    def get_source_terms(
+            self, temperature: DOFArray, tau: DOFArray, rhoY_o2: DOFArray,  # noqa N803
+            rhoY_o: Optional[DOFArray] = None) -> Union[DOFArray, tuple]:
+        r"""Source terms of fiber oxidation.
+
+        Parameters
+        ----------
+        temperature:
+            gas and solid temperature assuming local-equilibrium
+        tau:
+            the progress ratio of the oxidation: 1 for virgin, 0 for fully oxidized
+        rhoY_o2:
+            the mass fraction of molecular oxygen
+        rhoY_o:
+            the mass fraction of atomic oxygen
+        """
         raise NotImplementedError()
 
 
@@ -67,17 +80,8 @@ def _get_wall_effective_surface_area_fiber(self, tau: DOFArray) -> DOFArray:
         return 1.1794e5*(1.0 - 0.0547736137*progress - 0.9336950992*progress**2)
 
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-                         rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
-        """Return the effective source terms for fiber oxidation.
-
-        Parameters
-        ----------
-        temperature:
-        tau:
-            the progress ratio of the oxidation: 1 for virgin, 0 for fully oxidized
-        rhoY_o2:
-            the mass fraction of oxygen
-        """
+                         rhoY_o2: DOFArray, rhoY_o=None) -> DOFArray:  # noqa N803
+        """Return the effective source terms for fiber oxidation."""
         actx = temperature.array_context
 
         mw_o = 15.999
@@ -102,6 +106,25 @@ class Y3_Oxidation_Model(Oxidation):  # noqa N801
     .. math::
         C_{(s)} + O_2 \to CO_2
 
+    The degradation is computed as
+
+    .. math::
+        \frac{\partial \rho_s}{\partial t} =
+        \rho_C \frac{\partial \epsilon}{\partial t} = \dot{\omega}_\text{het}
+
+    The fiber radius as a function of mass loss $\tau$ is given by
+
+    .. math::
+        \tau = \frac{m}{m_0} = \frac{\epsilon}{\epsilon_0} =
+        \frac{\pi r^2/L}{\pi r_0^2/L} = \frac{r^2}{r_0^2}
+
+    and leads to the wall recession as
+
+    .. math::
+        \frac{\partial \epsilon}{\partial t} =
+        \frac{\partial \epsilon}{\partial r} \frac{\partial r}{\partial t} =
+        2 \frac{\epsilon_0}{r_0^2} r = 2 \frac{\epsilon_0}{r_0} \sqrt{\tau}
+
     .. automethod:: get_source_terms
     """
 
@@ -110,37 +133,18 @@ def __init__(self, wall_material, initial_fiber_radius=5e-6):
         self._init_fiber_radius = initial_fiber_radius
 
     def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
-        r"""Evaluate the effective surface of the fibers.
-
-        The fiber radius as a function of mass loss $\tau$ is given by
-
-        .. math::
-            \tau = \frac{m}{m_0} = \frac{\pi r^2/L}{\pi r_0^2/L} = \frac{r^2}{r_0^2}
-        """
+        r"""Evaluate the effective surface of the fibers."""
         actx = tau.array_context
         zeros = actx.np.zeros_like(tau)
-
-        # half the diameter
-        fiber_radius = \
-            self._init_fiber_radius * actx.np.sqrt(actx.np.minimum(tau, zeros))
-
-        epsilon_0 = self._material.volume_fraction(tau=1.0)
-
-        # TODO simplify this...
-        return 2.0 * epsilon_0 * fiber_radius/self._init_fiber_radius**2
+        return (
+            2.0 * self._material.volume_fraction(tau=1.0)/self._init_fiber_radius
+            * actx.np.sqrt(actx.np.maximum(tau, zeros))  # only use positive tau
+        )
 
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: DOFArray):  # noqa N803
+            rhoY_o2: DOFArray, rhoY_o=None):  # noqa N803
         r"""Return the effective source terms for the oxidation.
 
-        Parameters
-        ----------
-        temperature:
-        tau:
-            the progress ratio of the oxidation
-        rhoY_o2:
-            the mass fraction of oxygen
-
         Returns
         -------
         sources: tuple
@@ -148,6 +152,9 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
         """
         actx = temperature.array_context
 
+        # ignore small but negative fractions of O2
+        rhoY_o2 = actx.np.maximum(rhoY_o2, actx.np.zeros_like(rhoY_o2))
+
         mw_c = 12.011
         mw_o = 15.999
         mw_o2 = mw_o*2
@@ -168,7 +175,7 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
 class OxidationModel(Oxidation):
     """Evaluate the source terms for the carbon fiber oxidation.
 
-    The user must specify in the driver the functions for oxidation.
+    The user must specify the functions for oxidation in the driver.
     (Tentatively) Generalizing this class makes it easier for adding new,
     more complex models and for UQ runs.
 
@@ -176,7 +183,7 @@ class OxidationModel(Oxidation):
     """
 
     def __init__(self, surface_area_func, oxidation_func):
-        """Initialize the general oxidation class.
+        """Initialize the user-specified oxidation class.
 
         Parameters
         ----------
@@ -188,22 +195,16 @@ def __init__(self, surface_area_func, oxidation_func):
         self._surface_func = surface_area_func
         self._oxidation_func = oxidation_func
 
-    # TODO we potentially have to include atomic oxygen as well
-    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
-        """Return the effective source terms for the oxidation.
-
-        Parameters
-        ----------
-        temperature:
-        tau:
-            the progress ratio of the oxidation
-        rhoY_o2:
-            the mass fraction of oxygen
-        """
+    def get_source_terms(
+            self, temperature: DOFArray, tau: DOFArray, rhoY_o2: DOFArray,  # noqa N803
+            rhoY_o: Optional[DOFArray] = None) -> DOFArray:
+        """Return the effective source terms for the oxidation."""
         area = self._surface_func(tau)
+        if rhoY_o is None:
+            return self._oxidation_func(temperature=temperature, fiber_area=area,
+                                        rhoY_o2=rhoY_o2)
         return self._oxidation_func(temperature=temperature, fiber_area=area,
-                                    rhoY_o2=rhoY_o2)
+                                    rhoY_o2=rhoY_o2, rhoY_o=rhoY_o)
 
 
 class FiberEOS(PorousWallEOS):
diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index be12c9ce9..07b7e9409 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -164,8 +164,6 @@ def test_tacot_decomposition():
 
     tacot_decomp = decomposition.get_decomposition_parameters()
 
-    print(tacot_decomp)
-
     # virgin_mass = tacot_decomp["virgin_mass"]
     # char_mass = tacot_decomp["char_mass"]
     # fiber_mass = tacot_decomp["fiber_mass"]

From 8917682c8f62e75b52e32faa9080b6d52ec72692 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 1 May 2024 11:20:45 -0500
Subject: [PATCH 39/48] Fix to fluxes and BCs and related stuff; Add oxidation
 models; Modify initializer to setup velocity

---
 mirgecom/boundary.py               | 24 +++--------------
 mirgecom/eos.py                    | 13 +++++++++
 mirgecom/gas_model.py              | 20 +++++++++-----
 mirgecom/inviscid.py               | 42 +++++++++++++++++++++++++++---
 mirgecom/materials/carbon_fiber.py | 14 ++++++----
 mirgecom/materials/initializer.py  |  6 +++--
 6 files changed, 81 insertions(+), 38 deletions(-)

diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py
index 26ea591e4..8c7fef056 100644
--- a/mirgecom/boundary.py
+++ b/mirgecom/boundary.py
@@ -1027,7 +1027,7 @@ def inviscid_divergence_flux(self, dcoll, dd_bdry, gas_model, state_minus,
                                  **kwargs):
         """Get the inviscid boundary flux for the divergence operator."""
         normal = state_minus.array_context.thaw(dcoll.normal(dd_bdry))
-        return inviscid_flux(state_minus)@normal
+        return inviscid_flux(state_minus, gas_model=gas_model)@normal
 
     def cv_gradient_flux(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Get the cv flux for *dd_bdry* for use in the gradient operator."""
@@ -1297,8 +1297,6 @@ class PressureOutflowBoundary(MengaldoBoundaryCondition):
     def __init__(self, boundary_pressure=101325):
         """Initialize the boundary condition object."""
         self._pressure = boundary_pressure
-        self._slip = _SlipBoundaryComponent()
-        self._adiabatic = _AdiabaticBoundaryComponent()
 
     def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
         """Get the exterior solution on the boundary.
@@ -1416,27 +1414,11 @@ def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs):
     def grad_cv_bc(self, dcoll, dd_bdry, gas_model, state_minus, grad_cv_minus,
                    normal, **kwargs):
         """Return grad(CV) to be used in the boundary calculation of viscous flux."""
-#        return grad_cv_minus
-
-        dd_bdry = as_dofdesc(dd_bdry)
-        normal = state_minus.array_context.thaw(dcoll.normal(dd_bdry))
-        state_bc = self.state_bc(
-            dcoll=dcoll, dd_bdry=dd_bdry, gas_model=gas_model,
-            state_minus=state_minus, **kwargs)
-
-        grad_v_bc = self._slip.grad_velocity_bc(
-            state_minus, state_bc, grad_cv_minus, normal)
-
-        grad_mom_bc = (
-            state_bc.mass_density * grad_v_bc
-            + np.outer(state_bc.velocity, grad_cv_minus.mass))
-
-        return grad_cv_minus.replace(momentum=grad_mom_bc)
+        return grad_cv_minus
 
     def grad_temperature_bc(self, dcoll, dd_bdry, grad_t_minus, normal, **kwargs):
         """Return grad(temperature) to be used in viscous flux at wall."""
-#        return grad_t_minus
-        return self._adiabatic.grad_temperature_bc(grad_t_minus, normal)
+        return grad_t_minus
 
 
 class RiemannInflowBoundary(MengaldoBoundaryCondition):
diff --git a/mirgecom/eos.py b/mirgecom/eos.py
index 9ad289767..d839d467f 100644
--- a/mirgecom/eos.py
+++ b/mirgecom/eos.py
@@ -833,3 +833,16 @@ def get_species_source_terms(self, cv: ConservedVars, temperature: DOFArray):
 
         return make_conserved(dim, rho_source, energy_source, mom_source,
                               species_sources)
+
+    def get_mole_fractions(self, species_mass_fractions):
+        """Get mole fractions from mass fractions."""
+        mix_mol_weight = \
+            self._pyrometheus_mech.get_mix_molecular_weight(species_mass_fractions)
+        return self._pyrometheus_mech.get_mole_fractions(mix_mol_weight,
+                                                         species_mass_fractions)
+
+    def get_mass_fractions(self, species_mole_fractions):
+        """Get mass fractions from mole fractions."""
+        mol_weights = self.get_species_molecular_weights()
+        mmw = sum(species_mole_fractions*mol_weights)
+        return species_mole_fractions*mol_weights/mmw
diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py
index 1e2ce46a5..1b95d4516 100644
--- a/mirgecom/gas_model.py
+++ b/mirgecom/gas_model.py
@@ -436,14 +436,22 @@ def make_fluid_state(cv, gas_model,
             tortuosity=gas_model.wall_eos.tortuosity(tau=tau)
         )
 
-        temperature = gas_model.get_temperature(cv=cv, wv=wv,
-            tseed=temperature_seed)
-
-        pressure = gas_model.get_pressure(cv, wv, temperature)
+        pressure = None
+        temperature = None
 
         if limiter_func:
-            cv = limiter_func(cv=cv, wv=wv, pressure=pressure,
-                              temperature=temperature, dd=limiter_dd)
+            rv = limiter_func(cv=cv, wv=wv, temperature_seed=temperature_seed,
+                              gas_model=gas_model, dd=limiter_dd)
+            if isinstance(rv, np.ndarray):
+                cv, pressure, temperature = rv
+            else:
+                cv = rv
+
+        if temperature is None:
+            temperature = gas_model.get_temperature(cv=cv, wv=wv,
+                                                    tseed=temperature_seed)
+        if pressure is None:
+            pressure = gas_model.get_pressure(cv, wv, temperature)
 
         from mirgecom.eos import MixtureEOS
         if isinstance(gas_model.eos, MixtureEOS):
diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py
index aa5d8420b..f5103a9ff 100644
--- a/mirgecom/inviscid.py
+++ b/mirgecom/inviscid.py
@@ -65,9 +65,18 @@
 def inviscid_flux(state, gas_model=None):
     r"""Compute the inviscid flux vectors from fluid conserved vars *cv*.
 
-    The inviscid fluxes are
-    $(\rho\vec{V},(\rho{E}+p)\vec{V},\rho(\vec{V}\otimes\vec{V})
-    +p\mathbf{I}, \rho{Y_s}\vec{V})$
+    The inviscid fluxes for pure-fluids and porous flows are, respectively,
+    $$
+    \begin{bmatrix}(\rho{v})_{j}\\
+    \left(\rho{E}+p\right){v}_{j}\\
+    \left((\rho{v})_{j}{v}_{i}+p\delta_{ij}\right)\\
+    (\rho{Y})_{\alpha}{v}_{j}\end{bmatrix}
+    \text{ and }
+    \begin{bmatrix}(\epsilon \rho{v})_{j}\\
+    \left(\epsilon \rho{E}+\epsilon p\right){v}_{j}\\
+    \left((\epsilon \rho{v})_{j}{v}_{i}+p\delta_{ij}\right)\\
+    (\epsilon \rho{Y})_{\alpha}{v}_{j}\end{bmatrix} \text{ .}
+    $$
 
     .. note::
 
@@ -166,6 +175,31 @@ def inviscid_facial_flux_rusanov(state_pair, gas_model, normal):
     actx = state_pair.int.array_context
     lam = actx.np.maximum(state_pair.int.wavespeed, state_pair.ext.wavespeed)
     from mirgecom.flux import num_flux_lfr
+
+    # Use only the gas energy in penalization terms
+    # TODO think of a better way to generalize this for different fluxes??
+    if isinstance(state_pair.int, PorousFlowFluidState):
+        gas_energy_int = state_pair.int.mass_density*(
+            gas_model.eos.get_internal_energy(
+                temperature=state_pair.int.temperature,
+                species_mass_fractions=state_pair.int.cv.species_mass_fractions)
+            + 0.5*np.dot(state_pair.int.velocity, state_pair.int.velocity)
+        )
+        q_minus = state_pair.int.cv.replace(energy=gas_energy_int)
+
+        gas_energy_ext = state_pair.ext.mass_density*(
+            gas_model.eos.get_internal_energy(
+                temperature=state_pair.ext.temperature,
+                species_mass_fractions=state_pair.ext.cv.species_mass_fractions)
+            + 0.5*np.dot(state_pair.ext.velocity, state_pair.ext.velocity)
+        )
+        q_plus = state_pair.ext.cv.replace(energy=gas_energy_ext)
+
+        return num_flux_lfr(
+            f_minus_normal=inviscid_flux(state_pair.int, gas_model=gas_model)@normal,
+            f_plus_normal=inviscid_flux(state_pair.ext, gas_model=gas_model)@normal,
+            q_minus=q_minus, q_plus=q_plus, lam=lam)
+
     return num_flux_lfr(
         f_minus_normal=inviscid_flux(state_pair.int, gas_model=gas_model)@normal,
         f_plus_normal=inviscid_flux(state_pair.ext, gas_model=gas_model)@normal,
@@ -176,7 +210,7 @@ def inviscid_facial_flux_rusanov(state_pair, gas_model, normal):
 def inviscid_facial_flux_hll(state_pair, gas_model, normal):
     r"""High-level interface for inviscid facial flux using HLL numerical flux.
 
-    The Harten, Lax, van Leer approximate Riemann numerical flux is calculated as:
+    The Harten, Lax, van Leer approximate riemann numerical flux is calculated as:
 
     .. math::
 
diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 60e3956fe..8e801bcbd 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -34,9 +34,9 @@
 from typing import Optional, Union
 from abc import abstractmethod
 import numpy as np
+from pytools.obj_array import make_obj_array
 from meshmode.dof_array import DOFArray
 from mirgecom.wall_model import PorousWallEOS
-from pytools.obj_array import make_obj_array
 
 
 class Oxidation:
@@ -125,12 +125,16 @@ class Y3_Oxidation_Model(Oxidation):  # noqa N801
         \frac{\partial \epsilon}{\partial r} \frac{\partial r}{\partial t} =
         2 \frac{\epsilon_0}{r_0^2} r = 2 \frac{\epsilon_0}{r_0} \sqrt{\tau}
 
+    .. automethod:: __init__
     .. automethod:: get_source_terms
     """
 
-    def __init__(self, wall_material, initial_fiber_radius=5e-6):
+    def __init__(self, wall_material, initial_fiber_radius=5.5e-6, arrhenius=1e5,
+                 activation_energy=-120000.0):
         self._material = wall_material
         self._init_fiber_radius = initial_fiber_radius
+        self._arrhenius = arrhenius
+        self._Ea = activation_energy
 
     def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
         r"""Evaluate the effective surface of the fibers."""
@@ -152,8 +156,8 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
         """
         actx = temperature.array_context
 
-        # ignore small but negative fractions of O2
-        rhoY_o2 = actx.np.maximum(rhoY_o2, actx.np.zeros_like(rhoY_o2))
+        # ignore negative amounts of O2
+        rhoY_o2 = actx.np.maximum(rhoY_o2, actx.np.zeros_like(rhoY_o2))  # noqa N806
 
         mw_c = 12.011
         mw_o = 15.999
@@ -163,7 +167,7 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
 
         eff_surf_area = self._get_wall_effective_surface_area_fiber(tau)
 
-        k_f = 1.0e5*actx.np.exp(-120000.0/(univ_gas_const*temperature))
+        k_f = self._arrhenius*actx.np.exp(self._Ea/(univ_gas_const*temperature))
 
         m_dot_c = - rhoY_o2/mw_o2 * mw_c * eff_surf_area * k_f
         m_dot_o2 = - rhoY_o2/mw_o2 * mw_o2 * eff_surf_area * k_f
diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 6f60a9441..8d601ba55 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -162,9 +162,11 @@ def __call__(self, x_vec, gas_model, return_wv=False):
         else:
             eps_rho_gas = eps_gas*self._mass
 
-        # start with zero velocity inside the material
+        # TODO gotta define better the velocity..
         if self._velocity is not None:
-            momentum = (1.0-self._porous_region) * (eps_rho_gas*self._velocity)
+            velocity = self._velocity/eps_gas
+            momentum = make_obj_array([eps_rho_gas*velocity[i] for i in range(dim)])
+            # momentum = (1.0-self._porous_region) * (eps_rho_gas*self._velocity)
         else:
             momentum = make_obj_array([zeros for _ in range(dim)])
 

From dc2ece6b793b7a78d350125c596eb478b62e4b54 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 13 May 2024 08:58:37 -0500
Subject: [PATCH 40/48] Fix float tau; Make alpha explicit in transport

---
 mirgecom/materials/carbon_fiber.py | 8 ++++++--
 mirgecom/transport.py              | 5 ++---
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index 8e801bcbd..b15610db3 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -140,8 +140,9 @@ def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
         r"""Evaluate the effective surface of the fibers."""
         actx = tau.array_context
         zeros = actx.np.zeros_like(tau)
+        ones = zeros + 1.0
         return (
-            2.0 * self._material.volume_fraction(tau=1.0)/self._init_fiber_radius
+            2.0 * self._material.volume_fraction(tau=ones)/self._init_fiber_radius
             * actx.np.sqrt(actx.np.maximum(tau, zeros))  # only use positive tau
         )
 
@@ -318,7 +319,10 @@ def thermal_conductivity(self, temperature, tau) -> np.ndarray:
     # ~~~~~~~~ other properties
     def volume_fraction(self, tau: DOFArray) -> DOFArray:
         r"""Fraction $\phi$ occupied by the solid."""
-        return 0.12*tau
+        # XXX check if this keeps tau positive and bounds everything else to physical
+        # values (it doesn't crash, but it is non-ideal to have unphysical values)
+        actx = tau.array_context
+        return actx.np.maximum(0.12*tau, 0.0*tau)
 
     def permeability(self, tau: DOFArray) -> np.ndarray:
         r"""Permeability $K$ of the porous material."""
diff --git a/mirgecom/transport.py b/mirgecom/transport.py
index 8d7a7137b..46581e0a8 100644
--- a/mirgecom/transport.py
+++ b/mirgecom/transport.py
@@ -315,9 +315,8 @@ def species_diffusivity(self, cv: ConservedVars,  # type: ignore[override]
             d_{\alpha} = \frac{\kappa}{\rho \; Le \; C_p}
         """
         if self._lewis is not None:
-            return (self._sigma * self.viscosity(cv, dv)/(
-                cv.mass*self._lewis*eos.gamma(cv, dv.temperature))
-            )
+            alpha = self._sigma * self.viscosity(cv, dv)/(cv.mass*eos.gamma(cv, dv.temperature))
+            return self._lewis*alpha
         return self._d_alpha*(0*cv.mass + 1.)
 
 

From fbf0140383fb802cd2aaa3451144a301aeb40cd0 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 13 May 2024 09:15:24 -0500
Subject: [PATCH 41/48] Fix test and flake8

---
 test/test_wallmodel.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index 20bfdd4b7..ad3c94464 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -139,8 +139,19 @@ def test_wall_eos(actx_factory, order, my_material):
         op.norm(dcoll, solid_state.temperature - 900.0, np.inf)) < tol
 
 
-def test_tacot_decomposition():
+def test_tacot_decomposition(actx_factory):
     """Check the wall degradation model."""
+    actx = actx_factory()
+
+    dim = 2
+    nelems = 2
+    order = 2
+    mesh = get_box_mesh(dim, -0.1, 0.1, nelems)
+    dcoll = create_discretization_collection(actx, mesh, order=order)
+
+    nodes = actx.thaw(dcoll.nodes())
+    zeros = actx.np.zeros_like(nodes[0])
+
     temperature = 900.0 + zeros
 
     from mirgecom.materials.tacot import Pyrolysis

From 3cf706d8665e255052862cc8de4ef340ca3b7c9a Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Mon, 13 May 2024 09:21:36 -0500
Subject: [PATCH 42/48] flake8

---
 mirgecom/transport.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/mirgecom/transport.py b/mirgecom/transport.py
index 46581e0a8..8d10b181c 100644
--- a/mirgecom/transport.py
+++ b/mirgecom/transport.py
@@ -315,7 +315,10 @@ def species_diffusivity(self, cv: ConservedVars,  # type: ignore[override]
             d_{\alpha} = \frac{\kappa}{\rho \; Le \; C_p}
         """
         if self._lewis is not None:
-            alpha = self._sigma * self.viscosity(cv, dv)/(cv.mass*eos.gamma(cv, dv.temperature))
+            alpha = (
+                self._sigma * self.viscosity(cv, dv)
+                / (cv.mass * eos.gamma(cv, dv.temperature))
+            )
             return self._lewis*alpha
         return self._d_alpha*(0*cv.mass + 1.)
 

From 16bfdf77b8da796d0e1fc282498e862019c2f283 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 14 May 2024 23:25:12 -0500
Subject: [PATCH 43/48] Fix production conflict

---
 test/test_wallmodel.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index ad3c94464..1a6a8befa 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -22,17 +22,17 @@
 THE SOFTWARE.
 """
 
-import pytest
-import cantera
 import numpy as np
 from pytools.obj_array import make_obj_array
-from grudge import op
-from meshmode.dof_array import DOFArray
 from meshmode.array_context import (  # noqa
     pytest_generate_tests_for_pyopencl_array_context
     as pytest_generate_tests)
-from mirgecom.simutil import get_box_mesh
+import grudge.op as op
 from mirgecom.discretization import create_discretization_collection
+from mirgecom.simutil import get_box_mesh
+import pytest
+import cantera
+from meshmode.dof_array import DOFArray
 from mirgecom.eos import PyrometheusMixture
 from mirgecom.transport import SimpleTransport
 from mirgecom.gas_model import make_fluid_state

From 0868c8b0ed036b5184becc5b6c129d9e27467f89 Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Tue, 14 May 2024 23:57:32 -0500
Subject: [PATCH 44/48] There is no need to have a PR prior to this one..
 Incorporate the specific changes here directly...

---
 mirgecom/materials/initializer.py |  13 +---
 mirgecom/wall_model.py            | 114 +++++++-----------------------
 2 files changed, 29 insertions(+), 98 deletions(-)

diff --git a/mirgecom/materials/initializer.py b/mirgecom/materials/initializer.py
index 8d601ba55..20f021640 100644
--- a/mirgecom/materials/initializer.py
+++ b/mirgecom/materials/initializer.py
@@ -42,9 +42,8 @@ class SolidWallInitializer:
     materials, and/or their combination, subject or not to ablation.
     """
 
-    def __init__(self, temperature, material_densities):
+    def __init__(self, temperature):
         self._temp = temperature
-        self._mass = material_densities
 
     def __call__(self, x_vec, wall_model):
         """Evaluate the wall+gas properties for porous materials.
@@ -62,14 +61,8 @@ def __call__(self, x_vec, wall_model):
             The conserved variables for heat-conduction only materials.
         """
         actx = x_vec[0].array_context
-
-        mass = self._mass + actx.np.zeros_like(x_vec[0])
-        solid_mass = wall_model.solid_density(mass)
-        tau = wall_model.decomposition_progress(mass)
-
-        temperature = self._temp + actx.np.zeros_like(x_vec[0])
-        energy = solid_mass * wall_model.enthalpy(temperature=temperature,
-                                                  tau=tau)
+        mass = wall_model.density() + actx.np.zeros_like(x_vec[0])
+        energy = mass * wall_model.enthalpy(self._temp)
         return SolidWallConservedVars(mass=mass, energy=energy)
 
 
diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 7f9de1a6c..779a6dfba 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -86,11 +86,6 @@
 )
 from mirgecom.transport import GasTransportVars, TransportModel
 
-from grudge.dt_utils import characteristic_lengthscales
-from mirgecom.viscous import get_local_max_species_diffusivity
-import grudge.op as op
-from functools import reduce
-
 
 @with_container_arithmetic(bcast_obj_array=False,
                            bcast_container_types=(DOFArray, np.ndarray),
@@ -116,8 +111,6 @@ def array_context(self):
 class SolidWallDependentVars:
     """Wall dependent variables for heat conduction only materials."""
 
-    tau: DOFArray
-    density: Union[DOFArray, np.ndarray]
     thermal_conductivity: DOFArray
     temperature: DOFArray
 
@@ -139,7 +132,7 @@ class SolidWallModel:
     driver-defined functions that prescribe the actual properties and their
     spatial distribution.
 
-    .. automethod:: solid_density
+    .. automethod:: density
     .. automethod:: heat_capacity
     .. automethod:: enthalpy
     .. automethod:: thermal_diffusivity
@@ -148,72 +141,55 @@ class SolidWallModel:
     .. automethod:: dependent_vars
     """
 
-    def __init__(self, enthalpy_func, heat_capacity_func,
-                 thermal_conductivity_func, decomposition_func=None):
+    def __init__(self, density_func, enthalpy_func, heat_capacity_func,
+                 thermal_conductivity_func):
+        self._density_func = density_func
         self._enthalpy_func = enthalpy_func
         self._heat_capacity_func = heat_capacity_func
         self._thermal_conductivity_func = thermal_conductivity_func
-        self._decomposition_func = decomposition_func
 
-    def solid_density(self, material_densities) -> DOFArray:
-        r"""Return the solid density $\rho_s$."""
-        if isinstance(material_densities, DOFArray):
-            return material_densities
-        return sum(material_densities)
+    def density(self) -> DOFArray:
+        """Return the wall density for all components."""
+        return self._density_func()
 
-    def decomposition_progress(self, material_densities) -> DOFArray:
-        r"""Evaluate the progress ratio $\tau$ of the decomposition."""
-        mass = self.solid_density(material_densities)
-        if self._decomposition_func:
-            # for materials that decompose when heated
-            return self._decomposition_func(mass)
-        actx = mass.array_context
-        return actx.np.zeros_like(mass) + 1.0
+    def heat_capacity(self, temperature: DOFArray = None) -> DOFArray:
+        """Return the wall heat_capacity for all components."""
+        return self._heat_capacity_func(temperature)
 
-    def heat_capacity(self, temperature: DOFArray = None, tau=None) -> DOFArray:
-        """Return the wall heat_capacity."""
-        return self._heat_capacity_func(temperature=temperature, tau=tau)
+    def enthalpy(self, temperature: DOFArray) -> DOFArray:
+        """Return the wall enthalpy for all components."""
+        return self._enthalpy_func(temperature)
 
-    def enthalpy(self, temperature: DOFArray, tau=None) -> DOFArray:
-        """Return the wall enthalpy."""
-        return self._enthalpy_func(temperature=temperature, tau=tau)
-
-    def thermal_diffusivity(self, solid_state: SolidWallState) -> DOFArray:
+    def thermal_diffusivity(self, mass: DOFArray, temperature: DOFArray,
+            thermal_conductivity: DOFArray = None) -> DOFArray:
         """Return the wall thermal diffusivity for all components."""
-        tau = self.decomposition_progress(solid_state.cv.mass)
-        mass = self.solid_density(solid_state.cv.mass)
-        kappa = self.thermal_conductivity(solid_state.dv.temperature, tau)
-        cp = self.heat_capacity(solid_state.dv.temperature, tau)
-        return kappa/(mass * cp)
+        if thermal_conductivity is None:
+            thermal_conductivity = self.thermal_conductivity(temperature)
+        return thermal_conductivity/(mass * self.heat_capacity(temperature))
 
-    def thermal_conductivity(self, temperature: DOFArray, tau=None) -> DOFArray:
-        """Return the wall thermal conductivity."""
-        return self._thermal_conductivity_func(temperature=temperature, tau=tau)
+    def thermal_conductivity(self, temperature: DOFArray) -> DOFArray:
+        """Return the wall thermal conductivity for all components."""
+        return self._thermal_conductivity_func(temperature)
 
     def get_temperature(self, wv: SolidWallConservedVars,
             tseed: DOFArray = None, niter: int = 3) -> DOFArray:
         """Evaluate the temperature based on the energy."""
         if tseed is not None:
-            tau = self.decomposition_progress(wv.mass)
             temp = tseed*1.0
             for _ in range(0, niter):
-                h = self.enthalpy(temperature=temp, tau=tau)
-                cp = self.heat_capacity(temperature=temp, tau=tau)
-                temp = temp - (h - wv.energy/self.solid_density(wv.mass))/cp
+                h = self.enthalpy(temp)
+                cp = self.heat_capacity(temp)
+                temp = temp - (h - wv.energy/wv.mass)/cp
             return temp
 
-        return wv.energy/(self.solid_density(wv.mass)*self.heat_capacity())
+        return wv.energy/(self.density()*self.heat_capacity())
 
     def dependent_vars(self, wv: SolidWallConservedVars,
             tseed: DOFArray = None, niter: int = 3) -> SolidWallDependentVars:
         """Return solid wall dependent variables."""
-        tau = self.decomposition_progress(wv.mass)
-        density = self.solid_density(wv.mass)
         temperature = self.get_temperature(wv, tseed, niter)
-        kappa = self.thermal_conductivity(temperature=temperature, tau=tau)
+        kappa = self.thermal_conductivity(temperature)
         return SolidWallDependentVars(
-            tau=tau,
-            density=density,
             thermal_conductivity=kappa,
             temperature=temperature)
 
@@ -438,12 +414,6 @@ def solid_density(self, material_densities) -> DOFArray:
         if isinstance(material_densities, DOFArray):
             return material_densities
         return sum(material_densities)
-        # if isinstance(material_densities, DOFArray):
-        #     actx = material_densities.array_context
-        #     return actx.np.absolute(material_densities)
-        # else:
-        #     actx = material_densities[0].array_context
-        #     return sum(actx.np.absolute(material_densities))
 
     def decomposition_progress(self, material_densities) -> DOFArray:
         r"""Evaluate the progress ratio $\tau$ of the decomposition.
@@ -522,35 +492,3 @@ def heat_capacity(self, cv: ConservedVars, wv: PorousWallVars,
         """
         return (cv.mass*self.eos.heat_capacity_cv(cv, temperature)
                 + wv.density*self.wall_eos.heat_capacity(temperature, wv.tau))
-
-    # def thermal_diffusivity(self, fluid_state) -> DOFArray:
-    #     """Return the wall thermal diffusivity for all components."""
-    #     tau = self.decomposition_progress(solid_state.cv.mass)
-    #     mass = self.solid_density(solid_state.cv.mass)
-    #     kappa = fluid_state.tv.thermal_conductivity
-    #     cp = self.heat_capacity(fluid_state.cv, fluid_state.wv,
-    #                             fluid_state.dv.temperature)
-    #     return kappa/(mass * cp)
-
-
-def get_porous_flow_timestep(dcoll, gas_model, state, cfl, dd):
-    r"""Routine returns the the node-local maximum stable viscous timestep."""
-    actx = state.array_context
-
-    length_scales = characteristic_lengthscales(
-        state.array_context, dcoll, dd=dd)
-
-    kappa = reduce(actx.np.maximum, state.thermal_conductivity)
-    mass_cp = gas_model.heat_capacity(state.cv, state.wv, state.dv.temperature)
-
-    nu = state.viscosity / state.mass_density
-    alpha = kappa / mass_cp
-    d_species_max = \
-        get_local_max_species_diffusivity(actx, state.species_diffusivity)
-    viscous_stuff = nu + alpha + d_species_max
-    vdt = op.elementwise_min(
-        dcoll, dd,
-        length_scales / (state.wavespeed + ((viscous_stuff) / length_scales))
-    )
-
-    return cfl * vdt

From 057b2f7654a0c7af3a9ecddec3eaf3454d18c1db Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 15 May 2024 00:08:41 -0500
Subject: [PATCH 45/48] Remove updates to carbon --  that will be added as a
 new PR

---
 doc/misc.rst                       |   2 -
 mirgecom/materials/carbon_fiber.py | 191 +++++------------------------
 mirgecom/materials/tacot.py        |   8 +-
 test/test_wallmodel.py             |   5 -
 4 files changed, 36 insertions(+), 170 deletions(-)

diff --git a/doc/misc.rst b/doc/misc.rst
index 33d1e466d..df4c005a0 100644
--- a/doc/misc.rst
+++ b/doc/misc.rst
@@ -96,5 +96,3 @@ References
 .. [Ern_2008] Daniele A. Di Pietro, Alexandre Ern, Jean-Luc Guermond, Discontinuous Galerkin Methods for \
    Anisotropic Semidefinite Diffusion with Advection, SIAM Journal on Numerical Analysis 46 2 \
    `(DOI) <https://www.jstor.org/stable/40233233?seq=12>`__
-.. [Martin_2013] Alexandre Martin, Volume averaged modeling of the oxidation of porous carbon fiber material,
-    AIAA Thermophysics conference 2013-2636
diff --git a/mirgecom/materials/carbon_fiber.py b/mirgecom/materials/carbon_fiber.py
index b15610db3..d1e984cde 100644
--- a/mirgecom/materials/carbon_fiber.py
+++ b/mirgecom/materials/carbon_fiber.py
@@ -2,8 +2,6 @@
 
 .. autoclass:: Oxidation
 .. autoclass:: Y2_Oxidation_Model
-.. autoclass:: Y3_Oxidation_Model
-.. autoclass:: OxidationModel
 .. autoclass:: FiberEOS
 """
 
@@ -31,57 +29,60 @@
 THE SOFTWARE.
 """
 
-from typing import Optional, Union
+from typing import Optional
 from abc import abstractmethod
 import numpy as np
-from pytools.obj_array import make_obj_array
 from meshmode.dof_array import DOFArray
 from mirgecom.wall_model import PorousWallEOS
+from pytools.obj_array import make_obj_array
 
 
 class Oxidation:
-    """Abstract interface for wall oxidation models.
+    """Abstract interface for wall oxidation model.
 
     .. automethod:: get_source_terms
     """
 
     @abstractmethod
-    def get_source_terms(
-            self, temperature: DOFArray, tau: DOFArray, rhoY_o2: DOFArray,  # noqa N803
-            rhoY_o: Optional[DOFArray] = None) -> Union[DOFArray, tuple]:
-        r"""Source terms of fiber oxidation.
-
-        Parameters
-        ----------
-        temperature:
-            gas and solid temperature assuming local-equilibrium
-        tau:
-            the progress ratio of the oxidation: 1 for virgin, 0 for fully oxidized
-        rhoY_o2:
-            the mass fraction of molecular oxygen
-        rhoY_o:
-            the mass fraction of atomic oxygen
-        """
+    def get_source_terms(self, temperature: DOFArray,
+            tau: DOFArray, rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+        r"""Source terms of fiber oxidation."""
         raise NotImplementedError()
 
 
+# TODO per MTC review, can we generalize the oxidation model?
+# should we keep this in the driver?
+# https://github.com/illinois-ceesd/mirgecom/pull/875#discussion_r1261414281
 class Y2_Oxidation_Model(Oxidation):  # noqa N801
-    r"""Evaluate the source terms for carbon fiber oxidation in Y2 model.
+    """Evaluate the source terms for the Y2 model of carbon fiber oxidation.
 
+    .. automethod:: puma_effective_surface_area
     .. automethod:: get_source_terms
     """
 
-    def _get_wall_effective_surface_area_fiber(self, tau: DOFArray) -> DOFArray:
-        """Evaluate the effective surface of the fibers."""
-        # Polynomial fit based on PUMA data:
+    def puma_effective_surface_area(self, tau: DOFArray) -> DOFArray:
+        """Polynomial fit based on PUMA data."""
         # Original fit function: -1.1012e5*x**2 - 0.0646e5*x + 1.1794e5
         # Rescale by x==0 value and rearrange
         progress = 1.0-tau
         return 1.1794e5*(1.0 - 0.0547736137*progress - 0.9336950992*progress**2)
 
+    def _get_wall_effective_surface_area_fiber(self, tau: DOFArray) -> DOFArray:
+        """Evaluate the effective surface of the fibers."""
+        return self.puma_effective_surface_area(tau)
+
     def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-                         rhoY_o2: DOFArray, rhoY_o=None) -> DOFArray:  # noqa N803
-        """Return the effective source terms for fiber oxidation."""
+                         rhoY_o2: DOFArray) -> DOFArray:  # noqa N803
+        """Return the effective source terms for fiber oxidation.
+
+        Parameters
+        ----------
+        temperature:
+        tau:
+            the progress ratio of the oxidation: 1 for virgin, 0 for fully oxidized
+        rhoY_o2:
+            the mass fraction of oxygen
+        """
         actx = temperature.array_context
 
         mw_o = 15.999
@@ -98,120 +99,6 @@ def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
         return (mw_co/mw_o2 + mw_o/mw_o2 - 1)*rhoY_o2*k*eff_surf_area
 
 
-class Y3_Oxidation_Model(Oxidation):  # noqa N801
-    r"""Evaluate the source terms for carbon fiber oxidation in Y3 model.
-
-    Follows [Martin_2013]_ by using a single reaction given by
-
-    .. math::
-        C_{(s)} + O_2 \to CO_2
-
-    The degradation is computed as
-
-    .. math::
-        \frac{\partial \rho_s}{\partial t} =
-        \rho_C \frac{\partial \epsilon}{\partial t} = \dot{\omega}_\text{het}
-
-    The fiber radius as a function of mass loss $\tau$ is given by
-
-    .. math::
-        \tau = \frac{m}{m_0} = \frac{\epsilon}{\epsilon_0} =
-        \frac{\pi r^2/L}{\pi r_0^2/L} = \frac{r^2}{r_0^2}
-
-    and leads to the wall recession as
-
-    .. math::
-        \frac{\partial \epsilon}{\partial t} =
-        \frac{\partial \epsilon}{\partial r} \frac{\partial r}{\partial t} =
-        2 \frac{\epsilon_0}{r_0^2} r = 2 \frac{\epsilon_0}{r_0} \sqrt{\tau}
-
-    .. automethod:: __init__
-    .. automethod:: get_source_terms
-    """
-
-    def __init__(self, wall_material, initial_fiber_radius=5.5e-6, arrhenius=1e5,
-                 activation_energy=-120000.0):
-        self._material = wall_material
-        self._init_fiber_radius = initial_fiber_radius
-        self._arrhenius = arrhenius
-        self._Ea = activation_energy
-
-    def _get_wall_effective_surface_area_fiber(self, tau) -> DOFArray:
-        r"""Evaluate the effective surface of the fibers."""
-        actx = tau.array_context
-        zeros = actx.np.zeros_like(tau)
-        ones = zeros + 1.0
-        return (
-            2.0 * self._material.volume_fraction(tau=ones)/self._init_fiber_radius
-            * actx.np.sqrt(actx.np.maximum(tau, zeros))  # only use positive tau
-        )
-
-    def get_source_terms(self, temperature: DOFArray, tau: DOFArray,
-            rhoY_o2: DOFArray, rhoY_o=None):  # noqa N803
-        r"""Return the effective source terms for the oxidation.
-
-        Returns
-        -------
-        sources: tuple
-            the tuple ($\omega_{C}$, $\omega_{O_2}$, $\omega_{CO_2}$)
-        """
-        actx = temperature.array_context
-
-        # ignore negative amounts of O2
-        rhoY_o2 = actx.np.maximum(rhoY_o2, actx.np.zeros_like(rhoY_o2))  # noqa N806
-
-        mw_c = 12.011
-        mw_o = 15.999
-        mw_o2 = mw_o*2
-        mw_co2 = 44.010
-        univ_gas_const = 8.31446261815324  # J/(K-mol)
-
-        eff_surf_area = self._get_wall_effective_surface_area_fiber(tau)
-
-        k_f = self._arrhenius*actx.np.exp(self._Ea/(univ_gas_const*temperature))
-
-        m_dot_c = - rhoY_o2/mw_o2 * mw_c * eff_surf_area * k_f
-        m_dot_o2 = - rhoY_o2/mw_o2 * mw_o2 * eff_surf_area * k_f
-        m_dot_co2 = + rhoY_o2/mw_o2 * mw_co2 * eff_surf_area * k_f
-
-        return m_dot_c, m_dot_o2, m_dot_co2
-
-
-class OxidationModel(Oxidation):
-    """Evaluate the source terms for the carbon fiber oxidation.
-
-    The user must specify the functions for oxidation in the driver.
-    (Tentatively) Generalizing this class makes it easier for adding new,
-    more complex models and for UQ runs.
-
-    .. automethod:: get_source_terms
-    """
-
-    def __init__(self, surface_area_func, oxidation_func):
-        """Initialize the user-specified oxidation class.
-
-        Parameters
-        ----------
-        surface_area_func:
-            Function prescribing how the fiber area changes during oxidation.
-        oxidation_func:
-            Reaction rate for the oxidation model.
-        """
-        self._surface_func = surface_area_func
-        self._oxidation_func = oxidation_func
-
-    def get_source_terms(
-            self, temperature: DOFArray, tau: DOFArray, rhoY_o2: DOFArray,  # noqa N803
-            rhoY_o: Optional[DOFArray] = None) -> DOFArray:
-        """Return the effective source terms for the oxidation."""
-        area = self._surface_func(tau)
-        if rhoY_o is None:
-            return self._oxidation_func(temperature=temperature, fiber_area=area,
-                                        rhoY_o2=rhoY_o2)
-        return self._oxidation_func(temperature=temperature, fiber_area=area,
-                                    rhoY_o2=rhoY_o2, rhoY_o=rhoY_o)
-
-
 class FiberEOS(PorousWallEOS):
     r"""Evaluate the properties of the solid state containing only fibers.
 
@@ -235,8 +122,7 @@ class FiberEOS(PorousWallEOS):
     .. automethod:: decomposition_progress
     """
 
-    def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass,
-                 timescale=1.0):
+    def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass):
         """Bulk density considering the porosity and intrinsic density.
 
         Parameters
@@ -250,15 +136,11 @@ def __init__(self, dim, anisotropic_direction, char_mass, virgin_mass,
             final mass when the decomposition is complete.
         virgin_mass: float
             initial mass of the material.
-        timescale: float
-            Modifies the thermal conductivity and the radiation emission to
-            increase/decrease the wall time-scale. Defaults to 1.0 (no changes).
         """
         self._char_mass = char_mass
         self._virgin_mass = virgin_mass
         self._dim = dim
         self._anisotropic_dir = anisotropic_direction
-        self._timescale = timescale
 
         if anisotropic_direction >= dim:
             raise ValueError("Anisotropic axis must be less than dim.")
@@ -308,21 +190,18 @@ def thermal_conductivity(self, temperature, tau) -> np.ndarray:
             + 1.93072961e-10*temperature**3 - 3.52595953e-07*temperature**2
             + 4.54935976e-04*temperature**1 + 5.08960039e-02)
 
-        # initialize with the in-plane value then modify the normal direction
+        # initialize with the in-plane value
         kappa = make_obj_array([kappa_ij for _ in range(self._dim)])
+        # modify only the normal direction
         kappa[self._anisotropic_dir] = kappa_k
 
         # account for fiber shrinkage via "tau"
-        # XXX check if here is the best place for timescale
-        return kappa*tau*self._timescale
+        return kappa*tau
 
     # ~~~~~~~~ other properties
     def volume_fraction(self, tau: DOFArray) -> DOFArray:
         r"""Fraction $\phi$ occupied by the solid."""
-        # XXX check if this keeps tau positive and bounds everything else to physical
-        # values (it doesn't crash, but it is non-ideal to have unphysical values)
-        actx = tau.array_context
-        return actx.np.maximum(0.12*tau, 0.0*tau)
+        return 0.12*tau
 
     def permeability(self, tau: DOFArray) -> np.ndarray:
         r"""Permeability $K$ of the porous material."""
@@ -332,14 +211,12 @@ def permeability(self, tau: DOFArray) -> np.ndarray:
         permeability = make_obj_array([5.57e-11 + actx.np.zeros_like(tau)
                                        for _ in range(0, self._dim)])
         permeability[self._anisotropic_dir] = 2.62e-11 + actx.np.zeros_like(tau)
-
         return permeability
 
     def emissivity(self, temperature: DOFArray,  # type: ignore[override]
                    tau: Optional[DOFArray] = None) -> DOFArray:
         """Emissivity for energy radiation."""
-        # XXX check if here is the best place for timescale
-        return self._timescale * (
+        return (
             + 2.26413679e-18*temperature**5 - 2.03008004e-14*temperature**4
             + 7.05300324e-11*temperature**3 - 1.22131715e-07*temperature**2
             + 1.21137817e-04*temperature**1 + 8.66656964e-01)
diff --git a/mirgecom/materials/tacot.py b/mirgecom/materials/tacot.py
index 25f027c51..01586a253 100644
--- a/mirgecom/materials/tacot.py
+++ b/mirgecom/materials/tacot.py
@@ -120,12 +120,8 @@ def get_source_terms(self, temperature, chi):
             actx.np.where(actx.np.less(temperature, self._Tcrit[1]),
                 0.0, (-w2 * self._pre_exp[1] * actx.np.exp(-20444.44/temperature))
                 ),
-            # TODO: fiber oxidation:
-            # Include a resin-dependent term (probably linear) to slow down
-            # oxidation when there is a lot of resin, but make it goes with
-            # the "right" rate when the resin is consumed.
-            actx.np.zeros_like(temperature)
-            ])
+            # fiber oxidation: include in the RHS but don't do anything with it.
+            actx.np.zeros_like(temperature)])
 
 
 class TacotEOS(PorousWallEOS):
diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index 441db76c6..eec3f49c3 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -205,8 +205,3 @@ def test_tacot_decomposition(actx_factory):
         op.norm(dcoll, solid_mass_rhs[1] + 2.03565420370596, np.inf)) < tol
     assert actx.to_numpy(
         op.norm(dcoll, sample_source_gas - 28.8032660577024, np.inf)) < tol
-
-
-# TODO
-def test_oxidation_model():
-    assert 1 == 1

From 1999b1488e1809e4703e1601436f6872719b0f5d Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 15 May 2024 08:45:41 -0500
Subject: [PATCH 46/48] Remove useless changes

---
 mirgecom/eos.py        | 13 -------------
 mirgecom/transport.py  |  6 ++----
 mirgecom/wall_model.py |  4 ++--
 3 files changed, 4 insertions(+), 19 deletions(-)

diff --git a/mirgecom/eos.py b/mirgecom/eos.py
index bdddee2a8..e9f447c93 100644
--- a/mirgecom/eos.py
+++ b/mirgecom/eos.py
@@ -870,16 +870,3 @@ def get_species_source_terms(self, cv: ConservedVars, temperature: DOFArray):
 
         return make_conserved(dim, rho_source, energy_source, mom_source,
                               species_sources)
-
-    def get_mole_fractions(self, species_mass_fractions):
-        """Get mole fractions from mass fractions."""
-        mix_mol_weight = \
-            self._pyrometheus_mech.get_mix_molecular_weight(species_mass_fractions)
-        return self._pyrometheus_mech.get_mole_fractions(mix_mol_weight,
-                                                         species_mass_fractions)
-
-    def get_mass_fractions(self, species_mole_fractions):
-        """Get mass fractions from mole fractions."""
-        mol_weights = self.get_species_molecular_weights()
-        mmw = sum(species_mole_fractions*mol_weights)
-        return species_mole_fractions*mol_weights/mmw
diff --git a/mirgecom/transport.py b/mirgecom/transport.py
index 8d10b181c..8d7a7137b 100644
--- a/mirgecom/transport.py
+++ b/mirgecom/transport.py
@@ -315,11 +315,9 @@ def species_diffusivity(self, cv: ConservedVars,  # type: ignore[override]
             d_{\alpha} = \frac{\kappa}{\rho \; Le \; C_p}
         """
         if self._lewis is not None:
-            alpha = (
-                self._sigma * self.viscosity(cv, dv)
-                / (cv.mass * eos.gamma(cv, dv.temperature))
+            return (self._sigma * self.viscosity(cv, dv)/(
+                cv.mass*self._lewis*eos.gamma(cv, dv.temperature))
             )
-            return self._lewis*alpha
         return self._d_alpha*(0*cv.mass + 1.)
 
 
diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 779a6dfba..ebb82710c 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -473,7 +473,7 @@ def get_pressure(self, cv: ConservedVars, wv: PorousWallVars,
         return cv.mass*self.eos.gas_const()*temperature/wv.void_fraction
 
     def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
-                        temperature: DOFArray) -> DOFArray:
+                 temperature: DOFArray) -> DOFArray:
         r"""Return the internal energy of the gas+solid material.
 
         .. math::
@@ -484,7 +484,7 @@ def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
                 + wv.density*self.wall_eos.enthalpy(temperature, wv.tau))
 
     def heat_capacity(self, cv: ConservedVars, wv: PorousWallVars,
-                      temperature: DOFArray) -> DOFArray:
+                 temperature: DOFArray) -> DOFArray:
         r"""Return the heat capacity of the gas+solid material.
 
         .. math::

From 5a6842228173ae703d7e26c35d9757f56b8848cd Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Wed, 2 Oct 2024 17:45:40 -0500
Subject: [PATCH 47/48] flake8, pylint, mypy

---
 mirgecom/wall_model.py |  3 ++-
 test/test_wallmodel.py | 10 ++++++----
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index 4336f88ae..a30b80443 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -474,7 +474,8 @@ def get_pressure(self, cv: ConservedVars, wv: PorousWallVars,
         """
         if isinstance(self.eos, MixtureEOS):
             return self.eos.pressure(cv, temperature)/wv.void_fraction
-        return cv.mass*self.eos.gas_const()*temperature/wv.void_fraction
+        return 1.0/wv.void_fraction*(
+            cv.mass*self.eos.gas_const()*temperature)  # type: ignore
 
     def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
                  temperature: DOFArray) -> DOFArray:
diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index 10392bbad..d0bf9495d 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -43,6 +43,10 @@
 from mirgecom.thermochemistry import get_pyrometheus_wrapper_class_from_cantera
 
 
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+    [PytestPyOpenCLArrayContextFactory])
+
+
 @pytest.mark.parametrize("order", [1, 4])
 @pytest.mark.parametrize("my_material", ["fiber", "composite"])
 @pytest.mark.parametrize("my_eos", ["ideal", "mixture"])
@@ -82,7 +86,7 @@ def test_wall_eos(actx_factory, order, my_material, my_eos, use_diffused_interfa
 
         eos = PyrometheusMixture(pyro_obj, temperature_guess=cantera_soln.T)
 
-    if my_eos == "ideal":
+    else:
         eos = IdealSingleGas()
         y = None
 
@@ -90,6 +94,7 @@ def test_wall_eos(actx_factory, order, my_material, my_eos, use_diffused_interfa
 
     # {{{ Initialize wall model
 
+    material = None
     if my_material == "fiber":
         import mirgecom.materials.carbon_fiber as material_sample
         material = material_sample.FiberEOS(
@@ -152,9 +157,6 @@ def test_wall_eos(actx_factory, order, my_material, my_eos, use_diffused_interfa
     assert actx.to_numpy(
         op.norm(dcoll, solid_state.temperature - 900.0, np.inf)) < tol
 
-pytest_generate_tests = pytest_generate_tests_for_array_contexts(
-    [PytestPyOpenCLArrayContextFactory])
-
 
 def test_tacot_decomposition(actx_factory):
     """Check the wall degradation model."""

From 0618d14d6f587b8e735df8e900f81ee462c3ac7a Mon Sep 17 00:00:00 2001
From: Tulio <tulioricci@gmail.com>
Date: Thu, 3 Oct 2024 08:26:24 -0500
Subject: [PATCH 48/48] mypy; fix np.max in test_wallmodel

---
 mirgecom/wall_model.py | 2 +-
 test/test_wallmodel.py | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/mirgecom/wall_model.py b/mirgecom/wall_model.py
index a30b80443..c31112dec 100644
--- a/mirgecom/wall_model.py
+++ b/mirgecom/wall_model.py
@@ -474,7 +474,7 @@ def get_pressure(self, cv: ConservedVars, wv: PorousWallVars,
         """
         if isinstance(self.eos, MixtureEOS):
             return self.eos.pressure(cv, temperature)/wv.void_fraction
-        return 1.0/wv.void_fraction*(
+        return 1.0/wv.void_fraction*(  # type: ignore
             cv.mass*self.eos.gas_const()*temperature)  # type: ignore
 
     def internal_energy(self, cv: ConservedVars, wv: PorousWallVars,
diff --git a/test/test_wallmodel.py b/test/test_wallmodel.py
index d0bf9495d..54e252ddc 100644
--- a/test/test_wallmodel.py
+++ b/test/test_wallmodel.py
@@ -145,12 +145,12 @@ def test_wall_eos(actx_factory, order, my_material, my_eos, use_diffused_interfa
     assert isinstance(solid_state.wv.density, DOFArray)
 
     if my_material == "fiber":
-        assert np.max(actx.to_numpy(wv.density - 168.0)) < tol
-        assert np.max(actx.to_numpy(wv.void_fraction)) - 1.00 < tol
+        assert actx.to_numpy(actx.np.max(wv.density - 168.0)) < tol
+        assert actx.to_numpy(actx.np.max(wv.void_fraction)) - 1.00 < tol
 
     if my_material == "composite":
-        assert np.max(actx.to_numpy(wv.density - 280.0)) < tol
-        assert np.max(actx.to_numpy(wv.void_fraction)) - 1.00 < tol
+        assert actx.to_numpy(actx.np.max(wv.density - 280.0)) < tol
+        assert actx.to_numpy(actx.np.max(wv.void_fraction)) - 1.00 < tol
 
     assert actx.to_numpy(
         op.norm(dcoll, solid_state.pressure - 101325.0, np.inf)) < tol