Skip to content

Complete form of the SG EoS #1057

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
JRChreim merged 9 commits into from
Nov 24, 2025
Merged

Complete form of the SG EoS #1057

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
7bd0efa
Complete form of the SG EoS
JRChreim Nov 22, 2025
62d430e
Update m_variables_conversion.fpp
JRChreim Nov 22, 2025
c00ba3d
Update m_data_output.fpp
JRChreim Nov 23, 2025
9607e33
Update m_variables_conversion.fpp
JRChreim Nov 23, 2025
50fdb51
fix debug build
wilfonba Nov 23, 2025
ec8834a
correct roe average of qv
wilfonba Nov 23, 2025
d7fd48f
fix acc build
wilfonba Nov 23, 2025
98ce554
Update m_riemann_solvers.fpp
JRChreim Nov 23, 2025
7e030b8
Update m_riemann_solvers.fpp
JRChreim Nov 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 4 additions & 10 deletions src/common/m_variables_conversion.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1601,12 +1601,12 @@ contains
end subroutine s_finalize_variables_conversion_module

#ifndef MFC_PRE_PROCESS
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c)
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv)
$:GPU_ROUTINE(function_name='s_compute_speed_of_sound', &
& parallelism='[seq]', cray_inline=True)

real(wp), intent(in) :: pres
real(wp), intent(in) :: rho, gamma, pi_inf
real(wp), intent(in) :: rho, gamma, pi_inf, qv
real(wp), intent(in) :: H
real(wp), dimension(num_fluids), intent(in) :: adv
real(wp), intent(in) :: vel_sum
Expand Down Expand Up @@ -1634,13 +1634,7 @@ contains
pi_infs(2))/gammas(2)
c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
elseif (model_eqns == 3) then
c = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids
c = c + adv(q)*(1._wp/gammas(q) + 1._wp)* &
(pres + pi_infs(q)/(gammas(q) + 1._wp))
end do
c = c/rho
c = sum(adv*gs_min*(pres + ps_inf))/rho
elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles_euler))) then
! Sound speed for bubble mmixture to order O(\alpha)

Expand All @@ -1654,7 +1648,7 @@ contains
(rho*(1._wp - adv(num_fluids)))
end if
else
c = ((H - 5.e-1*vel_sum)/gamma)
c = (H - 5.e-1*vel_sum - qv/rho)/gamma
end if

if (mixture_err .and. c < 0._wp) then
Expand Down
8 changes: 5 additions & 3 deletions src/post_process/m_data_output.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1582,7 +1582,7 @@ 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
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
Expand All @@ -1601,6 +1601,7 @@ contains
pres_av = 0._wp
pres = 0._wp
c = 0._wp
qv = 0._wp

do k = 0, p
do j = 0, n
Expand All @@ -1625,13 +1626,14 @@ contains
gamma = gamma + adv(l)*fluid_pp(l)%gamma
pi_inf = pi_inf + adv(l)*fluid_pp(l)%pi_inf
rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k)
qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*fluid_pp(l)%qv
end do

H = ((gamma + 1._wp)*pres + pi_inf)/rho
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)
H, adv, 0._wp, 0._wp, c, qv)

Ma = maxvel/c
if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then
Expand Down
2 changes: 1 addition & 1 deletion src/post_process/m_start_up.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -700,7 +700,7 @@ contains

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)
H, adv, 0._wp, 0._wp, c, qv_sf(i, j, k))

q_sf(i, j, k) = c
end do
Expand Down
2 changes: 1 addition & 1 deletion src/simulation/m_cbc.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -854,7 +854,7 @@ 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)
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv)

! First-Order Spatial Derivatives of Primitive Variables

Expand Down
13 changes: 7 additions & 6 deletions src/simulation/m_data_output.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -273,19 +273,20 @@ contains
real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction
real(wp) :: gamma !< Cell-avg. sp. heat ratio
real(wp) :: pi_inf !< Cell-avg. liquid stiffness function
real(wp) :: qv !< Cell-avg. internal energy reference value
real(wp) :: c !< Cell-avg. sound speed
real(wp) :: H !< Cell-avg. enthalpy
real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers
integer :: j, k, l

! 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]')
$:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l,vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv]')
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, j, k, l)
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)
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)

if (viscous) then
call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf)
Expand Down Expand Up @@ -1292,7 +1293,7 @@ contains

! 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)
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv)

accel = accel_mag(j - 2, k, l)
end if
Expand Down Expand Up @@ -1382,7 +1383,7 @@ 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)
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv)

end if
end if
Expand Down Expand Up @@ -1447,7 +1448,7 @@ contains

! 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)
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv)

accel = accel_mag(j - 2, k - 2, l - 2)
end if
Expand Down
40 changes: 21 additions & 19 deletions src/simulation/m_riemann_solvers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,6 +317,7 @@ contains

real(wp) :: rho_avg
real(wp) :: H_avg
real(wp) :: qv_avg
real(wp) :: gamma_avg
real(wp) :: c_avg

Expand Down Expand Up @@ -629,16 +630,16 @@ 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)
vel_L_rms, 0._wp, c_L, qv_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)
vel_R_rms, 0._wp, c_R, qv_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)
vel_avg_rms, c_sum_Yi_Phi, c_avg, qv_avg)

if (mhd) then
call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L)
Expand Down Expand Up @@ -1356,10 +1357,10 @@ contains
end if

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)
vel_L_rms, 0._wp, c_L, qv_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)
vel_R_rms, 0._wp, c_R, qv_R)

if (mhd) then
call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L)
Expand Down Expand Up @@ -1937,6 +1938,7 @@ contains
real(wp) :: rho_avg
real(wp) :: H_avg
real(wp) :: gamma_avg
real(wp) :: qv_avg
real(wp) :: c_avg

real(wp) :: s_L, s_R, s_M, s_P, s_S
Expand Down Expand Up @@ -2156,15 +2158,15 @@ 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)
vel_L_rms, 0._wp, c_L, qv_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)
vel_R_rms, 0._wp, c_R, qv_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, 0._wp, c_avg)
vel_avg_rms, 0._wp, c_avg, qv_avg)

if (viscous) then
$:GPU_LOOP(parallelism='[seq]')
Expand Down Expand Up @@ -2472,16 +2474,16 @@ 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)
vel_L_rms, 0._wp, c_L, qv_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)
vel_R_rms, 0._wp, c_R, qv_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, 0._wp, c_avg)
vel_avg_rms, 0._wp, c_avg, qv_avg)

if (wave_speeds == 1) then
s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R)
Expand Down Expand Up @@ -2858,15 +2860,15 @@ contains
end if

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)
vel_L_rms, 0._wp, c_L, qv_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)
vel_R_rms, 0._wp, c_R, qv_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, 0._wp, c_avg)
vel_avg_rms, 0._wp, c_avg, qv_avg)

if (viscous) then
$:GPU_LOOP(parallelism='[seq]')
Expand Down Expand Up @@ -3295,15 +3297,15 @@ 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)
vel_L_rms, 0._wp, c_L, qv_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)
vel_R_rms, 0._wp, c_R, qv_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)
vel_avg_rms, c_sum_Yi_Phi, c_avg, qv_avg)

if (viscous) then
if (chemistry) then
Expand Down Expand Up @@ -3750,8 +3752,8 @@ contains
H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)

! (2) Compute fast wave speeds
call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, 0._wp, c%L)
call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, 0._wp, c%R)
call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, 0._wp, c%L, qv%L)
call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, 0._wp, c%R, qv%R)
call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L)
call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R)

Expand Down
5 changes: 3 additions & 2 deletions src/simulation/m_sim_helpers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,19 +92,20 @@ contains
!! @param j x index
!! @param k y index
!! @param l z index
subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
$:GPU_ROUTINE(function_name='s_compute_enthalpy',parallelism='[seq]', &
& cray_inline=True)

type(scalar_field), intent(in), dimension(sys_size) :: q_prim_vf
real(wp), intent(inout), dimension(num_fluids) :: alpha
real(wp), intent(inout), dimension(num_vels) :: vel
real(wp), intent(inout) :: rho, gamma, pi_inf, vel_sum, H, pres
real(wp), intent(out) :: qv
integer, intent(in) :: j, k, l
real(wp), dimension(2), intent(inout) :: Re

real(wp), dimension(num_fluids) :: alpha_rho, Gs
real(wp) :: qv, E, G_local
real(wp) :: E, G_local

integer :: i

Expand Down
9 changes: 5 additions & 4 deletions src/simulation/m_time_steppers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -709,6 +709,7 @@ contains
real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction
real(wp) :: gamma !< Cell-avg. sp. heat ratio
real(wp) :: pi_inf !< Cell-avg. liquid stiffness function
real(wp) :: qv !< Cell-avg. fluid reference energy
real(wp) :: c !< Cell-avg. sound speed
real(wp) :: H !< Cell-avg. enthalpy
real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers
Expand All @@ -725,18 +726,18 @@ contains
idwint)
end if

$:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H]')
$:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv]')
do l = 0, p
do k = 0, n
do j = 0, m
if (igr) then
call s_compute_enthalpy(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
call s_compute_enthalpy(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
else
call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
end if

! Compute mixture sound speed
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c)
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)

call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, j, k, l)
end do
Expand Down
Loading