Skip to content

Commit a85d354

Browse files
sbryngelsonsbryngelson
committed
Evaluate Mie-Gruneisen coefficients per cell along the 5-equation path
The mixture loop calls s_eos_coefficients per phase when any fluid's EOS is state dependent; the stiffened-gas arithmetic is untouched (27/27 bit-identical on the gate). The sound speed takes the partial densities as an optional argument and mixes frozen per-phase moduli, with the derivative term applied at that one site. Callers that cannot supply them are refused by the validator: 5-equation, HLL/HLLC/LF only, no bubbles, hypoelasticity, IBM, IGR, relativity, MHD, chemistry, acoustic source, probes, post-process c, characteristic BCs or Wood's law. Validation: symmetric-impact shock speed and plateau density match the Hugoniot to 1.4e-3, and a small pulse travels at the analytic sound speed to 1.9e-4 (the missing derivative term would give 24%).
1 parent ebdc35e commit a85d354

18 files changed

Lines changed: 599 additions & 71 deletions

‎examples/1D_mg_acoustic/case.py‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
"""
2+
Right-moving acoustic pulse in a single Mie-Gruneisen fluid at its reference state.
3+
The pulse is a simple wave (drho, dp = c^2 drho, du = c drho/rho0) so only the right-going
4+
characteristic carries it; the harness measures its speed against the general analytic c.
5+
"""
6+
7+
import argparse
8+
import json
9+
import math
10+
11+
parser = argparse.ArgumentParser(description="1D Mie-Gruneisen acoustic pulse")
12+
parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
13+
parser.add_argument("-N", type=int, default=200)
14+
parser.add_argument("--cfl", type=float, default=0.4)
15+
args = parser.parse_args()
16+
17+
rho0, p0, c0, s, gruneisen = 1.0, 1.0, 1.0, 1.5, 0.4
18+
c = math.sqrt(c0**2 + (1.0 + gruneisen) * p0 / rho0) # the frozen speed at the reference state
19+
amp, x0, width = 1.0e-4, 0.3, 0.05
20+
N, L, T_end = args.N, 1.0, 0.4
21+
dt = args.cfl * (L / N) / c
22+
Nt = math.ceil(T_end / dt)
23+
dt = T_end / Nt
24+
pulse = f"{amp}*exp(-((x - {x0})/{width})**2)"
25+
26+
print(
27+
json.dumps(
28+
{
29+
"run_time_info": "F",
30+
"x_domain%beg": 0.0,
31+
"x_domain%end": L,
32+
"m": N - 1,
33+
"n": 0,
34+
"p": 0,
35+
"dt": dt,
36+
"t_step_start": 0,
37+
"t_step_stop": Nt,
38+
"t_step_save": Nt,
39+
"num_patches": 1,
40+
"model_eqns": 2,
41+
"num_fluids": 1,
42+
"time_stepper": 3,
43+
"recon_type": "weno",
44+
"weno_order": 5,
45+
"weno_eps": 1.0e-16,
46+
"mapped_weno": "T",
47+
"riemann_solver": 2,
48+
"wave_speeds": 1,
49+
"avg_state": 2,
50+
"bc_x%beg": -3,
51+
"bc_x%end": -3,
52+
"format": 1,
53+
"precision": 2,
54+
"prim_vars_wrt": "T",
55+
"parallel_io": "F",
56+
"patch_icpp(1)%geometry": 1,
57+
"patch_icpp(1)%x_centroid": 0.5,
58+
"patch_icpp(1)%length_x": L,
59+
"patch_icpp(1)%alpha_rho(1)": f"{rho0} + {pulse}",
60+
"patch_icpp(1)%alpha(1)": 1.0,
61+
"patch_icpp(1)%vel(1)": f"{c}/{rho0}*{pulse}",
62+
"patch_icpp(1)%pres": f"{p0} + {c}**2*{pulse}",
63+
"fluid_pp(1)%eos": "mie_gruneisen",
64+
"fluid_pp(1)%mg_rho0": rho0,
65+
"fluid_pp(1)%mg_c0": c0,
66+
"fluid_pp(1)%mg_s": s,
67+
"fluid_pp(1)%mg_gruneisen": gruneisen,
68+
}
69+
)
70+
)

‎examples/1D_mg_impact/case.py‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
"""
2+
Symmetric impact of two Mie-Gruneisen slabs approaching at relative speed U.
3+
Each slab is brought to rest by a shock with particle-velocity jump U/2, so the shock state
4+
lies exactly on the Hugoniot u_s = c0 + s u_p; the harness checks the shock speed and the
5+
plateau density against that relation.
6+
"""
7+
8+
import argparse
9+
import json
10+
import math
11+
12+
parser = argparse.ArgumentParser(description="1D Mie-Gruneisen symmetric impact")
13+
parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
14+
parser.add_argument("-N", type=int, default=800)
15+
parser.add_argument("--U", type=float, default=1.0, help="closing speed of the two slabs")
16+
parser.add_argument("--cfl", type=float, default=0.4)
17+
args = parser.parse_args()
18+
19+
rho0, p0, c0, s, gruneisen = 1.0, 1.0e-3, 1.0, 1.5, 0.4
20+
N, L, T_end = args.N, 1.0, 0.2
21+
dt = args.cfl * (L / N) / (c0 + (s + 1.0) * args.U)
22+
Nt = math.ceil(T_end / dt)
23+
dt = T_end / Nt
24+
25+
case = {
26+
"run_time_info": "F",
27+
"x_domain%beg": 0.0,
28+
"x_domain%end": L,
29+
"m": N - 1,
30+
"n": 0,
31+
"p": 0,
32+
"dt": dt,
33+
"t_step_start": 0,
34+
"t_step_stop": Nt,
35+
"t_step_save": Nt,
36+
"num_patches": 2,
37+
"model_eqns": 2,
38+
"num_fluids": 1,
39+
"time_stepper": 3,
40+
"recon_type": "weno",
41+
"weno_order": 5,
42+
"weno_eps": 1.0e-16,
43+
"mapped_weno": "T",
44+
"riemann_solver": 2,
45+
"wave_speeds": 1,
46+
"avg_state": 2,
47+
"bc_x%beg": -3,
48+
"bc_x%end": -3,
49+
"format": 1,
50+
"precision": 2,
51+
"prim_vars_wrt": "T",
52+
"parallel_io": "F",
53+
"fluid_pp(1)%eos": "mie_gruneisen",
54+
"fluid_pp(1)%mg_rho0": rho0,
55+
"fluid_pp(1)%mg_c0": c0,
56+
"fluid_pp(1)%mg_s": s,
57+
"fluid_pp(1)%mg_gruneisen": gruneisen,
58+
}
59+
for pid, (x_c, vel) in enumerate([(0.25, 0.5 * args.U), (0.75, -0.5 * args.U)], start=1):
60+
case.update(
61+
{
62+
f"patch_icpp({pid})%geometry": 1,
63+
f"patch_icpp({pid})%x_centroid": x_c,
64+
f"patch_icpp({pid})%length_x": 0.5 * L,
65+
f"patch_icpp({pid})%alpha_rho(1)": rho0,
66+
f"patch_icpp({pid})%alpha(1)": 1.0,
67+
f"patch_icpp({pid})%vel(1)": vel,
68+
f"patch_icpp({pid})%pres": p0,
69+
}
70+
)
71+
print(json.dumps(case))

‎src/common/m_variables_conversion.fpp‎

Lines changed: 38 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1206,6 +1206,7 @@ contains
12061206
real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K
12071207
#:endif
12081208
real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K
1209+
real(wp) :: gamma_i, pi_inf_i, dpi_i
12091210
integer :: i !< Loop iterator over fluids
12101211
12111212
! The bubbly closure is written for one carrier liquid, which keeps its own coefficients
@@ -1227,8 +1228,14 @@ contains
12271228
$:GPU_LOOP(parallelism='[seq]')
12281229
do i = 1, num_fluids
12291230
rho_K = rho_K + alpha_rho_K(i)
1230-
gamma_K = gamma_K + alpha_K(i)*gammas(i)
1231-
pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i)
1231+
if (any_state_dependent_eos) then
1232+
call s_eos_coefficients(alpha_rho_K(i)/max(alpha_K(i), sgm_eps), i, gamma_i, pi_inf_i, dpi_i)
1233+
else
1234+
gamma_i = gammas(i)
1235+
pi_inf_i = pi_infs(i)
1236+
end if
1237+
gamma_K = gamma_K + alpha_K(i)*gamma_i
1238+
pi_inf_K = pi_inf_K + alpha_K(i)*pi_inf_i
12321239
qv_K = qv_K + alpha_rho_K(i)*qvs(i)
12331240
end do
12341241
end if
@@ -1473,7 +1480,7 @@ contains
14731480
14741481
!> Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H, |u|^2 and qv all cancel out of c^2
14751482
!! = ((Gamma + 1)p + Pi)/(Gamma rho). Averaged states, whose enthalpy is a free input, use the _avg variant.
1476-
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
1483+
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
14771484
14781485
$:GPU_ROUTINE(parallelism='[seq]')
14791486
@@ -1484,8 +1491,14 @@ contains
14841491
real(wp), dimension(num_fluids), intent(in) :: adv
14851492
#:endif
14861493
real(wp), intent(out) :: c
1487-
real(wp) :: alf !< Subgrid void fraction; dilute by construction
1488-
integer :: q
1494+
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
1495+
real(wp), dimension(3), intent(in), optional :: alpha_rho
1496+
#:else
1497+
real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
1498+
#:endif
1499+
real(wp) :: alf !< Subgrid void fraction; dilute by construction
1500+
real(wp) :: rho_q, gamma_q, pi_inf_q, dpi_q
1501+
integer :: q
14891502
14901503
if (chemistry) then ! Reacting mixture sound speed
14911504
c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
@@ -1494,7 +1507,16 @@ contains
14941507
else
14951508
! Every case below is a bulk modulus over a density. The equation of state enters
14961509
! only through f_bulk_modulus; the cases differ in how the phases are mixed.
1497-
if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean
1510+
if (any_state_dependent_eos .and. present(alpha_rho)) then ! frozen mixing: each phase's modulus at its own density
1511+
c = 0._wp
1512+
$:GPU_LOOP(parallelism='[seq]')
1513+
do q = 1, num_fluids
1514+
rho_q = alpha_rho(q)/max(adv(q), sgm_eps)
1515+
call s_eos_coefficients(rho_q, q, gamma_q, pi_inf_q, dpi_q)
1516+
c = c + adv(q)*(f_bulk_modulus(pres, gamma_q, pi_inf_q) - rho_q*dpi_q/gamma_q)
1517+
end do
1518+
c = c/rho
1519+
else if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean
14981520
c = 1._wp/(rho*(adv(1)/f_bulk_modulus(pres, gammas(1), pi_infs(1)) + adv(2)/f_bulk_modulus(pres, gammas(2), &
14991521
& pi_infs(2))))
15001522
else if (model_eqns == model_eqns_6eq) then ! volume-weighted arithmetic mean
@@ -1528,7 +1550,7 @@ contains
15281550
!> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its
15291551
!! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from
15301552
!! s_compute_speed_of_sound; keep the condition below in step with the branch list there.
1531-
subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c)
1553+
subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho)
15321554
15331555
$:GPU_ROUTINE(parallelism='[seq]')
15341556
@@ -1539,17 +1561,23 @@ contains
15391561
real(wp), dimension(num_fluids), intent(in) :: adv
15401562
#:endif
15411563
real(wp), intent(out) :: c
1564+
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
1565+
real(wp), dimension(3), intent(in), optional :: alpha_rho
1566+
#:else
1567+
real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
1568+
#:endif
15421569
15431570
if (chemistry) then ! Reacting mixture sound speed
15441571
if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
15451572
c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H))
15461573
else
1547-
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
1574+
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
15481575
end if
15491576
else if (relativity) then ! Relativistic sound speed
15501577
c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
1551-
else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler)) then
1552-
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
1578+
else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) &
1579+
& .or. any_state_dependent_eos) then
1580+
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
15531581
else ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives
15541582
c = (H - 5.e-1*vel_sum - qv/rho)/gamma
15551583

‎src/post_process/m_data_output.fpp‎

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

12821282
call s_compute_mixture_coefficients(alpha_rho, adv, rho, gamma, pi_inf, qv)
12831283

1284-
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
1284+
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
12851285

12861286
Ma = maxvel/c
12871287
if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then

‎src/simulation/m_data_output.fpp‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -163,11 +163,11 @@ contains
163163
real(wp) :: rho !< Cell-avg. density
164164

165165
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
166-
real(wp), dimension(3) :: alpha !< Cell-avg. volume fraction
167-
real(wp), dimension(3) :: vel !< Cell-avg. velocity
166+
real(wp), dimension(3) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density
167+
real(wp), dimension(3) :: vel !< Cell-avg. velocity
168168
#:else
169-
real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction
170-
real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity
169+
real(wp), dimension(num_fluids) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density
170+
real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity
171171
#:endif
172172
real(wp) :: vel_sum !< Cell-avg. velocity sum
173173
real(wp) :: pres !< Cell-avg. pressure
@@ -189,15 +189,15 @@ contains
189189
ccfl_max_loc = 0._wp
190190
Rc_min_loc = huge(1.0_wp)
191191
! Computing Stability Criteria at Current Time-step
192-
$:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, icfl, vcfl, &
193-
& Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', &
192+
$:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, alpha_rho, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, &
193+
& icfl, vcfl, Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', &
194194
& reductionOp='[max, min]')
195195
do l = 0, p
196196
do k = 0, n
197197
do j = 0, m
198-
call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, vel, vel_sum, qv, j, k, l)
198+
call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, k, l)
199199

200-
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c)
200+
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho)
201201

202202
if (any_non_newtonian) then
203203
Re(1) = 0._wp

‎src/simulation/m_riemann_solver_hll.fpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -333,13 +333,13 @@ contains
333333
end if
334334
end if
335335

336-
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L)
336+
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L)
337337

338-
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R)
338+
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R)
339339

340340
if (wave_speeds == wave_speeds_pressure) then
341341
call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, &
342-
& H_avg, c_sum_Yi_Phi, alpha_R, c_avg)
342+
& H_avg, c_sum_Yi_Phi, alpha_R, c_avg, alpha_rho_R)
343343
end if
344344

345345
if (mhd) then

‎src/simulation/m_riemann_solver_hllc.fpp‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -276,15 +276,15 @@ contains
276276
& qv_R, rho_avg, vel_avg_rms, H_avg, gamma_avg, qv_avg)
277277
end if
278278

279-
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L)
279+
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L)
280280

281-
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R)
281+
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R)
282282

283283
! Only the pressure-based wave-speed estimate reads the averaged state, and building it
284284
! costs eight square roots per face under the Roe average.
285285
if (wave_speeds == wave_speeds_pressure) then
286286
call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, &
287-
& H_avg, c_sum_Yi_Phi, alpha_R, c_avg)
287+
& H_avg, c_sum_Yi_Phi, alpha_R, c_avg, alpha_rho_R)
288288
end if
289289

290290
if (viscous) then
@@ -630,17 +630,17 @@ contains
630630
end do
631631
end if
632632

633-
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L)
633+
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L)
634634

635-
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R)
635+
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R)
636636

637637
! Only the pressure-based wave-speed estimate reads the averaged state, and building it
638638
! costs eight square roots per face under the Roe average.
639639
if (wave_speeds == wave_speeds_pressure) then
640640
! Zero, not c_sum_Yi_Phi: this loop never forms the chemistry average, and
641641
! chemistry with bubbles_euler/qbmm is prohibited, so the branch is unreachable.
642642
call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, &
643-
& H_avg, 0._wp, alpha_R, c_avg)
643+
& H_avg, 0._wp, alpha_R, c_avg, alpha_rho_R)
644644
end if
645645

646646
if (viscous) then
@@ -1058,15 +1058,16 @@ contains
10581058
end if
10591059
end if
10601060
1061-
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L)
1061+
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L)
10621062
1063-
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R)
1063+
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R)
10641064
10651065
! Only the pressure-based wave-speed estimate reads the averaged state, and building it
10661066
! costs eight square roots per face under the Roe average.
10671067
if (wave_speeds == wave_speeds_pressure) then
10681068
call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, &
1069-
& vel_avg_rms, H_avg, c_sum_Yi_Phi, alpha_R, c_avg)
1069+
& vel_avg_rms, H_avg, c_sum_Yi_Phi, alpha_R, c_avg, &
1070+
& alpha_rho_R)
10701071
end if
10711072
10721073
if (viscous) then

‎src/simulation/m_riemann_solver_lf.fpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -211,9 +211,9 @@ contains
211211
call s_compute_energy(pres_R, alpha_rho_R, alpha_R, vel_R_rms, E_R)
212212
end if
213213

214-
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L)
214+
call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L)
215215

216-
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R)
216+
call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R)
217217

218218
s_L = 0._wp; s_R = 0._wp
219219

0 commit comments

Comments
 (0)