Privatize c_sum_Yi_Phi in the pure-fluid HLLC kernel, and evaluate the JWL and Vinet isentrope in closed form

Author: sbryngelson

src/common/m_variables_conversion.fpp

@@ -1421,6 +1421,20 @@ contains
 
     end function f_is_state_dependent
 
+    !> True when fluid i's reference curve is itself an isentrope (de_ref = -p_ref d(1/rho), which holds for JWL and Vinet but not
+    !! for the Mie-Gruneisen Hugoniot) and its Gruneisen coefficient is constant. Those two together make the isentrope through any
+    !! state closed-form, so it never has to be integrated.
+    function f_has_isentropic_reference(i) result(yes)
+
+        $:GPU_ROUTINE(function_name='f_has_isentropic_reference', parallelism='[seq]', cray_inline=True)
+
+        integer, intent(in) :: i
+        logical             :: yes
+
+        yes = (eoss(i) == eos_jwl .or. eoss(i) == eos_vinet) .and. gruneisen_as(i) == 0._wp
+
+    end function f_has_isentropic_reference
+
     !> Gamma, Pi, dPi/drho and dGamma/drho of fluid i at density rho, the coefficients of rho e = Gamma p + Pi(rho). Stiffened and
     !! ideal gas keep the constants resolved at init, bit for bit.
     subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
@@ -1588,11 +1602,19 @@ contains
         real(wp), intent(in)  :: pres, rho, xi
         integer, intent(in)   :: i
         real(wp), intent(out) :: p_isen
+        real(wp)              :: p_ref_from, p_ref_to, e_ref, dp_drho, de_drho, G0, dG0
 
-        if (f_is_state_dependent(i)) then
-            call s_rk4(ode_isentrope, i, rho, pres, xi*rho, p_isen)
-        else
+        if (.not. f_is_state_dependent(i)) then
             p_isen = (pres + isentrope_B(i))*xi**isentrope_n(i) - isentrope_B(i)
+        else if (f_has_isentropic_reference(i)) then
+            ! Exact: the offset from an isentropic reference obeys dDelta/Delta = Gamma drho/rho, so
+            ! p - p_ref scales as (rho'/rho)**(1 + Gamma). Integrating it instead costs a decimal per
+            ! doubling of the expansion and turns the pressure negative past roughly twentyfold.
+            call s_reference_curve(rho, i, p_ref_from, e_ref, dp_drho, de_drho, G0, dG0)
+            call s_reference_curve(xi*rho, i, p_ref_to, e_ref, dp_drho, de_drho, G0, dG0)
+            p_isen = p_ref_to + (pres - p_ref_from)*xi**(1._wp + gruneisen0s(i))
+        else
+            call s_rk4(ode_isentrope, i, rho, pres, xi*rho, p_isen)
         end if
 
     end subroutine s_phase_pressure_on_isentrope

src/simulation/m_pressure_relaxation.fpp

@@ -194,6 +194,11 @@ contains
             end if
         end do
 
+        ! A solve that ran out of iterations leaves rho_K_s holding whatever the last iterate produced, and writing
+        ! that into the volume fractions is how one bad cell becomes a NaN a few steps later. Leave the cell
+        ! unrelaxed instead. Written as .not. (<=) so a NaN residual takes the same path as a diverged one.
+        if (.not. (abs(f_pres) <= TOLERANCE)) return
+
         ! Update volume fractions
         $:GPU_LOOP(parallelism='[seq]')
         do i = 1, num_fluids

src/simulation/m_riemann_solver_hllc.fpp

@@ -857,8 +857,10 @@ contains
                             #: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
                         #: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]'
+                            ! c_sum_Yi_Phi is written by the chemistry block inside this kernel, so it has to be
+                            ! private here as it already is in the hypoelastic list above. Shared, every thread wrote
+                            ! the same scalar and the reacting Roe sound speed came out different run to run.
+                            #: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, 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, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, R_species, h_iL, h_iR]'
                         #: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