Skip to content

Commit 76d6a5c

Browse files
JRChreimJRChreim
authored
Merge pull request MFlowCode#1057 from JRChreim/SGEoS
Complete form of the SG EoS
2 parents 22af239 + 7e030b8 commit 76d6a5c

File tree

9 files changed

+57
-45
lines changed

9 files changed

+57
-45
lines changed

‎src/common/m_variables_conversion.fpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1601,12 +1601,12 @@ contains
16011601
end subroutine s_finalize_variables_conversion_module
16021602

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

16081608
real(wp), intent(in) :: pres
1609-
real(wp), intent(in) :: rho, gamma, pi_inf
1609+
real(wp), intent(in) :: rho, gamma, pi_inf, qv
16101610
real(wp), intent(in) :: H
16111611
real(wp), dimension(num_fluids), intent(in) :: adv
16121612
real(wp), intent(in) :: vel_sum
@@ -1642,7 +1642,7 @@ contains
16421642
end do
16431643
c = c/rho
16441644
elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles_euler))) then
1645-
! Sound speed for bubble mmixture to order O(\alpha)
1645+
! Sound speed for bubble mixture to order O(\alpha)
16461646

16471647
if (mpp_lim .and. (num_fluids > 1)) then
16481648
c = (1._wp/gamma + 1._wp)* &
@@ -1654,7 +1654,7 @@ contains
16541654
(rho*(1._wp - adv(num_fluids)))
16551655
end if
16561656
else
1657-
c = ((H - 5.e-1*vel_sum)/gamma)
1657+
c = (H - 5.e-1*vel_sum - qv/rho)/gamma
16581658
end if
16591659

16601660
if (mixture_err .and. c < 0._wp) then

‎src/post_process/m_data_output.fpp‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1582,7 +1582,7 @@ contains
15821582
impure subroutine s_write_energy_data_file(q_prim_vf, q_cons_vf)
15831583
type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf, q_cons_vf
15841584
real(wp) :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et
1585-
real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H
1585+
real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H, qv
15861586
real(wp), dimension(num_vels) :: vel
15871587
real(wp), dimension(num_fluids) :: adv
15881588
integer :: i, j, k, l, s !looping indices
@@ -1610,6 +1610,7 @@ contains
16101610
rho = 0._wp
16111611
gamma = 0._wp
16121612
pi_inf = 0._wp
1613+
qv = 0._wp
16131614
pres = q_prim_vf(E_idx)%sf(i, j, k)
16141615
Egint = Egint + q_prim_vf(E_idx + 2)%sf(i, j, k)*(fluid_pp(2)%gamma*pres)*dV
16151616
do s = 1, num_vels
@@ -1625,13 +1626,14 @@ contains
16251626
gamma = gamma + adv(l)*fluid_pp(l)%gamma
16261627
pi_inf = pi_inf + adv(l)*fluid_pp(l)%pi_inf
16271628
rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k)
1629+
qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*fluid_pp(l)%qv
16281630
end do
16291631
1630-
H = ((gamma + 1._wp)*pres + pi_inf)/rho
1632+
H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho
16311633
16321634
call s_compute_speed_of_sound(pres, rho, &
16331635
gamma, pi_inf, &
1634-
H, adv, 0._wp, 0._wp, c)
1636+
H, adv, 0._wp, 0._wp, c, qv)
16351637
16361638
Ma = maxvel/c
16371639
if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then

‎src/post_process/m_start_up.fpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -700,7 +700,7 @@ contains
700700

701701
call s_compute_speed_of_sound(pres, rho_sf(i, j, k), &
702702
gamma_sf(i, j, k), pi_inf_sf(i, j, k), &
703-
H, adv, 0._wp, 0._wp, c)
703+
H, adv, 0._wp, 0._wp, c, qv_sf(i, j, k))
704704

705705
q_sf(i, j, k) = c
706706
end do

‎src/simulation/include/inline_riemann.fpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88

99
H_avg = 5.e-1_wp*(H_L + H_R)
1010
gamma_avg = 5.e-1_wp*(gamma_L + gamma_R)
11+
qv_avg = 5.e-1_wp*(qv_L + qv_R)
1112

1213
#:enddef arithmetic_avg
1314

@@ -32,6 +33,9 @@
3233
vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2._wp/ &
3334
(sqrt(rho_L) + sqrt(rho_R))**2._wp
3435

36+
qv_avg = (sqrt(rho_L)*qv_L + sqrt(rho_R)*qv_R)/ &
37+
(sqrt(rho_L) + sqrt(rho_R))
38+
3539
if (chemistry) then
3640
eps = 0.001_wp
3741
call get_species_enthalpies_rt(T_L, h_iL)

‎src/simulation/m_cbc.fpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -854,7 +854,7 @@ contains
854854
H = (E + pres)/rho
855855
856856
! Compute mixture sound speed
857-
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c)
857+
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv)
858858
859859
! First-Order Spatial Derivatives of Primitive Variables
860860

‎src/simulation/m_data_output.fpp‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -273,19 +273,20 @@ contains
273273
real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction
274274
real(wp) :: gamma !< Cell-avg. sp. heat ratio
275275
real(wp) :: pi_inf !< Cell-avg. liquid stiffness function
276+
real(wp) :: qv !< Cell-avg. internal energy reference value
276277
real(wp) :: c !< Cell-avg. sound speed
277278
real(wp) :: H !< Cell-avg. enthalpy
278279
real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers
279280
integer :: j, k, l
280281

281282
! Computing Stability Criteria at Current Time-step
282-
$:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l,vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H]')
283+
$:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l,vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv]')
283284
do l = 0, p
284285
do k = 0, n
285286
do j = 0, m
286-
call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
287+
call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
287288

288-
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c)
289+
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)
289290

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

12931294
! Compute mixture sound Speed
12941295
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, &
1295-
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c)
1296+
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv)
12961297

12971298
accel = accel_mag(j - 2, k, l)
12981299
end if
@@ -1382,7 +1383,7 @@ contains
13821383
end if
13831384
! Compute mixture sound speed
13841385
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, &
1385-
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c)
1386+
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv)
13861387

13871388
end if
13881389
end if
@@ -1447,7 +1448,7 @@ contains
14471448

14481449
! Compute mixture sound speed
14491450
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, &
1450-
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c)
1451+
((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv)
14511452

14521453
accel = accel_mag(j - 2, k - 2, l - 2)
14531454
end if

0 commit comments

Comments
 (0)