Add chem_params%transport_dt: limit dt by the local chemistry transport

docs/documentation/case.md

@@ -1256,6 +1256,7 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 | `chem_params%%transport_model` | Integer | Methodology for calculating the diffusion coefficients   |
 | `chem_params%%reaction_substeps` | Integer | Sub-steps for operator-split reaction integration (0 = off) |
 | `chem_params%%adap_substeps`   | Logical | Per-rank adaptive sub-step count driven by local stiffness  |
+| `chem_params%%transport_dt`    | Logical | Limit the adaptive dt by the local chemistry transport      |
 | `chem_params%%reaction_substeps_max` | Integer | Sub-step ceiling when `adap_substeps` is enabled       |
 | `cantera_file`                | String  | Cantera-format mechanism file (e.g., .yaml)              |
 | `surface_cantera_file`        | String  | Cantera-format mechanism file for heterogeneous surface chemistry |
@@ -1265,6 +1266,8 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 - `chem_params%%reaction_substeps` controls how the reaction source is integrated. With `0` (default) the net production rates are added to the flow right-hand side and advanced by the flow time stepper (fine for hydrogen). With a value `> 0`, the reaction is instead integrated by operator splitting after each flow update: every cell's constant-density, constant-internal-energy reactor is advanced over the timestep with that many sub-steps of an **α-QSS** (quasi-steady-state) integrator — a matrix-free, Jacobian-free predictor–corrector (Mott/CHEMEQ2) that splits the net rate into creation/destruction parts and applies a Padé α-weighting, so it stays stable on stiff mechanisms where an explicit source diverges. This decouples the (often much faster) chemical timescale from the flow timestep and is required for stiff mechanisms — e.g. hydrocarbons such as GRI-Mech methane, which otherwise diverge on the first step
 - `chem_params%%adap_substeps` (default `F`) makes each rank choose its α-QSS sub-step count per flow step from the largest chemical stiffness among its own cells: the count sits at `reaction_substeps` (the floor) in inert or burned gas and rises toward `reaction_substeps_max` (the ceiling) only across the reaction front. It uses no MPI collectives. When enabled, `reaction_substeps >= 1` and `reaction_substeps_max >= reaction_substeps` are required
 
+- `chem_params%%transport_dt` (default `F`) replaces the constant `fluid_pp%%Re` in the viscous time-step limit with the local mixture transport: dt <= `cfl_target`·Δx²/(d·D), with d the number of dimensions (explicit diffusion is stable while the d per-direction diffusion numbers sum below the limit) and D the largest of μ/ρ (with `viscous`), and λ/(ρ c_v) and the mixture-averaged species diffusivities (with `chem_params%%diffusion`); c_v, not c_p, because a cell cannot expand within a step. Mixture transport grows steeply with temperature, so near hot walls and flames a dt limited by the free-stream viscosity can exceed the explicit diffusion limit. Requires `cfl_adap_dt` or `cfl_const_dt`, and `viscous` or `chem_params%%diffusion`. A non-Newtonian `mu_max` or bulk viscosity (`fluid_pp%%Re(2)`) still bounds dt when it is tighter. Not supported with `igr`
+
 - `cantera_file` specifies the chemical mechanism file. If the file is part of the standard Cantera library, only the filename is required. Otherwise, the file must be located in the same directory as your `case.py` file
 
 - `surface_cantera_file` and `surface_phase` specify the Cantera mechanism file and interface phase used for heterogeneous surface chemistry. These parameters must be specified together when a surface mechanism is used.

src/common/m_chemistry.fpp

@@ -13,7 +13,7 @@ module m_chemistry
         & get_creation_destruction_rates, get_mole_fractions, get_species_binary_mass_diffusivities, &
         & get_species_mass_diffusivities_mixavg, gas_constant, get_mixture_molecular_weight, get_mixture_energy_mass, &
         & get_mixture_thermal_conductivity_mixavg, get_species_enthalpies_rt, get_mixture_viscosity_mixavg, &
-        & get_mixture_specific_heat_cp_mass, get_mixture_enthalpy_mass
+        & get_mixture_specific_heat_cp_mass, get_mixture_specific_heat_cv_mass, get_mixture_enthalpy_mass
     use m_thermochem_state, only: get_mixavg_transport_state, get_species_enthalpies_mass
     use m_global_parameters
 
@@ -42,6 +42,39 @@ contains
 
     end subroutine compute_viscosity_and_inversion
 
+    !> Effective Reynolds numbers 1/(rho*D) from the largest local transport diffusivity D -- mu/rho with viscous, and
+    !! lambda/(rho*cv) and the mixture-averaged species diffusivities with chem_params%diffusion -- so the viscous dt limit
+    !! (dx^2*rho*Re) tracks the transport actually integrated, summed over the num_dims directions.
+    subroutine s_compute_transport_dt_re(pres, rho, T, Ys, Re)
+
+        $:GPU_ROUTINE(function_name='s_compute_transport_dt_re',parallelism='[seq]', cray_inline=True)
+
+        real(wp), intent(in)                         :: pres, rho, T
+        real(wp), dimension(num_species), intent(in) :: Ys
+        real(wp), dimension(2), intent(out)          :: Re
+        real(wp)                                     :: D_max, prop, cv
+        real(wp), dimension(num_species)             :: D_k
+
+        D_max = 0._wp
+        if (viscous) then
+            call get_mixture_viscosity_mixavg(T, Ys, prop)
+            D_max = prop/rho
+        end if
+        if (chem_params%diffusion) then
+            ! cv, not cp: a cell cannot expand within a step, so the grid-scale mode of T diffuses at lambda/(rho*cv)
+            call get_mixture_thermal_conductivity_mixavg(T, Ys, prop)
+            call get_mixture_specific_heat_cv_mass(T, Ys, cv)
+            D_max = max(D_max, prop/(rho*cv))
+            if (chem_params%transport_model == 1) then
+                call get_species_mass_diffusivities_mixavg(pres, T, Ys, D_k)
+                D_max = max(D_max, maxval(D_k))
+            end if
+        end if
+        ! The diffusion numbers of the num_dims directions add up
+        Re(:) = 1._wp/(rho*real(num_dims, wp)*max(D_max, sgm_eps))
+
+    end subroutine s_compute_transport_dt_re
+
     !> Initialize the temperature field from conservative variables by inverting the energy equation.
     subroutine s_compute_q_T_sf(q_T_sf, q_cons_vf, bounds)
 

src/common/m_derived_types.fpp

@@ -567,6 +567,8 @@ module m_derived_types
         !> Zero MPI: each rank sizes its own work from its own cells. Default F = fixed reaction_substeps.
         logical :: adap_substeps
         integer :: reaction_substeps_max
+        !> Limit dt by the mixture's local viscous, thermal and species diffusivities rather than fluid_pp%Re.
+        logical :: transport_dt
     end type chemistry_parameters
 
     !> Condensed-phase reactive-burn (programmed pressure detonation) parameters. The rate is

src/post_process/m_global_parameters.fpp

@@ -205,6 +205,7 @@ contains
         chem_params%reaction_substeps = 0
         chem_params%adap_substeps = .false.
         chem_params%reaction_substeps_max = 0
+        chem_params%transport_dt = .false.
 
         ! Fluids physical parameters (post-specific; G = dflt_real differs from pre/sim)
         do i = 1, num_fluids_max

src/pre_process/m_global_parameters.fpp

@@ -422,6 +422,7 @@ contains
         chem_params%reaction_substeps = 0
         chem_params%adap_substeps = .false.
         chem_params%reaction_substeps_max = 0
+        chem_params%transport_dt = .false.
 
         ! Fluids physical parameters
         do i = 1, num_fluids_max

src/simulation/m_global_parameters.fpp

@@ -423,6 +423,7 @@ contains
         chem_params%reaction_substeps = 0
         chem_params%adap_substeps = .false.
         chem_params%reaction_substeps_max = 0
+        chem_params%transport_dt = .false.
 
         num_bc_patches = 0
         bc_io = .false.

src/simulation/m_sim_helpers.fpp

@@ -254,7 +254,7 @@ contains
         end if
 
         ! Viscous calculations
-        if (viscous) then
+        if (viscous .or. chem_params%transport_dt) then
             if (p > 0) then
                 if (grid_geometry == 3) then
                     fltr_dtheta = f_compute_filtered_dtheta(k, l)

src/simulation/m_time_steppers.fpp

@@ -671,22 +671,24 @@ contains
             real(wp), dimension(num_vels)   :: vel               !< Cell-avg. velocity
             real(wp), dimension(num_fluids) :: alpha, alpha_rho  !< Cell-avg. volume fraction, partial density
         #:endif
-        real(wp)               :: vel_sum            !< Cell-avg. velocity sum
-        real(wp)               :: pres               !< Cell-avg. pressure
-        real(wp)               :: gamma              !< Cell-avg. sp. heat ratio
-        real(wp)               :: pi_inf             !< Cell-avg. liquid stiffness function
-        real(wp)               :: qv                 !< Cell-avg. fluid reference energy
-        real(wp)               :: c                  !< Cell-avg. sound speed
-        real(wp), dimension(2) :: Re                 !< Cell-avg. Reynolds numbers
-        real(wp), dimension(4) :: max_dt             !< Cell dt candidates (inviscid, viscous, capillary, thermal)
-        real(wp)               :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local, tcfl_dt_local, coll_dt_local
-        real(wp), dimension(5) :: dt_candidates_loc  !< Rank-local dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap)
-        real(wp), dimension(5) :: dt_candidates_glb  !< Global dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap)
-        real(wp)               :: dt_prev
-        logical                :: is_fluid_cell      !< Cell lies outside every immersed boundary
-        integer                :: bad_cell           !< 1 + linear index of a fluid cell with no valid dt, 0 if none
-        integer                :: j, k, l            !< Generic loop iterators
-        integer                :: fl                 !< Fluid loop iterator
+        real(wp)                         :: vel_sum            !< Cell-avg. velocity sum
+        real(wp)                         :: pres               !< Cell-avg. pressure
+        real(wp)                         :: gamma              !< Cell-avg. sp. heat ratio
+        real(wp)                         :: pi_inf             !< Cell-avg. liquid stiffness function
+        real(wp)                         :: qv                 !< Cell-avg. fluid reference energy
+        real(wp)                         :: c                  !< Cell-avg. sound speed
+        real(wp), dimension(2)           :: Re                 !< Cell-avg. Reynolds numbers
+        real(wp), dimension(4)           :: max_dt             !< Cell dt candidates (inviscid, viscous, capillary, thermal)
+        real(wp)                         :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local, tcfl_dt_local, coll_dt_local
+        real(wp), dimension(5)           :: dt_candidates_loc  !< Rank-local dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap)
+        real(wp), dimension(5)           :: dt_candidates_glb  !< Global dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap)
+        real(wp)                         :: dt_prev
+        real(wp), dimension(num_species) :: Ys                 !< Cell mass fractions, for chem_params%transport_dt
+        real(wp), dimension(2)           :: Re_chem            !< Transport-limited Reynolds numbers, for chem_params%transport_dt
+        logical                          :: is_fluid_cell      !< Cell lies outside every immersed boundary
+        integer                          :: bad_cell           !< 1 + linear index of a fluid cell with no valid dt, 0 if none
+        integer                          :: j, k, l            !< Generic loop iterators
+        integer                          :: fl                 !< Fluid loop iterator
 
         if (.not. igr) then
             call s_convert_conservative_to_primitive_variables(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, idwint)
@@ -700,7 +702,7 @@ contains
         coll_dt_local = huge(1.0_wp)
         bad_cell = 0
         $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, alpha_rho, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, fl, &
-                            & max_dt, is_fluid_cell]', reduction='[[icfl_dt_local, vcfl_dt_local, ccfl_dt_local, &
+                            & max_dt, is_fluid_cell, Ys, Re_chem]', reduction='[[icfl_dt_local, vcfl_dt_local, ccfl_dt_local, &
                             & tcfl_dt_local], [bad_cell]]', reductionOp='[min, max]')
         do l = 0, p
             do k = 0, n
@@ -733,6 +735,21 @@ contains
                             Re(1) = 1._wp/max(Re(1), sgm_eps)
                         end if
 
+                        if (chemistry .and. chem_params%transport_dt) then
+                            $:GPU_LOOP(parallelism='[seq]')
+                            do fl = 1, num_species
+                                Ys(fl) = q_prim_vf(eqn_idx%species%beg + fl - 1)%sf(j, k, l)
+                            end do
+                            call s_compute_transport_dt_re(pres, rho, q_T_sf%sf(j, k, l), Ys, Re_chem)
+                            ! Keep the bounds Cantera does not replace: the HB mu_max in Re(1) and bulk viscosity in Re(2).
+                            ! Without viscous, Re holds dflt_real, not a bound
+                            if (viscous) then
+                                if (any_non_newtonian) Re_chem(1) = min(Re_chem(1), Re(1))
+                                Re_chem(2) = min(Re_chem(2), Re(2))
+                            end if
+                            Re = Re_chem
+                        end if
+
                         call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, alpha, alpha_rho, j, k, l)
 
                         ! min() drops a NaN candidate, so a cell with a NaN or non-positive dt would otherwise vanish from