Merge master into ci-phoenix-bench-single-node

Conflict in submit-slurm-job.sh: master (#1832) split partition selection by job_type, keeping the single-partition bench selector behind BENCH_GPU_PARTITION. This branch removes that selector entirely (bench-pair.sh benches master and PR in one job on the same GPUs), so the job_type branch is dead code here. Took this branch's unconditional partition list; select-gpu-partition.sh stays deleted.

Claude-Session: https://claude.ai/code/session_01G77jhrA4JPDz5TqJzt8ACC

Author: sbryngelson

src/common/m_derived_types.fpp

@@ -526,11 +526,12 @@ module m_derived_types
     !> Condensed-phase reactive-burn (programmed pressure detonation) parameters. The rate is
     !> dlambda/dt = k (1 - lambda) ((p - pign)/pref)^n, optionally scaled by exp(-ta/T) when ta > 0.
     type reactive_burn_parameters
-        real(wp) :: k     !< Rate coefficient [1/s]
-        real(wp) :: pign  !< Ignition pressure threshold [Pa]
-        real(wp) :: pref  !< Reference pressure for the pressure drive [Pa]
-        real(wp) :: n     !< Pressure-drive exponent
-        real(wp) :: ta    !< Activation temperature [K] (0 = pure pressure-driven; > 0 adds exp(-ta/T))
+        real(wp) :: k         !< Rate coefficient [1/s]
+        real(wp) :: pign      !< Ignition pressure threshold [Pa]
+        real(wp) :: pref      !< Reference pressure for the pressure drive [Pa]
+        real(wp) :: n         !< Pressure-drive exponent
+        real(wp) :: ta        !< Activation temperature [K] (0 = pure pressure-driven; > 0 adds exp(-ta/T))
+        integer  :: substeps  !< Operator-split sub-steps per time step (0 = source added to the flow RHS)
     end type reactive_burn_parameters
 
     !> Lagrangian bubble parameters

src/common/m_global_parameters_common.fpp

@@ -339,6 +339,7 @@ contains
         rburn%pref = dflt_real
         rburn%n = dflt_real
         rburn%ta = 0._wp
+        rburn%substeps = 0
 
         ! Case-optimization params: under case-opt these are compile-time constants in sim (skip assignment); in pre/post
         ! MFC_CASE_OPTIMIZATION is always False so the block always executes there.

src/simulation/m_reactive_burn.fpp

@@ -15,14 +15,45 @@
 module m_reactive_burn
 
     use m_global_parameters
-    use m_variables_conversion, only: f_sg_thermal
+    use m_variables_conversion, only: f_sg_thermal, s_compute_mixture_coefficients, f_pressure
 
     implicit none
 
-    private; public :: s_compute_reactive_burn
+    private; public :: s_compute_reactive_burn, s_reactive_burn_substep
 
 contains
 
+    !> Programmed-burn rate dlambda/dt for one cell state. Both the RHS source and the operator-split integrator call this, so the
+    !! rate law is stated once.
+    !! @param pres             Mixture pressure
+    !! @param lambda           Reaction progress, i.e. the product volume fraction
+    !! @param alpha_rho_react  Reactant partial density, for the optional Arrhenius factor
+    !! @param alpha_react      Reactant volume fraction, for the optional Arrhenius factor
+    !! @param rate             dlambda/dt; zero below the ignition pressure and once the reactant is spent
+    subroutine s_burn_rate(pres, lambda, alpha_rho_react, alpha_react, rate)
+
+        $:GPU_ROUTINE(function_name='s_burn_rate', parallelism='[seq]', cray_inline=True)
+
+        real(wp), intent(in)  :: pres, lambda, alpha_rho_react, alpha_react
+        real(wp), intent(out) :: rate
+        real(wp)              :: drive
+
+        ! Pressure-driven programmed burn: fires only behind the shock (p > rburn%pign).
+        drive = (pres - rburn%pign)/rburn%pref
+        if (drive > 0._wp .and. lambda < 1._wp) then
+            rate = rburn%k*(1._wp - lambda)*drive**rburn%n
+            ! Optional Arrhenius dependence on the reactant phasic temperature from the stiffened-gas
+            ! EOS. rburn%ta = 0, the default, leaves the pure pressure-driven rate unchanged.
+            if (rburn%ta > 0._wp) then
+                rate = rate*exp(-rburn%ta/f_sg_thermal(pres, alpha_rho_react/max(alpha_react, sgm_eps), isentrope_n(1), &
+                                & isentrope_B(1), cvs(1)))
+            end if
+        else
+            rate = 0._wp
+        end if
+
+    end subroutine s_burn_rate
+
     !> Add the programmed-burn reaction source to the continuity and volume-fraction RHS.
     !! @param rhs_vf     Right-hand-side accumulator (inout)
     !! @param q_cons_vf  Conserved variables (partial densities live here)
@@ -34,9 +65,9 @@ contains
         type(scalar_field), dimension(sys_size), intent(in)    :: q_cons_vf, q_prim_vf
         type(int_bounds_info), dimension(1:3), intent(in)      :: bounds
         integer                                                :: x, y, z
-        real(wp)                                               :: rho, pres, lambda, rate, mdot, drive, T_r
+        real(wp)                                               :: rho, pres, lambda, rate, mdot
 
-        $:GPU_PARALLEL_LOOP(collapse=3, private='[rho, pres, lambda, rate, mdot, drive, T_r]', copyin='[bounds]')
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[rho, pres, lambda, rate, mdot]', copyin='[bounds]')
         do z = bounds(3)%beg, bounds(3)%end
             do y = bounds(2)%beg, bounds(2)%end
                 do x = bounds(1)%beg, bounds(1)%end
@@ -45,20 +76,9 @@ contains
                     pres = q_prim_vf(eqn_idx%E)%sf(x, y, z)
                     lambda = q_prim_vf(eqn_idx%adv%beg + 1)%sf(x, y, z)  ! reaction progress = product volume fraction
 
-                    ! pressure-driven programmed burn: fires only behind the shock (p > rburn%pign)
-                    drive = (pres - rburn%pign)/rburn%pref
-                    if (drive > 0._wp .and. lambda < 1._wp) then
-                        rate = rburn%k*(1._wp - lambda)*drive**rburn%n  ! dlambda/dt
-
-                        ! Optional Arrhenius temperature dependence: rate *= exp(-rburn%ta/T_r), with T_r the
-                        ! reactant phasic temperature from the stiffened-gas EOS T = (p + pi_inf)/((Gamma-1) rho cv).
-                        ! rburn%ta = 0 (default) leaves the pure pressure-driven rate unchanged.
-                        if (rburn%ta > 0._wp) then
-                            T_r = f_sg_thermal(pres, q_cons_vf(eqn_idx%cont%beg)%sf(x, y, z)/q_prim_vf(eqn_idx%adv%beg)%sf(x, y, &
-                                               & z), isentrope_n(1), isentrope_B(1), cvs(1))
-                            rate = rate*exp(-rburn%ta/T_r)
-                        end if
-
+                    call s_burn_rate(pres, lambda, q_cons_vf(eqn_idx%cont%beg)%sf(x, y, z), q_prim_vf(eqn_idx%adv%beg)%sf(x, y, &
+                                     & z), rate)
+                    if (rate > 0._wp) then
                         mdot = rho*rate  ! mass reactant -> product
 
                         ! continuity: reactant loses mass, product gains it
@@ -76,4 +96,83 @@ contains
 
     end subroutine s_compute_reactive_burn
 
+    !> Operator-split alternative to s_compute_reactive_burn, used when rburn%substeps > 0. The flow is frozen and the burn ODE is
+    !! integrated over one time step in equal sub-steps, so the reaction time scale is decoupled from the acoustic CFL. The mixture
+    !! pressure is re-evaluated from the frozen internal energy each sub-step, which is what carries the rate's own feedback: the
+    !! coefficients move as the reactant becomes product.
+    !! @param q_cons_vf  Conserved variables, updated in place
+    !! @param dtime      Time step to integrate across
+    !! @param bounds     Interior cell bounds
+    subroutine s_reactive_burn_substep(q_cons_vf, dtime, bounds)
+
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
+        real(wp), intent(in)                                   :: dtime
+        type(int_bounds_info), dimension(1:3), intent(in)      :: bounds
+        integer                                                :: x, y, z, i, sub
+        real(wp)                                               :: rho, pres, lambda, rate
+        real(wp)                                               :: dt_sub, e_int, gamma_mix, pi_inf_mix, qv_mix
+        real(wp)                                               :: rho_mix, dlambda, dmass
+
+        ! Bounded by num_fluids_max, not num_fluids: under case optimization num_fluids is a compile-time
+        ! constant that can be 1, and the reactant/product indices below are literal. The validator holds
+        ! reactive_burn to two fluids, so only the first two entries are ever used.
+        real(wp), dimension(num_fluids_max) :: alpha_rho, alpha
+
+        dt_sub = dtime/real(rburn%substeps, wp)
+
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_rho, alpha, rho, pres, lambda, rate, e_int, gamma_mix, pi_inf_mix, &
+                            & qv_mix, rho_mix, dlambda, dmass, i, sub]', copyin='[bounds, dt_sub]')
+        do z = bounds(3)%beg, bounds(3)%end
+            do y = bounds(2)%beg, bounds(2)%end
+                do x = bounds(1)%beg, bounds(1)%end
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do i = 1, num_fluids
+                        alpha_rho(i) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(x, y, z)
+                        alpha(i) = q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(x, y, z)
+                    end do
+                    rho = alpha_rho(1) + alpha_rho(2)
+
+                    ! Internal energy per unit volume is what the burn conserves: it moves mass and volume
+                    ! between the phases, and the qv difference surfaces as pressure through the mixture EOS.
+                    e_int = q_cons_vf(eqn_idx%E)%sf(x, y, z)
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do i = eqn_idx%mom%beg, eqn_idx%mom%end
+                        e_int = e_int - 0.5_wp*q_cons_vf(i)%sf(x, y, z)**2/rho
+                    end do
+
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do sub = 1, rburn%substeps
+                        lambda = alpha(2)
+                        call s_compute_mixture_coefficients(alpha_rho, alpha, rho_mix, gamma_mix, pi_inf_mix, qv_mix)
+                        pres = f_pressure(e_int, gamma_mix, pi_inf_mix, qv_mix)
+                        call s_burn_rate(pres, lambda, alpha_rho(1), alpha(1), rate)
+                        if (rate <= 0._wp) exit
+                        ! A sub-step longer than the reaction time would carry the progress variable past one.
+                        ! Stop it there and hand over the reactant's remaining mass in the same sub-step: capping
+                        ! the two independently strands mass at zero volume, and the EOS divides one by the other.
+                        if (rate*dt_sub >= 1._wp - lambda) then
+                            dlambda = 1._wp - lambda
+                            dmass = alpha_rho(1)
+                        else
+                            dlambda = rate*dt_sub
+                            dmass = min(rho*dlambda, alpha_rho(1))
+                        end if
+                        alpha(1) = alpha(1) - dlambda
+                        alpha(2) = alpha(2) + dlambda
+                        alpha_rho(1) = alpha_rho(1) - dmass
+                        alpha_rho(2) = alpha_rho(2) + dmass
+                    end do
+
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do i = 1, num_fluids
+                        q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(x, y, z) = alpha_rho(i)
+                        q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(x, y, z) = alpha(i)
+                    end do
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+    end subroutine s_reactive_burn_substep
+
 end module m_reactive_burn

src/simulation/m_rhs.fpp

@@ -869,7 +869,9 @@ contains
             call nvtxEndRange
         end if
 
-        if (reactive_burn) then
+        ! With rburn%substeps > 0 the burn is integrated by operator splitting after the flow
+        ! update (s_reactive_burn_substep), not added to the flow RHS here.
+        if (reactive_burn .and. rburn%substeps == 0) then
             call nvtxStartRange("RHS-REACTIVE-BURN")
             call s_compute_reactive_burn(rhs_vf, q_cons_qp%vf, q_prim_qp%vf, idwint)
             call nvtxEndRange

src/simulation/m_riemann_solver_hllc.fpp

@@ -839,17 +839,41 @@ contains
                         ! statement and private variable from the pure-fluid emission, keeping its body and directive
                         ! identical to the single kernel a build without hypoelasticity would compile. Sharing one kernel
                         ! pinned it at the GPU register ceiling for every HLLC user.
+                        ! One source of truth for this kernel's private variables: both emissions of the shared body take
+                        ! _hllc_s*, and only the hypoelastic one adds _hllc_e*. Two hand-written lists drifted apart once --
+                        ! c_sum_Yi_Phi was private in one and shared in the other, which races under OpenMP offload.
+                        ! Names are lists joined once, so no fragment carries a trailing separator to get wrong.
+                        #:set _hllc_s1 = ['i', 'j', 'k', 'l', 'q', 'T_L', 'T_R', 'vel_L_rms', 'vel_R_rms', 'pres_L', 'pres_R', &
+                            & 'rho_L', 'gamma_L', 'pi_inf_L', 'qv_L', 'rho_R', 'gamma_R']
+                        #:set _hllc_s2 = ['pi_inf_R', 'qv_R', 'alpha_L_sum', 'alpha_R_sum', 'E_L', 'E_R', 'MW_L', 'MW_R', &
+                            & 'R_gas_L', 'R_gas_R', 'Cp_L', 'Cp_R', 'Cv_L', 'Cv_R', 'c_sum_Yi_Phi']
+                        #:set _hllc_s3 = ['Gamm_L', 'Gamm_R', 'Y_L', 'Y_R', 'H_L', 'H_R', 'qv_avg', 'rho_avg', 'gamma_avg', &
+                            & 'H_avg', 'c_L', 'c_R', 'c_avg', 's_P', 's_M', 'xi_P', 'xi_M', 'xi_L']
+                        #:set _hllc_s4 = ['xi_R', 'xi_L_m1', 'xi_R_m1', 'Ms_L', 'Ms_R', 'pres_SL', 'pres_SR', 'vel_L', 'vel_R', &
+                            & 'Re_L', 'Re_R', 'alpha_L', 'alpha_R', 'alpha_rho_L', 'alpha_rho_R']
+                        #:set _hllc_s5 = ['alpha_lim_L', 'alpha_lim_R', 's_L', 's_R', 's_S', 'vel_avg_rms', 'pcorr', 'Ys_L', &
+                            & 'Ys_R', 'Xs_L', 'Xs_R', 'Gamma_iL', 'Gamma_iR', 'Cp_iL', 'Cp_iR']
+                        #:set _hllc_s6 = ['R_species', 'h_iL', 'h_iR']
+                        #:set _hllc_e1 = ['ptilde_L', 'ptilde_R', 'tau_e_L', 'tau_e_R', 'G_L', 'G_R', 'damage_L', 'damage_R', &
+                            & 'U_L', 'U_R', 'F_L', 'F_R', 'F_star_L', 'F_star_R', 'F_HLLC']
+                        #:set _hllc_e2 = ['u_n_HLLC', 'u_t_HLLC', 'u_t2_HLLC', 'pres_tot_L', 'pres_tot_R', 'u_n_L', 'u_n_R', &
+                            & 'u_t_L', 'u_t_R', 'u_t2_L', 'u_t2_R', 'tau_nn_L', 'tau_nn_R']
+                        #:set _hllc_e3 = ['tau_nt_L', 'tau_nt_R', 'tau_tt_L', 'tau_tt_R', 'tau_nt2_L', 'tau_nt2_R', 'tau_t2t2_L', &
+                            & 'tau_t2t2_R', 'tau_t1t2_L', 'tau_t1t2_R', 'tau_qq_L', 'tau_qq_R']
+                        #:set _hllc_e4 = ['p_face', 'tau_qq_face', 'A_L', 'A_R', 'denom_A', 'u_t_star', 'tau_nt_star', &
+                            & 'u_t2_star', 'tau_nt2_star', 'pres_tot_star', 'F_HLL', 'u_n_HLL_trace']
+                        #:set _hllc_e5 = ['u_t_HLL_trace', 'u_t2_HLL_trace', 'p_face_HLL', 'tau_qq_face_HLL', 'tau_nn_HLL', &
+                            & 'phi', 'Sigma_L', 'Sigma_R', 'dSigma', 'Sigma_ref', 'a_L_ref']
+                        #:set _hllc_e6 = ['a_R_ref', 'a_ref', 'du_t', 'dtau_nt', 'du_t2', 'dtau_nt2', 'sensor_ptot', 'sensor_vt', &
+                            & 'sensor_tnt', 'sensor_combined', 'idx_phys']
                         #:if HYPO
-                            ! Private list split across _hllc_p1/p2/p3 for Fypp line-length limits
-                            #:set _hllc_p1 = '[i, j, k, l, q, T_L, T_R, vel_L_rms, vel_R_rms, pres_L, pres_R, rho_L, gamma_L, pi_inf_L, qv_L, rho_R, gamma_R, pi_inf_R, qv_R, alpha_L_sum, alpha_R_sum, E_L, E_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, c_sum_Yi_Phi, Gamm_L, Gamm_R, Y_L, Y_R, H_L, H_R, qv_avg, rho_avg, gamma_avg, H_avg, c_L, c_R, c_avg, s_P, s_M, xi_P, xi_M, xi_L, xi_R, xi_L_m1, xi_R_m1, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, alpha_rho_L, alpha_rho_R, alpha_lim_L, alpha_lim_R, s_L, s_R, s_S, vel_avg_rms, pcorr, ptilde_L, ptilde_R, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, R_species, h_iL, h_iR, tau_e_L, tau_e_R, G_L, G_R, damage_L, damage_R,'
-                            #:set _hllc_p2 = 'U_L, U_R, F_L, F_R, F_star_L, F_star_R, F_HLLC, u_n_HLLC, u_t_HLLC, u_t2_HLLC, pres_tot_L, pres_tot_R, u_n_L, u_n_R, u_t_L, u_t_R, u_t2_L, u_t2_R, tau_nn_L, tau_nn_R, tau_nt_L, tau_nt_R, tau_tt_L, tau_tt_R, tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R, tau_qq_L, tau_qq_R, p_face, tau_qq_face, A_L, A_R, denom_A, u_t_star, tau_nt_star, u_t2_star, tau_nt2_star, pres_tot_star,'
-                            #:set _hllc_p3 = 'F_HLL, u_n_HLL_trace, u_t_HLL_trace, u_t2_HLL_trace, p_face_HLL, tau_qq_face_HLL, tau_nn_HLL, phi, Sigma_L, Sigma_R, dSigma, Sigma_ref, a_L_ref, a_R_ref, a_ref, du_t, dtau_nt, du_t2, dtau_nt2, sensor_ptot, sensor_vt, sensor_tnt, sensor_combined, idx_phys]'
-                            #:set _hllc_priv = _hllc_p1 + _hllc_p2 + _hllc_p3
+                            #:set _hllc_priv = '[' + ', '.join(_hllc_s1 + _hllc_s2 + _hllc_s3 + _hllc_s4 + _hllc_s5 + _hllc_s6 &
+                                                               & + _hllc_e1 + _hllc_e2 + _hllc_e3 + _hllc_e4 + _hllc_e5 &
+                                                               & + _hllc_e6) + ']'
                         #:else
-                            ! Master's pure-fluid private list, unchanged
-                            #:set _hllc_priv = '[i, T_L, T_R, vel_L_rms, vel_R_rms, pres_L, pres_R, rho_L, gamma_L, pi_inf_L, qv_L, rho_R, gamma_R, pi_inf_R, qv_R, alpha_L_sum, alpha_R_sum, E_L, E_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, Y_L, Y_R, H_L, H_R, qv_avg, rho_avg, gamma_avg, H_avg, c_L, c_R, c_avg, s_P, s_M, xi_P, xi_M, xi_L, xi_R, xi_L_m1, xi_R_m1, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, alpha_rho_L, alpha_rho_R, alpha_lim_L, alpha_lim_R, s_L, s_R, s_S, vel_avg_rms, pcorr, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, R_species, h_iL, h_iR]'
+                            #:set _hllc_priv = '[' + ', '.join(_hllc_s1 + _hllc_s2 + _hllc_s3 + _hllc_s4 + _hllc_s5 + _hllc_s6) &
+                                                               & + ']'
                         #:endif
-                        ! The two calls below are identical on purpose. An offload kernel is named
                         ! after the .fpp line of its GPU_PARALLEL_LOOP, so one shared call would give
                         ! both emissions the same name; amdflang then launches the wrong one and a
                         ! hypoelastic run faults inside the pure-fluid kernel. Two call sites are what

src/simulation/m_riemann_solver_lf.fpp

@@ -94,8 +94,8 @@ contains
                                     & Re_L, Re_R, s_L, s_R, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, pcorr, &
                                     & vel_grad_L, vel_grad_R, idx_right_phys, vel_L_rms, vel_R_rms, alpha_L_sum, alpha_R_sum, &
                                     & pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, &
-                                    & E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, &
-                                    & MW_R, T_L, T_R, Y_L, Y_R]', firstprivate='[Re_size_loc1, Re_size_loc2]')
+                                    & Gamm_L, Gamm_R, E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, &
+                                    & R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, Y_R]', firstprivate='[Re_size_loc1, Re_size_loc2]')
                 do l = ${Z_BND}$%beg, ${Z_BND}$%end
                     do k = ${Y_BND}$%beg, ${Y_BND}$%end
                         do j = ${X_BND}$%beg, ${X_BND}$%end