Pass an eos_state to s_compute_speed_of_sound

src/common/m_derived_types.fpp

@@ -121,6 +121,29 @@ module m_derived_types
         type(int_bounds_info) :: x, y, z
     end type bc_xyz_info
 
+    !> Thermodynamic state handed to the equation-of-state routines.
+    !>
+    !> Carries scalars only. The volume-fraction array `adv` stays a separate argument: a derived-type
+    !> component cannot have a runtime extent, and `num_fluids` is a parameter only under case
+    !> optimization, so `dimension(num_fluids)` will not compile in a general build. Padding to
+    !> num_fluids_max would compile but put ten reals into a per-cell private struct on device.
+    !>
+    !> `H` is the specific total enthalpy and must include `qv`, because the sound-speed relation
+    !> subtracts `qv/rho`. Build states with s_eos_state so that invariant holds by construction
+    !> rather than by convention; s_eos_state_roe exists for the Roe-averaged paths, which supply an
+    !> `H` that is deliberately not the exact state enthalpy.
+    !> Contains no allocatable members - safe to use inside device routines.
+    type eos_state
+        real(wp) :: rho      !< Mixture density
+        real(wp) :: pres     !< Pressure
+        real(wp) :: gamma    !< Stiffened-gas gamma (1/(Gamma-1))
+        real(wp) :: pi_inf   !< Stiffened-gas stiffness
+        real(wp) :: qv       !< Heat of formation (volumetric)
+        real(wp) :: vel_sum  !< |u|^2
+        real(wp) :: H        !< Specific total enthalpy, including qv
+        real(wp) :: c_c      !< Roe-averaged chemistry sound-speed term (0 when unused)
+    end type eos_state
+
     !> QBMM moment index mappings - separate from bub beg/end so eqn_idx contains no allocatables.
     type qbmm_idx_info
         integer, dimension(:), allocatable     :: rs       !< R moment indices per bubble bin

src/common/m_variables_conversion.fpp

@@ -25,7 +25,7 @@ module m_variables_conversion
         & s_convert_mixture_to_mixture_variables, s_convert_species_to_mixture_variables, &
         & s_convert_species_to_mixture_variables_kernel, s_convert_conservative_to_primitive_variables, &
         & s_convert_primitive_to_conservative_variables, s_convert_primitive_to_flux_variables, s_compute_pressure, &
-        & s_compute_species_fraction, s_compute_speed_of_sound, s_compute_fast_magnetosonic_speed, &
+        & s_compute_species_fraction, s_compute_speed_of_sound, s_compute_fast_magnetosonic_speed, s_eos_state, s_eos_state_roe, &
         & s_finalize_variables_conversion_module, gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
 
     real(wp), allocatable, dimension(:)   :: Gs_vc
@@ -1202,55 +1202,92 @@ contains
 
     end subroutine s_finalize_variables_conversion_module
 
+    !> Build an exact thermodynamic state. The specific total enthalpy is derived here rather than supplied, so it cannot disagree
+    !! with qv: H = (E + p)/rho with E = gamma*p + pi_inf + qv + rho*|u|^2/2. Every caller that previously open-coded the closed
+    !! form should use this, which makes the defect in #1707 unrepresentable.
+    subroutine s_eos_state(s, pres, rho, gamma, pi_inf, qv, vel_sum)
+
+        $:GPU_ROUTINE(function_name='s_eos_state', parallelism='[seq]', cray_inline=True)
+
+        type(eos_state), intent(out) :: s
+        real(wp), intent(in)         :: pres, rho, gamma, pi_inf, qv, vel_sum
+
+        s%pres = pres
+        s%rho = rho
+        s%gamma = gamma
+        s%pi_inf = pi_inf
+        s%qv = qv
+        s%vel_sum = vel_sum
+        s%c_c = 0._wp
+        s%H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho + 5.e-1_wp*vel_sum
+
+    end subroutine s_eos_state
+
+    !> Build a state whose enthalpy is supplied by the caller. The Roe-averaged Riemann paths, the chemistry Roe branch and the
+    !! relativistic branch all pass an H that is deliberately not the exact enthalpy of the state, so they cannot use s_eos_state.
+    subroutine s_eos_state_roe(s, pres, rho, gamma, pi_inf, qv, vel_sum, H)
+
+        $:GPU_ROUTINE(function_name='s_eos_state_roe', parallelism='[seq]', cray_inline=True)
+
+        type(eos_state), intent(out) :: s
+        real(wp), intent(in)         :: pres, rho, gamma, pi_inf, qv, vel_sum, H
+
+        s%pres = pres
+        s%rho = rho
+        s%gamma = gamma
+        s%pi_inf = pi_inf
+        s%qv = qv
+        s%vel_sum = vel_sum
+        s%H = H
+        s%c_c = 0._wp
+
+    end subroutine s_eos_state_roe
+
     !> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models.
-    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv)
+    subroutine s_compute_speed_of_sound(s, adv, c)
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
-        real(wp), intent(in) :: pres
-        real(wp), intent(in) :: rho, gamma, pi_inf, qv
-        real(wp), intent(in) :: H
+        type(eos_state), intent(in) :: s
         #:if not MFC_CASE_OPTIMIZATION and USING_AMD
             real(wp), dimension(3), intent(in) :: adv
         #:else
             real(wp), dimension(num_fluids), intent(in) :: adv
         #:endif
-        real(wp), intent(in)  :: vel_sum
-        real(wp), intent(in)  :: c_c
         real(wp), intent(out) :: c
         real(wp)              :: blkmod1, blkmod2
         integer               :: q
 
         if (chemistry) then  ! Reacting mixture sound speed
-            if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
-                c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H))
+            if (avg_state == avg_state_roe .and. abs(s%c_c) > verysmall) then
+                c = sqrt(s%c_c - (s%gamma - 1.0_wp)*(s%vel_sum - s%H))
             else
-                c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
+                c = sqrt((1.0_wp + 1.0_wp/s%gamma)*s%pres/s%rho)
             end if
         else if (relativity) then  ! Relativistic sound speed
-            c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
+            c = sqrt((1._wp + 1._wp/s%gamma)*s%pres/s%rho/s%H)
         else
             if (alt_soundspeed) then  ! Wood's mixture sound speed via bulk moduli
-                blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1)
-                blkmod2 = ((gammas(2) + 1._wp)*pres + pi_infs(2))/gammas(2)
-                c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
+                blkmod1 = ((gammas(1) + 1._wp)*s%pres + pi_infs(1))/gammas(1)
+                blkmod2 = ((gammas(2) + 1._wp)*s%pres + pi_infs(2))/gammas(2)
+                c = (1._wp/(s%rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
             else if (model_eqns == model_eqns_6eq) then  ! Six-equation model sound speed
                 c = 0._wp
                 $:GPU_LOOP(parallelism='[seq]')
                 do q = 1, num_fluids
-                    c = c + adv(q)*gs_min(q)*(pres + pi_infs(q)/(gammas(q) + 1._wp))
+                    c = c + adv(q)*gs_min(q)*(s%pres + pi_infs(q)/(gammas(q) + 1._wp))
                 end do
-                c = c/rho
+                c = c/s%rho
             else if (model_eqns == model_eqns_5eq .and. bubbles_euler) then
                 ! Sound speed for bubble mixture to order O(\alpha)
 
                 if (mpp_lim .and. (num_fluids > 1)) then
-                    c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/rho
+                    c = (1._wp/s%gamma + 1._wp)*(s%pres + s%pi_inf/(s%gamma + 1._wp))/s%rho
                 else
-                    c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids)))
+                    c = (1._wp/s%gamma + 1._wp)*(s%pres + s%pi_inf/(s%gamma + 1._wp))/(s%rho*(1._wp - adv(num_fluids)))
                 end if
             else
-                c = (H - 5.e-1*vel_sum - qv/rho)/gamma
+                c = (s%H - 5.e-1*s%vel_sum - s%qv/s%rho)/s%gamma
             end if
 
             if (mixture_err .and. c < 0._wp) then

src/post_process/m_data_output.fpp

@@ -1243,6 +1243,7 @@ contains
         real(wp), dimension(num_vels) :: vel
         real(wp), dimension(num_fluids) :: adv
         integer :: i, j, k, l, s  !< looping indices
+        type(eos_state) :: eos_s
 
         Egk = 0._wp
         Elp = 0._wp
@@ -1288,7 +1289,8 @@ contains
 
                     H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho
 
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, 0._wp, 0._wp, c, qv)
+                    call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                    call s_compute_speed_of_sound(eos_s, adv, c)
 
                     Ma = maxvel/c
                     if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then

src/post_process/m_start_up.fpp

@@ -187,10 +187,11 @@ contains
              & -offset_z%beg:p + offset_z%end) :: liutex_mag
         real(wp), dimension(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end, &
              & 3) :: liutex_axis
-        integer       :: i, j, k, l, kx, ky, kz, kf, j_glb, k_glb, l_glb
-        character(50) :: filename
-        logical       :: file_exists
-        integer       :: x_beg, x_end, y_beg, y_end, z_beg, z_end
+        integer         :: i, j, k, l, kx, ky, kz, kf, j_glb, k_glb, l_glb
+        character(50)   :: filename
+        logical         :: file_exists
+        integer         :: x_beg, x_end, y_beg, y_end, z_beg, z_end
+        type(eos_state) :: eos_s
 
         if (output_partial_domain) then
             call s_define_output_region
@@ -533,8 +534,9 @@ contains
 
                         H = ((gamma_sf(i, j, k) + 1._wp)*pres + pi_inf_sf(i, j, k) + qv_sf(i, j, k))/rho_sf(i, j, k)
 
-                        call s_compute_speed_of_sound(pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), H, adv, &
-                                                      & 0._wp, 0._wp, c, qv_sf(i, j, k))
+                        call s_eos_state_roe(eos_s, pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), qv_sf(i, j, k), &
+                                             & 0._wp, H)
+                        call s_compute_speed_of_sound(eos_s, adv, c)
 
                         out%q_sf(i, j, k) = c
                     end do

src/simulation/m_cbc.fpp

@@ -511,6 +511,7 @@ contains
         real(wp)               :: Cv, Cp, e_mix, Mw, R_gas
         real(wp)               :: vel_K_sum, vel_dv_dt_sum
         integer                :: i, j, k, r  !< Generic loop iterators
+        type(eos_state)        :: eos_s
         ! Reshaping of inputted data and association of the FD and PI coefficients, or CBC coefficients, respectively, hinging on
         ! selected CBC coordinate direction
 
@@ -597,7 +598,7 @@ contains
                                     & dalpha_rho_ds, dpres_ds, dvel_dt, dadv_dt, dalpha_rho_dt, L, lambda, Ys, dYs_dt, dYs_ds, &
                                     & h_k, Cp_i, Gamma_i, Xs, drho_dt, dpres_dt, dpi_inf_dt, dqv_dt, dgamma_dt, rho, pres, E, H, &
                                     & gamma, pi_inf, qv, c, Ma, T, sum_Enthalpies, Cv, Cp, e_mix, Mw, R_gas, vel_K_sum, &
-                                    & vel_dv_dt_sum, i, j]', copyin='[dir_idx]')
+                                    & vel_dv_dt_sum, i, j, eos_s]', copyin='[dir_idx]')
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         ! Transferring the Primitive Variables
@@ -661,7 +662,8 @@ contains
                         H = (E + pres)/rho
 
                         ! Compute mixture sound speed
-                        call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv)
+                        call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_K_sum, H)
+                        call s_compute_speed_of_sound(eos_s, adv_local, c)
 
                         ! First-Order Spatial Derivatives of Primitive Variables
 

src/simulation/m_data_output.fpp

@@ -193,21 +193,23 @@ contains
         real(wp)               :: Rc_min_loc, Rc_min_glb      !< Rc stability extrema on local and global grids
         real(wp)               :: icfl, vcfl, ccfl, Rc
         integer                :: fl                          !< Fluid loop iterator
+        type(eos_state)        :: eos_s
 
         icfl_max_loc = 0._wp
         vcfl_max_loc = 0._wp
         ccfl_max_loc = 0._wp
         Rc_min_loc = huge(1.0_wp)
         ! Computing Stability Criteria at Current Time-step
         $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, icfl, &
-                            & vcfl, Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', &
-                            & reductionOp='[max, min]')
+                            & vcfl, Rc, ccfl, fl, eos_s]', reduction='[[icfl_max_loc, vcfl_max_loc, &
+                            & ccfl_max_loc], [Rc_min_loc]]', reductionOp='[max, min]')
         do l = 0, p
             do k = 0, n
                 do j = 0, m
                     call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
 
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)
+                    call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_sum, H)
+                    call s_compute_speed_of_sound(eos_s, alpha, c)
 
                     if (any_non_newtonian) then
                         Re(1) = 0._wp
@@ -1162,6 +1164,7 @@ contains
         real(wp)                        :: rad, thickness    !< For integral quantities
         logical                         :: trigger           !< For integral quantities
         real(wp)                        :: rhoYks(1:num_species)
+        type(eos_state)                 :: eos_s
 
         T = dflt_T_guess
 
@@ -1299,8 +1302,8 @@ contains
                     end if
 
                     ! Compute mixture sound Speed
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, &
-                                                  & 0._wp, c, qv)
+                    call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                    call s_compute_speed_of_sound(eos_s, alpha, c)
 
                     accel = accel_mag(j - 2, k, l)
                 end if
@@ -1382,8 +1385,8 @@ contains
                             Rdot(:) = nRdot(:)/nbub
                         end if
                         ! Compute mixture sound speed
-                        call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
-                                                      & 0._wp, 0._wp, c, qv)
+                        call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                        call s_compute_speed_of_sound(eos_s, alpha, c)
                     end if
                 end if
             else
@@ -1445,8 +1448,8 @@ contains
                             end if
 
                             ! Compute mixture sound speed
-                            call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
-                                                          & 0._wp, 0._wp, c, qv)
+                            call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                            call s_compute_speed_of_sound(eos_s, alpha, c)
 
                             accel = accel_mag(j - 2, k - 2, l - 2)
                         end if

src/simulation/m_riemann_solver_hll.fpp

@@ -97,6 +97,7 @@ contains
         type(riemann_states_vec3) :: cm        !< Conservative momentum variables
         integer :: i, j, k, l                  !< Generic loop iterators
         integer :: Re_size_loc1, Re_size_loc2  !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU
+        type(eos_state) :: eos_s_L, eos_s_R, eos_s_avg
         ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions
 
         call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
@@ -120,8 +121,8 @@ contains
                                     & Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, &
                                     & gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, damage_R, &
                                     & rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, &
-                                    & Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R]', &
-                                    & copyin='[norm_dir]', firstprivate='[Re_size_loc1, Re_size_loc2]')
+                                    & Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R, eos_s_L, &
+                                    & eos_s_R, eos_s_avg]', copyin='[norm_dir]', 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
@@ -320,17 +321,17 @@ contains
 
                             @:compute_average_state()
 
-                            call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, c_L, &
-                                                          & qv_L)
+                            eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp)
+                            call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                            call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, c_R, &
-                                                          & qv_R)
+                            eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp)
+                            call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                             !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                             ! variables are placeholders to call the subroutine.
 
-                            call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, &
-                                                          & c_sum_Yi_Phi, c_avg, qv_avg)
+                            eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, c_sum_Yi_Phi)
+                            call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                             if (mhd) then
                                 call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L)