Keep species mass fractions in [0, 1] by construction instead of clipping

docs/documentation/case.md

@@ -1262,6 +1262,7 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 | `surface_phase`               | String  | Cantera interface phase name for heterogeneous surface chemistry |
 
 - `chem_params%%transport_model` specifies the methodology for calculating diffusion coefficients and other transport properties, `1` for mixture-average, `2` for Unity-Lewis
+- With ``chemistry = 'T'``, species mass fractions are kept in [0, 1] by construction: reconstructed face values are limited onto the simplex (Zhang & Shu 2010), so the species fluxes sum to the mass flux (Larrouturou 1991), and with HLLC they are the upwind mass fractions times the mass flux (other Riemann solvers keep only the sum consistency); the diffusion correction velocity carries upwind mass fractions; and the α-QSS reaction substep keeps non-negative compositions non-negative without clipping. As a last resort, any cell still outside the bounds at the end of an RK stage is cleaned as in PeleC (clean_massfrac): partial densities clipped to [0, ρ], ρ reset to their sum, velocity and specific energy kept, with the number of fluid-cell cleanings (per cell per RK stage) printed. Cell averages provably stay non-negative when, in every cell, dt times the outgoing mass flux per unit length, summed over directions, is at most the reconstruction's Gauss–Lobatto endpoint weight (1/12 for WENO5, 1/6 for WENO3, 1/2 for MUSCL) times ρ; at usual CFL numbers, and for the diffusion terms, the limiter and upwinding are standard practice rather than a proof, and the fallback cleanup is the safety net.
 - `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
 

docs/module_categories.json

@@ -37,7 +37,8 @@
             "m_acoustic_src",
             "m_body_forces",
             "m_pressure_relaxation",
-            "m_collisions"
+            "m_collisions",
+            "m_species_bounds"
         ]
     },
     {

src/common/m_chemistry.fpp

@@ -160,7 +160,7 @@ contains
         type(int_bounds_info), dimension(1:3), intent(in)      :: bounds
         integer                                                :: x, y, z, eqn, s, nsub
         real(wp)                                               :: rho, energy, T, T_new, dt_sub, Ysum
-        real(wp)                                               :: r, r2, wr, loss_i, prod_p, loss_p, Lbar, pbar
+        real(wp)                                               :: wr, loss_i, prod_p, loss_p, Lbar, pbar
         real(wp)                                               :: stiff_max, cell_stiff
         real(wp), parameter                                    :: y_floor = 1.e-16_wp
         ! stiff_target: fractional net composition change per sub-step targeted when sizing the adaptive
@@ -219,8 +219,8 @@ contains
         end if
         dt_sub = dtime/real(nsub, wp)
 
-        $:GPU_PARALLEL_LOOP(collapse=3, private='[Ys, cdot, ddot, y0, prod0, Lloss, alp, eqn, s, rho, energy, T, T_new, Ysum, r, &
-                            & r2, wr, loss_i, prod_p, loss_p, Lbar, pbar]', copyin='[bounds, dt_sub, nsub]')
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[Ys, cdot, ddot, y0, prod0, Lloss, alp, eqn, s, rho, energy, T, T_new, Ysum, wr, &
+                            & loss_i, prod_p, loss_p, Lbar, pbar]', copyin='[bounds, dt_sub, nsub]')
         do z = bounds(3)%beg, bounds(3)%end
             do y = bounds(2)%beg, bounds(2)%end
                 do x = bounds(1)%beg, bounds(1)%end
@@ -254,10 +254,8 @@ contains
                             prod0(eqn) = wr*cdot(eqn)  ! mass-fraction production
                             loss_i = wr*ddot(eqn)  ! mass-fraction loss
                             Lloss(eqn) = loss_i/max(Ys(eqn), y_floor)  ! pseudo-first-order loss rate
-                            r = dt_sub*Lloss(eqn); r2 = r*r
-                            alp(eqn) = (180._wp + 60._wp*r + 11._wp*r2 + r2*r)/(360._wp + 60._wp*r + 12._wp*r2 + r2*r)
-                            Ys(eqn) = y0(eqn) + dt_sub*(prod0(eqn) - loss_i)/(1._wp + alp(eqn)*dt_sub*Lloss(eqn))
-                            if (Ys(eqn) < 0._wp) Ys(eqn) = 0._wp
+                            alp(eqn) = f_qss_alpha(dt_sub*Lloss(eqn))
+                            Ys(eqn) = y0(eqn) + dt_sub*(prod0(eqn) - Lloss(eqn)*y0(eqn))/(1._wp + alp(eqn)*dt_sub*Lloss(eqn))
                         end do
                         ! corrector: re-evaluate rates at the predicted state (T is the Newton guess; T_new is intent(out))
                         call get_temperature(energy, T, Ys, .true., T_new)
@@ -269,12 +267,9 @@ contains
                             prod_p = wr*cdot(eqn)
                             loss_p = wr*ddot(eqn)
                             Lbar = 0.5_wp*(Lloss(eqn) + loss_p/max(Ys(eqn), y_floor))
-                            ! reuse the predictor's alp(eqn) here (and in the denominator) rather than
-                            ! recomputing alpha from the averaged loss Lbar -- a deliberate CHEMEQ2
-                            ! simplification, stability-neutral (alpha in [0.5, 1]) and mitigated by sub-stepping.
+                            alp(eqn) = f_qss_alpha(dt_sub*Lbar)
                             pbar = alp(eqn)*prod_p + (1._wp - alp(eqn))*prod0(eqn)
                             Ys(eqn) = y0(eqn) + dt_sub*(pbar - Lbar*y0(eqn))/(1._wp + alp(eqn)*dt_sub*Lbar)
-                            if (Ys(eqn) < 0._wp) Ys(eqn) = 0._wp
                             Ysum = Ysum + Ys(eqn)
                         end do
                         if (Ysum > y_floor) then
@@ -283,8 +278,7 @@ contains
                                 Ys(eqn) = Ys(eqn)/Ysum
                             end do
                         else
-                            ! Degenerate corrector (every species clipped to zero): fall back to the
-                            ! sub-step's starting composition rather than dividing by a vanishing sum.
+                            ! Degenerate corrector (every species vanished): keep the sub-step's starting composition
                             $:GPU_LOOP(parallelism='[seq]')
                             do eqn = 1, num_species
                                 Ys(eqn) = y0(eqn)
@@ -306,6 +300,18 @@ contains
 
     end subroutine s_chemistry_reaction_substep
 
+    !> alpha-QSS weight for a sub-step with dt*(loss rate) = r. (1 - alpha)*r < 1 for every r >= 0, so a step from non-negative mass
+    !! fractions with non-negative production stays non-negative.
+    pure function f_qss_alpha(r) result(alp)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in) :: r
+        real(wp)             :: alp
+
+        alp = (180._wp + r*(60._wp + r*(11._wp + r)))/(360._wp + r*(60._wp + r*(12._wp + r)))
+
+    end function f_qss_alpha
+
     !> Compute species mass diffusion fluxes at cell interfaces using mixture-averaged diffusivities.
     subroutine s_compute_chemistry_diffusion_flux(idir, q_prim_qp, flux_src_vf, irx, iry, irz, q_T_sf)
 
@@ -420,6 +426,13 @@ contains
                                     & - eqn_idx%species%beg + 1)
                             end do
 
+                            ! The correction flux Y_k*rho_Vic is advective: upwind Y, one side for all species so sum Y = 1
+                            $:GPU_LOOP(parallelism='[seq]')
+                            do i = eqn_idx%species%beg, eqn_idx%species%end
+                                Ys_cell(i - eqn_idx%species%beg + 1) = merge(Ys_L(i - eqn_idx%species%beg + 1), &
+                                        & Ys_R(i - eqn_idx%species%beg + 1), rho_Vic >= 0._wp)
+                            end do
+
                             ! Apply corrections for mass conservation
                             $:GPU_LOOP(parallelism='[seq]')
                             do eqn = eqn_idx%species%beg, eqn_idx%species%end

src/common/m_mpi_common.fpp

@@ -373,6 +373,8 @@ contains
         integer :: ierr  !< Generic flag used to identify and report MPI errors
 
         call MPI_ALLREDUCE(var_loc, var_glb, 1, mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr)
+#else
+        var_glb = var_loc
 #endif
 
     end subroutine s_mpi_allreduce_sum

src/simulation/m_rhs.fpp

@@ -38,6 +38,7 @@ module m_rhs
     use m_igr
     use m_thinc
     use m_pressure_relaxation
+    use m_species_bounds
 
     implicit none
 
@@ -997,6 +998,8 @@ contains
             end if
         end if
 
+        if (chemistry) call s_bound_species_faces(q_prim_qp%vf, qL_rsx_vf, qR_rsx_vf, id)
+
         ! Reconstruct viscous derivatives for viscosity
         if (weno_Re_flux) then
             iv%beg = eqn_idx%mom%beg; iv%end = eqn_idx%mom%end

src/simulation/m_species_bounds.fpp

@@ -0,0 +1,177 @@
+!>
+!! @file
+!! @brief Contains module m_species_bounds
+
+#:include 'macros.fpp'
+
+!> @brief Keeps reconstructed species mass fractions on the simplex (Y >= 0, sum Y = 1). Each cell's face values are first shifted
+!! to sum to 1, then scaled toward the cell average by one theta shared by all species and both faces (Zhang & Shu, JCP 2010). With
+!! face Y on the simplex the HLLC species fluxes are Y_upwind times the mass flux (Larrouturou, JCP 1991), so they sum to it and
+!! keep cell averages non-negative under the Zhang-Shu CFL.
+module m_species_bounds
+
+    use m_derived_types
+    use m_global_parameters
+    use m_mpi_common, only: s_mpi_allreduce_sum
+    use m_ibm, only: ib_markers
+
+    implicit none
+
+    private
+    public :: s_bound_species_faces, s_clean_species, s_report_species_cleanup
+
+    !> Violations below this fraction of rho are roundoff: cleaned but not counted
+    real(wp), parameter :: clean_report_tol = 1.e-12_wp
+    !> Fluid-cell cleanings, and sum of |rho_new - rho_old|, since the last report
+    real(wp) :: clean_cells = 0._wp, clean_mass = 0._wp
+
+contains
+
+    !> Bound the species face values vL/vR of direction id, over the range the reconstruction filled.
+    subroutine s_bound_species_faces(q_prim_vf, vL, vR, id)
+
+        type(scalar_field), dimension(sys_size), intent(in)                                    :: q_prim_vf
+        real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: vL, vR
+        integer, intent(in)                                                                    :: id
+        type(int_bounds_info), dimension(3)                                                    :: b
+        real(wp)                                                                               :: w, theta, sL, sR, ybar, yL, yR, xi
+        integer                                                                                :: i, j, k, l, polyn
+
+        ! Same range as s_reconstruct_cell_boundary_values (m_rhs.fpp)
+
+        polyn = merge(weno_polyn, muscl_polyn, recon_type == recon_type_weno)
+        b = idwbuff
+        b(id)%beg = b(id)%beg + polyn; b(id)%end = b(id)%end - polyn
+        w = merge(f_lobatto_weight(polyn), 0.5_wp, recon_type == recon_type_weno)  ! MUSCL is linear: no interior point
+
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, theta, sL, sR, ybar, yL, yR, xi]', copyin='[b, w]')
+        do l = b(3)%beg, b(3)%end
+            do k = b(2)%beg, b(2)%end
+                do j = b(1)%beg, b(1)%end
+                    sL = 0._wp; sR = 0._wp
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do i = eqn_idx%species%beg, eqn_idx%species%end
+                        sL = sL + vL(j, k, l, i); sR = sR + vR(j, k, l, i)
+                    end do
+                    theta = 1._wp
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do i = eqn_idx%species%beg, eqn_idx%species%end
+                        ybar = q_prim_vf(i)%sf(j, k, l)
+                        yL = vL(j, k, l, i) + ybar*(1._wp - sL)
+                        yR = vR(j, k, l, i) + ybar*(1._wp - sR)
+                        vL(j, k, l, i) = yL; vR(j, k, l, i) = yR
+                        theta = min(theta, f_theta(ybar, yL), f_theta(ybar, yR))
+                        if (w < 0.5_wp) then
+                            xi = (ybar - w*(yL + yR))/(1._wp - 2._wp*w)
+                            theta = min(theta, f_theta(ybar, xi))
+                        end if
+                    end do
+                    if (theta < 1._wp) then
+                        $:GPU_LOOP(parallelism='[seq]')
+                        do i = eqn_idx%species%beg, eqn_idx%species%end
+                            ybar = q_prim_vf(i)%sf(j, k, l)
+                            vL(j, k, l, i) = ybar + theta*(vL(j, k, l, i) - ybar)
+                            vR(j, k, l, i) = ybar + theta*(vR(j, k, l, i) - ybar)
+                        end do
+                    end if
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+    end subroutine s_bound_species_faces
+
+    !> Fallback (as PeleC's clean_massfrac): clip each rho*Y_k to [0, rho], set rho = sum rho*Y_k, and scale momentum and energy by
+    !! rho_new/rho_old so velocity and specific energy are kept. A cell whose clipped sum is not positive is left untouched. Every
+    !! cell is cleaned, but only fluid cells (not IB cells) off by > clean_report_tol*rho, or left untouched, are counted as cell
+    !! cleanings (per cell per RK stage) for s_report_species_cleanup.
+    subroutine s_clean_species(q_cons_vf)
+
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
+        real(wp)                                               :: rho_old, rho_new, f, n_cells, d_mass, viol, rhoY
+        logical                                                :: is_fluid
+        integer                                                :: i, j, k, l
+
+        n_cells = 0._wp; d_mass = 0._wp
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, rho_old, rho_new, f, viol, rhoY, is_fluid]', reduction='[[n_cells, &
+                            & d_mass]]', reductionOp='[+]')
+        do l = 0, p
+            do k = 0, n
+                do j = 0, m
+                    rho_old = real(q_cons_vf(eqn_idx%cont%beg)%sf(j, k, l), wp)
+                    viol = 0._wp
+                    rho_new = 0._wp
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do i = eqn_idx%species%beg, eqn_idx%species%end
+                        rhoY = real(q_cons_vf(i)%sf(j, k, l), wp)
+                        viol = max(viol, -rhoY, rhoY - rho_old)
+                        rho_new = rho_new + min(max(rhoY, 0._wp), rho_old)
+                    end do
+                    if (viol > 0._wp .and. rho_old > 0._wp) then
+                        is_fluid = .true.
+                        if (ib) is_fluid = (ib_markers%sf(j, k, l) == 0)
+                        if (rho_new > 0._wp) then
+                            $:GPU_LOOP(parallelism='[seq]')
+                            do i = eqn_idx%species%beg, eqn_idx%species%end
+                                rhoY = min(max(real(q_cons_vf(i)%sf(j, k, l), wp), 0._wp), rho_old)
+                                q_cons_vf(i)%sf(j, k, l) = real(rhoY, stp)
+                            end do
+                            f = rho_new/rho_old
+                            q_cons_vf(eqn_idx%cont%beg)%sf(j, k, l) = real(rho_new, stp)
+                            $:GPU_LOOP(parallelism='[seq]')
+                            do i = eqn_idx%mom%beg, eqn_idx%E
+                                q_cons_vf(i)%sf(j, k, l) = real(f*real(q_cons_vf(i)%sf(j, k, l), wp), stp)
+                            end do
+                            if (is_fluid) then
+                                if (viol > clean_report_tol*rho_old) n_cells = n_cells + 1._wp
+                                d_mass = d_mass + abs(rho_new - rho_old)
+                            end if
+                        else if (is_fluid) then
+                            n_cells = n_cells + 1._wp
+                        end if
+                    end if
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+        clean_cells = clean_cells + n_cells; clean_mass = clean_mass + d_mass
+
+    end subroutine s_clean_species
+
+    !> Print, once per time step and only if nonzero, how many cell cleanings s_clean_species made across all ranks.
+    impure subroutine s_report_species_cleanup(t_step)
+
+        integer, intent(in) :: t_step
+        real(wp)            :: cells_glb, mass_glb
+
+        call s_mpi_allreduce_sum(clean_cells, cells_glb)
+        call s_mpi_allreduce_sum(clean_mass, mass_glb)
+        if (proc_rank == 0 .and. cells_glb > 0._wp) print '(A,I0,A,I0,A,ES10.3)', 'Species cleanup at step ', t_step, ': ', &
+            & nint(cells_glb), ' cell cleanings, sum |d rho| = ', mass_glb
+        clean_cells = 0._wp; clean_mass = 0._wp
+
+    end subroutine s_report_species_cleanup
+
+    !> Largest theta in [0, 1] keeping ybar + theta*(y - ybar) >= 0, given ybar >= 0.
+    pure function f_theta(ybar, y) result(theta)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in) :: ybar, y
+        real(wp)             :: theta
+
+        theta = 1._wp
+        if (y < 0._wp) theta = max(0._wp, ybar)/max(ybar - y, sgm_eps)
+
+    end function f_theta
+
+    !> Endpoint weight of the Gauss-Lobatto rule exact for the reconstruction's degree 2*polyn: 1/(N(N-1)), N = polyn + 2.
+    pure function f_lobatto_weight(polyn) result(w)
+
+        integer, intent(in) :: polyn
+        real(wp)            :: w
+
+        w = 1._wp/real((polyn + 2)*(polyn + 1), wp)
+
+    end function f_lobatto_weight
+
+end module m_species_bounds

src/simulation/m_time_steppers.fpp

@@ -12,6 +12,7 @@ module m_time_steppers
     use m_global_parameters
     use m_rhs
     use m_chemistry
+    use m_species_bounds
     use m_reactive_burn, only: s_reactive_burn_substep
     use m_pressure_relaxation
     use m_hypoelastic, only: s_enforce_cont_damage_bounds
@@ -572,6 +573,7 @@ contains
             end if
 
             if (cont_damage) call s_enforce_cont_damage_bounds(q_cons_ts(1)%vf)
+            if (chemistry) call s_clean_species(q_cons_ts(1)%vf)
 
             ! Grind: minimum wall-clock time of a full RK stage (compute + halo H2D/D2H +
             ! update + IBM correction, aside from I/O) over steady-state stages. Wall clock
@@ -601,6 +603,8 @@ contains
             call nvtxEndRange
         end if
 
+        if (chemistry) call s_report_species_cleanup(t_step)
+
         ! Operator-split condensed-phase burn: integrate the progress variable per cell after the flow
         ! update, with sub-stepping, instead of adding the source to the flow RHS (rburn%substeps > 0).
         if (reactive_burn .and. rburn%substeps > 0) then