Centralise the equation-of-state expressions in the Riemann solvers

docs/documentation/gpuParallelization.md

@@ -614,15 +614,14 @@ Does not do anything for OpenMP currently
 pure, sequential per-thread helpers is:
 
 ```fortran
-subroutine s_accumulate_mixture_properties(nf, alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K)
+subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K)
 
-    $:GPU_ROUTINE(function_name='s_accumulate_mixture_properties', parallelism='[seq]', cray_inline=True)
+    $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients', parallelism='[seq]', cray_inline=True)
 
-    integer, intent(in)                 :: nf
-    real(wp), dimension(nf), intent(in) :: alpha_rho_K, alpha_K
-    real(wp), intent(out)               :: rho_K, gamma_K, pi_inf_K, qv_K
+    real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K
+    real(wp), intent(out)                       :: rho_K, gamma_K, pi_inf_K, qv_K
     ...
-end subroutine s_accumulate_mixture_properties
+end subroutine s_compute_mixture_coefficients
 ```
 
 **When to use it.** Extract a block into a `GPU_ROUTINE` helper when:

src/common/m_variables_conversion.fpp

@@ -25,8 +25,9 @@ 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_finalize_variables_conversion_module, gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
+        & s_compute_species_fraction, s_compute_mixture_coefficients, s_compute_energy, s_compute_speed_of_sound, &
+        & s_compute_speed_of_sound_avg, s_compute_fast_magnetosonic_speed, s_finalize_variables_conversion_module, gammas, &
+        & gs_min, pi_infs, ps_inf, cvs, qvs, qvps
 
     real(wp), allocatable, dimension(:)   :: Gs_vc
     integer, allocatable, dimension(:)    :: bubrs_vc
@@ -44,7 +45,6 @@ module m_variables_conversion
     real(wp), allocatable, dimension(:,:,:), public :: rho_sf     !< Scalar density function
     real(wp), allocatable, dimension(:,:,:), public :: gamma_sf   !< Scalar sp. heat ratio function
     real(wp), allocatable, dimension(:,:,:), public :: pi_inf_sf  !< Scalar liquid stiffness function
-    real(wp), allocatable, dimension(:,:,:), public :: qv_sf      !< Scalar liquid energy reference function
 
 contains
 
@@ -155,7 +155,6 @@ contains
             rho_sf(i, j, k) = rho
             gamma_sf(i, j, k) = gamma
             pi_inf_sf(i, j, k) = pi_inf
-            qv_sf(i, j, k) = qv
         end if
 
     end subroutine s_convert_mixture_to_mixture_variables
@@ -189,7 +188,6 @@ contains
             rho_sf(k, l, r) = rho
             gamma_sf(k, l, r) = gamma
             pi_inf_sf(k, l, r) = pi_inf
-            qv_sf(k, l, r) = qv
         end if
 
     end subroutine s_convert_species_to_mixture_variables
@@ -220,29 +218,16 @@ contains
         if (present(G_K)) G_K = 0._wp
 
         ! Constrain partial densities and volume fractions within physical bounds
-        if (num_fluids == 1 .and. bubbles_euler) then
-            rho_K = alpha_rho_K(1)
-            gamma_K = gammas(1)
-            pi_inf_K = pi_infs(1)
-            qv_K = qvs(1)
-        else
-            if (mpp_lim) then
-                alpha_K_sum = 0._wp
-                do i = 1, num_fluids
-                    alpha_rho_K(i) = max(0._wp, alpha_rho_K(i))
-                    alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp)
-                    alpha_K_sum = alpha_K_sum + alpha_K(i)
-                end do
-                alpha_K = alpha_K/max(alpha_K_sum, sgm_eps)
-            end if
-            rho_K = 0._wp; gamma_K = 0._wp; pi_inf_K = 0._wp; qv_K = 0._wp
+        if (mpp_lim) then
+            alpha_K_sum = 0._wp
             do i = 1, num_fluids
-                rho_K = rho_K + alpha_rho_K(i)
-                gamma_K = gamma_K + alpha_K(i)*gammas(i)
-                pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i)
-                qv_K = qv_K + alpha_rho_K(i)*qvs(i)
+                alpha_rho_K(i) = max(0._wp, alpha_rho_K(i))
+                alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp)
+                alpha_K_sum = alpha_K_sum + alpha_K(i)
             end do
+            alpha_K = alpha_K/max(alpha_K_sum, sgm_eps)
         end if
+        call s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K)
 
         if (present(G_K)) then
             G_K = 0._wp
@@ -339,18 +324,15 @@ contains
                     allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size))
                     allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size))
                     allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size))
-                    allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size))
                 else
                     allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
                     allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
                     allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
-                    allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
                 end if
             else
                 allocate (rho_sf(-buff_size:m + buff_size,0:0,0:0))
                 allocate (gamma_sf(-buff_size:m + buff_size,0:0,0:0))
                 allocate (pi_inf_sf(-buff_size:m + buff_size,0:0,0:0))
-                allocate (qv_sf(-buff_size:m + buff_size,0:0,0:0))
             end if
         end if
 
@@ -1192,7 +1174,7 @@ contains
     !> Deallocate fluid property arrays and post-processing fields allocated during module initialization.
     impure subroutine s_finalize_variables_conversion_module()
 
-        if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf, qv_sf)
+        if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf)
 
         @:DEALLOCATE(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs_vc)
         if (allocated(bubrs_vc)) then
@@ -1204,33 +1186,88 @@ contains
 
     end subroutine s_finalize_variables_conversion_module
 
-    !> 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)
+    !> Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot aliases the void fraction
+    !! (eqn_idx%alf == eqn_idx%adv%end), so alpha is not a composition there and the coefficients are the liquid's. Clipping stays
+    !! with callers; it differs between solvers and cannot coincide with that case, as mpp_lim requires num_fluids > 1.
+    subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K)
+
+        $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients', parallelism='[seq]', cray_inline=True)
+
+        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+            real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K
+        #:else
+            real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K
+        #:endif
+        real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K
+        integer               :: i  !< Loop iterator over fluids
+
+        if (num_fluids == 1 .and. bubbles_euler) then
+            rho_K = alpha_rho_K(1)
+            gamma_K = gammas(1)
+            pi_inf_K = pi_infs(1)
+            qv_K = qvs(1)
+        else
+            rho_K = 0._wp
+            gamma_K = 0._wp
+            pi_inf_K = 0._wp
+            qv_K = 0._wp
+
+            $:GPU_LOOP(parallelism='[seq]')
+            do i = 1, num_fluids
+                rho_K = rho_K + alpha_rho_K(i)
+                gamma_K = gamma_K + alpha_K(i)*gammas(i)
+                pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i)
+                qv_K = qv_K + alpha_rho_K(i)*qvs(i)
+            end do
+        end if
+
+    end subroutine s_compute_mixture_coefficients
+
+    !> Total energy per unit volume, thermodynamic terms only. Callers add magnetic and elastic energy, which are not
+    !! equation-of-state terms. The chemistry and relativistic branches use a different relation and stay open-coded.
+    subroutine s_compute_energy(pres, alpha_rho_K, alpha_K, vel_sum, E)
+
+        $:GPU_ROUTINE(function_name='s_compute_energy', parallelism='[seq]', cray_inline=True)
+
+        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+            real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K
+        #:else
+            real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K
+        #:endif
+        real(wp), intent(in)  :: pres, vel_sum
+        real(wp), intent(out) :: E
+        real(wp)              :: rho, gamma, pi_inf, qv
+
+        call s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho, gamma, pi_inf, qv)
+
+        E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_sum + qv
+
+    end subroutine s_compute_energy
+
+    !> Compute the speed of sound of a thermodynamic state.
+    !!
+    !! Enthalpy is not an argument: substituting H = ((Gamma + 1)p + Pi + qv)/rho + |u|^2/2 into
+    !! c^2 = (H - |u|^2/2 - qv/rho)/Gamma leaves c^2 = ((Gamma + 1)p + Pi)/(Gamma rho), so H, |u|^2
+    !! and qv all cancel. Averaged states, whose enthalpy is a free input, use
+    !! s_compute_speed_of_sound_avg.
+    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, 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
+        real(wp), intent(in) :: pres, rho, gamma, pi_inf
         #: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))
-            else
-                c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
-            end if
-        else if (relativity) then  ! Relativistic sound speed
-            c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
+            c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
+        else if (relativity) then  ! Relativistic sound speed, whose enthalpy is 1 + (Gamma + 1)p/rho
+            c = sqrt((1._wp + 1._wp/gamma)*pres/rho/(1._wp + (gamma + 1._wp)*pres/rho))
         else
             if (alt_soundspeed) then  ! Wood's mixture sound speed via bulk moduli
                 blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1)
@@ -1251,8 +1288,8 @@ contains
                 else
                     c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids)))
                 end if
-            else
-                c = (H - 5.e-1*vel_sum - qv/rho)/gamma
+            else  ! Stiffened-gas mixture, with H, |u|^2 and qv cancelled out
+                c = ((gamma + 1._wp)*pres + pi_inf)/(gamma*rho)
             end if
 
             if (mixture_err .and. c < 0._wp) then
@@ -1264,6 +1301,46 @@ contains
 
     end subroutine s_compute_speed_of_sound
 
+    !> Compute the speed of sound of an interface-averaged state.
+    !!
+    !! An average of two states is not a state: its enthalpy is not the one its pressure and density
+    !! imply, so the caller supplies it, along with |u|^2 and qv. Only the enthalpy-reading branches
+    !! differ; the rest defer to s_compute_speed_of_sound. Keep the condition below in step with the
+    !! branch list there.
+    subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c
+        #: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(out) :: c
+
+        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))
+            else
+                call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
+            end if
+        else if (relativity) then  ! Relativistic sound speed
+            c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
+        else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler)) then
+            call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
+        else  ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives
+            c = (H - 5.e-1*vel_sum - qv/rho)/gamma
+
+            if (mixture_err .and. c < 0._wp) then
+                c = 100._wp*sgm_eps
+            else
+                c = sqrt(c)
+            end if
+        end if
+
+    end subroutine s_compute_speed_of_sound_avg
+
     !> Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components.
     subroutine s_compute_fast_magnetosonic_speed(rho, c, B, norm, c_fast, h)
 

src/post_process/m_data_output.fpp

@@ -1238,11 +1238,11 @@ contains
     impure subroutine s_write_energy_data_file(q_prim_vf, q_cons_vf)
 
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf, q_cons_vf
-        real(wp) :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et
-        real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H, qv
-        real(wp), dimension(num_vels) :: vel
-        real(wp), dimension(num_fluids) :: adv
-        integer :: i, j, k, l, s  !< looping indices
+        real(wp)                                            :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et
+        real(wp)                                            :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma
+        real(wp), dimension(num_vels)                       :: vel
+        real(wp), dimension(num_fluids)                     :: adv
+        integer                                             :: i, j, k, l, s  !< looping indices
 
         Egk = 0._wp
         Elp = 0._wp
@@ -1267,7 +1267,6 @@ contains
                     rho = 0._wp
                     gamma = 0._wp
                     pi_inf = 0._wp
-                    qv = 0._wp
                     pres = q_prim_vf(eqn_idx%E)%sf(i, j, k)
                     Egint = Egint + q_prim_vf(eqn_idx%E + 2)%sf(i, j, k)*(gammas(2)*pres)*dV
                     do s = 1, num_vels
@@ -1283,12 +1282,9 @@ contains
                         gamma = gamma + adv(l)*gammas(l)
                         pi_inf = pi_inf + adv(l)*pi_infs(l)
                         rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k)
-                        qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*qvs(l)
                     end do
 
-                    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_compute_speed_of_sound(pres, rho, gamma, pi_inf, 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

@@ -177,11 +177,11 @@ contains
     end subroutine s_perform_time_step
 
     !> Derive requested flow quantities from primitive variables and write them to the formatted database files.
-    impure subroutine s_save_data(t_step, varname, pres, c, H)
+    impure subroutine s_save_data(t_step, varname, pres, c)
 
         integer, intent(inout)                 :: t_step
         character(LEN=name_len), intent(inout) :: varname
-        real(wp), intent(inout)                :: pres, c, H
+        real(wp), intent(inout)                :: pres, c
 
         real(wp), dimension(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end, &
              & -offset_z%beg:p + offset_z%end) :: liutex_mag
@@ -531,10 +531,7 @@ contains
 
                         pres = q_prim_vf(eqn_idx%E)%sf(i, j, k)
 
-                        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_compute_speed_of_sound(pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), adv, c)
 
                         out%q_sf(i, j, k) = c
                     end do

src/post_process/p_main.fpp

@@ -15,7 +15,6 @@ program p_main
     character(LEN=name_len) :: varname
     real(wp)                :: pres
     real(wp)                :: c
-    real(wp)                :: H
     real(wp)                :: start, finish
 
     call s_initialize_mpi_domain()
@@ -41,7 +40,7 @@ program p_main
 
         call s_perform_time_step(t_step)
 
-        call s_save_data(t_step, varname, pres, c, H)
+        call s_save_data(t_step, varname, pres, c)
 
         call cpu_time(finish)
 

src/simulation/m_cbc.fpp

@@ -501,7 +501,6 @@ contains
         real(wp)               :: rho         !< Cell averaged density
         real(wp)               :: pres        !< Cell averaged pressure
         real(wp)               :: E           !< Cell averaged energy
-        real(wp)               :: H           !< Cell averaged enthalpy
         real(wp)               :: gamma       !< Cell averaged specific heat ratio
         real(wp)               :: pi_inf      !< Cell averaged liquid stiffness
         real(wp)               :: qv          !< Cell averaged fluid reference energy
@@ -595,7 +594,7 @@ contains
                 ! FD2 or FD4 of RHS at j = 0
                 $:GPU_PARALLEL_LOOP(collapse=2, private='[r, k, alpha_rho, vel, adv_local, mf, dvel_ds, dadv_ds, Re_cbc, &
                                     & 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, &
+                                    & h_k, Cp_i, Gamma_i, Xs, drho_dt, dpres_dt, dpi_inf_dt, dqv_dt, dgamma_dt, rho, pres, E, &
                                     & 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]')
                 do r = is3%beg, is3%end
@@ -655,13 +654,10 @@ contains
                                 gamma = 1.0_wp/(Cp/Cv - 1.0_wp)
                             end if
                         else
-                            E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_K_sum
                         end if
 
-                        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_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv_local, c)
 
                         ! First-Order Spatial Derivatives of Primitive Variables