Skip to content

Commit 7bd0efa

Browse files
JRChreimJRChreim
committed
Complete form of the SG EoS
This PR is to expand the SG EoS to its general form. Such a change is needed for the PC model to work on both 5 and 6 equation models
1 parent 22af239 commit 7bd0efa

File tree

8 files changed

+47
-46
lines changed

8 files changed

+47
-46
lines changed

‎src/common/m_variables_conversion.fpp‎

Lines changed: 4 additions & 10 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
@@ -1634,13 +1634,7 @@ contains
16341634
pi_infs(2))/gammas(2)
16351635
c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
16361636
elseif (model_eqns == 3) then
1637-
c = 0._wp
1638-
$:GPU_LOOP(parallelism='[seq]')
1639-
do q = 1, num_fluids
1640-
c = c + adv(q)*(1._wp/gammas(q) + 1._wp)* &
1641-
(pres + pi_infs(q)/(gammas(q) + 1._wp))
1642-
end do
1643-
c = c/rho
1637+
c = sum( adv * gs_min * ( pres + ps_inf ) ) / rho
16441638
elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles_euler))) then
16451639
! Sound speed for bubble mmixture to order O(\alpha)
16461640

@@ -1654,7 +1648,7 @@ contains
16541648
(rho*(1._wp - adv(num_fluids)))
16551649
end if
16561650
else
1657-
c = ((H - 5.e-1*vel_sum)/gamma)
1651+
c = (H - 5.e-1*vel_sum - qv/rho)/gamma
16581652
end if
16591653

16601654
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
@@ -1601,6 +1601,7 @@ contains
16011601
pres_av = 0._wp
16021602
pres = 0._wp
16031603
c = 0._wp
1604+
qv = 0._wp
16041605
16051606
do k = 0, p
16061607
do j = 0, n
@@ -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/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

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 21 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -317,6 +317,7 @@ contains
317317

318318
real(wp) :: rho_avg
319319
real(wp) :: H_avg
320+
real(wp) :: qv_avg
320321
real(wp) :: gamma_avg
321322
real(wp) :: c_avg
322323

@@ -629,16 +630,16 @@ contains
629630
@:compute_average_state()
630631

631632
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, &
632-
vel_L_rms, 0._wp, c_L)
633+
vel_L_rms, 0._wp, c_L, qv_L)
633634

634635
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, &
635-
vel_R_rms, 0._wp, c_R)
636+
vel_R_rms, 0._wp, c_R, qv_R)
636637

637638
!> The computation of c_avg does not require all the variables, and therefore the non '_avg'
638639
! variables are placeholders to call the subroutine.
639640

640641
call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, &
641-
vel_avg_rms, c_sum_Yi_Phi, c_avg)
642+
vel_avg_rms, c_sum_Yi_Phi, c_avg, qv_avg)
642643

643644
if (mhd) then
644645
call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L)
@@ -1356,10 +1357,10 @@ contains
13561357
end if
13571358

13581359
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, &
1359-
vel_L_rms, 0._wp, c_L)
1360+
vel_L_rms, 0._wp, c_L, qv_L)
13601361

13611362
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, &
1362-
vel_R_rms, 0._wp, c_R)
1363+
vel_R_rms, 0._wp, c_R, qv_R)
13631364

13641365
if (mhd) then
13651366
call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L)
@@ -1937,6 +1938,7 @@ contains
19371938
real(wp) :: rho_avg
19381939
real(wp) :: H_avg
19391940
real(wp) :: gamma_avg
1941+
real(wp) :: qv_avg
19401942
real(wp) :: c_avg
19411943

19421944
real(wp) :: s_L, s_R, s_M, s_P, s_S
@@ -2156,15 +2158,15 @@ contains
21562158
@:compute_average_state()
21572159

21582160
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, &
2159-
vel_L_rms, 0._wp, c_L)
2161+
vel_L_rms, 0._wp, c_L, qv_L)
21602162

21612163
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, &
2162-
vel_R_rms, 0._wp, c_R)
2164+
vel_R_rms, 0._wp, c_R, qv_R)
21632165

21642166
!> The computation of c_avg does not require all the variables, and therefore the non '_avg'
21652167
! variables are placeholders to call the subroutine.
21662168
call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, &
2167-
vel_avg_rms, 0._wp, c_avg)
2169+
vel_avg_rms, 0._wp, c_avg, qv_avg)
21682170

21692171
if (viscous) then
21702172
$:GPU_LOOP(parallelism='[seq]')
@@ -2472,16 +2474,16 @@ contains
24722474
@:compute_average_state()
24732475

24742476
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, &
2475-
vel_L_rms, 0._wp, c_L)
2477+
vel_L_rms, 0._wp, c_L, qv_L)
24762478

24772479
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, &
2478-
vel_R_rms, 0._wp, c_R)
2480+
vel_R_rms, 0._wp, c_R, qv_R)
24792481

24802482
!> The computation of c_avg does not require all the variables, and therefore the non '_avg'
24812483
! variables are placeholders to call the subroutine.
24822484

24832485
call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, &
2484-
vel_avg_rms, 0._wp, c_avg)
2486+
vel_avg_rms, 0._wp, c_avg, qv_avg)
24852487

24862488
if (wave_speeds == 1) then
24872489
s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R)
@@ -2858,15 +2860,15 @@ contains
28582860
end if
28592861

28602862
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, &
2861-
vel_L_rms, 0._wp, c_L)
2863+
vel_L_rms, 0._wp, c_L, qv_L)
28622864

28632865
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, &
2864-
vel_R_rms, 0._wp, c_R)
2866+
vel_R_rms, 0._wp, c_R, qv_R)
28652867

28662868
!> The computation of c_avg does not require all the variables, and therefore the non '_avg'
28672869
! variables are placeholders to call the subroutine.
28682870
call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, &
2869-
vel_avg_rms, 0._wp, c_avg)
2871+
vel_avg_rms, 0._wp, c_avg, qv_avg)
28702872

28712873
if (viscous) then
28722874
$:GPU_LOOP(parallelism='[seq]')
@@ -3295,15 +3297,15 @@ contains
32953297
@:compute_average_state()
32963298

32973299
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, &
3298-
vel_L_rms, 0._wp, c_L)
3300+
vel_L_rms, 0._wp, c_L, qv_L)
32993301

33003302
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, &
3301-
vel_R_rms, 0._wp, c_R)
3303+
vel_R_rms, 0._wp, c_R, qv_R)
33023304

33033305
!> The computation of c_avg does not require all the variables, and therefore the non '_avg'
33043306
! variables are placeholders to call the subroutine.
33053307
call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, &
3306-
vel_avg_rms, c_sum_Yi_Phi, c_avg)
3308+
vel_avg_rms, c_sum_Yi_Phi, c_avg, qv_avg)
33073309

33083310
if (viscous) then
33093311
if (chemistry) then
@@ -3750,8 +3752,8 @@ contains
37503752
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)
37513753

37523754
! (2) Compute fast wave speeds
3753-
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)
3754-
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)
3755+
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)
3756+
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)
37553757
call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L)
37563758
call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R)
37573759

‎src/simulation/m_sim_helpers.fpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -92,19 +92,20 @@ contains
9292
!! @param j x index
9393
!! @param k y index
9494
!! @param l z index
95-
subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
95+
subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
9696
$:GPU_ROUTINE(function_name='s_compute_enthalpy',parallelism='[seq]', &
9797
& cray_inline=True)
9898

9999
type(scalar_field), intent(in), dimension(sys_size) :: q_prim_vf
100100
real(wp), intent(inout), dimension(num_fluids) :: alpha
101101
real(wp), intent(inout), dimension(num_vels) :: vel
102102
real(wp), intent(inout) :: rho, gamma, pi_inf, vel_sum, H, pres
103+
real(wp), intent(out) :: qv
103104
integer, intent(in) :: j, k, l
104105
real(wp), dimension(2), intent(inout) :: Re
105106

106107
real(wp), dimension(num_fluids) :: alpha_rho, Gs
107-
real(wp) :: qv, E, G_local
108+
real(wp) :: E, G_local
108109

109110
integer :: i
110111

‎src/simulation/m_time_steppers.fpp‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -709,6 +709,7 @@ contains
709709
real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction
710710
real(wp) :: gamma !< Cell-avg. sp. heat ratio
711711
real(wp) :: pi_inf !< Cell-avg. liquid stiffness function
712+
real(wp) :: qv !< Cell-avg. fluid reference energy
712713
real(wp) :: c !< Cell-avg. sound speed
713714
real(wp) :: H !< Cell-avg. enthalpy
714715
real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers
@@ -725,18 +726,18 @@ contains
725726
idwint)
726727
end if
727728

728-
$:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H]')
729+
$:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv]')
729730
do l = 0, p
730731
do k = 0, n
731732
do j = 0, m
732733
if (igr) then
733-
call s_compute_enthalpy(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
734+
call s_compute_enthalpy(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
734735
else
735-
call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
736+
call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
736737
end if
737738

738739
! Compute mixture sound speed
739-
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c)
740+
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)
740741

741742
call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, j, k, l)
742743
end do

0 commit comments

Comments
 (0)