Skip to content

Commit 5e9d629

Browse files
Fahad NabidFahad Nabid
Fahad Nabid
authored and
Fahad Nabid
committed
JWL solver wiring: cons-prim conversions, sound speed, Riemann solvers (5-eq model)
Wire the JWL EOS into the solver hot path for the five-equation model (model_eqns = 2, Allaire et al. JCP 2002). All changes are guarded by jwl_idx > 0 and model_eqns == model_eqns_5eq so non-JWL cases are unaffected. m_variables_conversion: - s_convert_conservative_to_primitive_variables: after s_compute_pressure gives the stiffened-gas pressure, override it with s_jwl_mix_pressure_er using Y = alpha_rho(jwl_idx)/rho and alpha = alpha(jwl_idx). Without this, the solver uses the wrong EOS to recover pressure from conserved variables, giving thermodynamically inconsistent primitive fields. - s_convert_primitive_to_conservative_variables: add JWL branch that computes E = rho*e_mix(rho,p,Y,alpha) + 0.5*rho|u|^2 + qv via s_jwl_mix_energy_pr instead of the stiffened-gas E = gamma*p + pi_inf formula, keeping prim->cons consistent with cons->prim. - s_compute_speed_of_sound: add JWL branch using s_jwl_mixture_sound_speed_squared (frozen mass-weighted rule). Accepts an optional alpha_rho_j argument for the correct Y_j = alpha_rho_j/rho; falls back to the rho0 proxy alpha_j*rho0/rho when not supplied (e.g. avg-state calls). - Export s_jwl_mixture_sound_speed_squared via the module public list. m_riemann_solver_hll / m_riemann_solver_hllc: - Add use m_jwl to both modules. - In every E_L/E_R reconstruction block: when jwl_idx > 0 and 5-eq, compute e_L/e_R from pressure via s_jwl_mix_energy_pr and form E = rho*e + 0.5*rho|u|^2 instead of the stiffened-gas gamma*p + pi_inf formula. This keeps the Riemann wave-speed estimates thermodynamically consistent with the JWL EOS for both HLL and HLLC solvers. - In m_riemann_solver_hll: pass alpha_rho_j to both s_compute_speed_of_sound calls so the sound speed uses the actual phasic partial density for Y_j rather than the rho0 proxy.
1 parent d83068f commit 5e9d629

3 files changed

Lines changed: 130 additions & 18 deletions

File tree

‎src/common/m_variables_conversion.fpp‎

Lines changed: 49 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ module m_variables_conversion
3535
s_compute_pressure, &
3636
s_jwl_pcold, &
3737
s_jwl_sound_speed_squared, &
38+
s_jwl_mixture_sound_speed_squared, &
3839
s_jwl_energy_pr, &
3940
s_jwl_mix_pressure_er, &
4041
s_jwl_mix_energy_pr, &
@@ -693,6 +694,18 @@ contains
693694
call s_compute_pressure(qK_cons_vf(eqn_idx%E)%sf(j, k, l), qK_cons_vf(eqn_idx%alf)%sf(j, k, l), dyn_pres_K, &
694695
& pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres, T, pres_mag=pres_mag)
695696

697+
! JWL override: when a JWL fluid is present in a 5-equation model, replace the stiffened-gas
698+
! pressure with the correct JWL mixture inversion. alpha_rho_K(jwl_idx)/rho_K gives the true
699+
! products mass fraction Y; alpha_K(jwl_idx) gives the JWL volume fraction. Both are already
700+
! assembled (and clamped via mpp_lim) by s_convert_species_to_mixture_variables_acc above.
701+
#:if not chemistry
702+
if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq .and. (.not. bubbles_euler)) then
703+
call s_jwl_mix_pressure_er(rho_K, (real(qK_cons_vf(eqn_idx%E)%sf(j, k, l), &
704+
& wp) - dyn_pres_K)/max(rho_K, sgm_eps), alpha_rho_K(jwl_idx)/max(rho_K, &
705+
& sgm_eps), alpha_K(jwl_idx), jwl_idx, pres)
706+
end if
707+
#:endif
708+
696709
qK_prim_vf(eqn_idx%E)%sf(j, k, l) = pres
697710

698711
if (chemistry) then
@@ -939,6 +952,16 @@ contains
939952
! MHD energy includes magnetic pressure contribution
940953
q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, &
941954
& l) + dyn_pres + pres_mag + pi_inf + qv
955+
else if (jwl_idx > 0 .and. (model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then
956+
! JWL five-equation model: E = rho*e_mix(rho,p,Y,alpha) + 0.5*rho*|u|^2 + qv.
957+
! Y = alpha_rho_j/rho, alpha_j from the primitive adv slot.
958+
block
959+
real(wp) :: e_mix_jwl, Y_j, alpha_j
960+
alpha_j = q_prim_vf(eqn_idx%adv%beg + jwl_idx - 1)%sf(j, k, l)
961+
Y_j = q_prim_vf(jwl_idx)%sf(j, k, l)/max(rho, sgm_eps)
962+
call s_jwl_mix_energy_pr(rho, q_prim_vf(eqn_idx%E)%sf(j, k, l), Y_j, alpha_j, jwl_idx, e_mix_jwl)
963+
q_cons_vf(eqn_idx%E)%sf(j, k, l) = rho*e_mix_jwl + dyn_pres + qv
964+
end block
942965
else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then
943966
! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv
944967
q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + pi_inf + qv
@@ -1268,7 +1291,7 @@ contains
12681291

12691292
#ifndef MFC_PRE_PROCESS
12701293
!> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models.
1271-
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv)
1294+
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv, alpha_rho_j)
12721295

12731296
$:GPU_ROUTINE(parallelism='[seq]')
12741297

@@ -1280,11 +1303,12 @@ contains
12801303
#:else
12811304
real(wp), dimension(num_fluids), intent(in) :: adv
12821305
#:endif
1283-
real(wp), intent(in) :: vel_sum
1284-
real(wp), intent(in) :: c_c
1285-
real(wp), intent(out) :: c
1286-
real(wp) :: blkmod1, blkmod2
1287-
integer :: q
1306+
real(wp), intent(in) :: vel_sum
1307+
real(wp), intent(in) :: c_c
1308+
real(wp), intent(out) :: c
1309+
real(wp), intent(in), optional :: alpha_rho_j !< JWL partial density alpha_j*rho_j for correct Y_j
1310+
real(wp) :: blkmod1, blkmod2
1311+
integer :: q
12881312

12891313
if (chemistry) then ! Reacting mixture sound speed
12901314
if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
@@ -1294,6 +1318,25 @@ contains
12941318
end if
12951319
else if (relativity) then ! Relativistic sound speed
12961320
c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
1321+
else if (jwl_idx > 0) then
1322+
! JWL mixture sound speed: frozen (mass-weighted) rule.
1323+
! Y_j = alpha_rho_j / rho (correct: uses the actual phasic partial density).
1324+
! If alpha_rho_j is not passed (e.g. avg-state call), fall back to the rho0 proxy
1325+
! alpha_j*rho0/rho which is exact at initial conditions and only a wave-speed bound.
1326+
block
1327+
real(wp) :: c2, Y_j, alpha_j
1328+
alpha_j = min(max(adv(jwl_idx), 0._wp), 1._wp)
1329+
if (present(alpha_rho_j)) then
1330+
Y_j = min(max(alpha_rho_j/max(rho, sgm_eps), 0._wp), 1._wp)
1331+
else
1332+
Y_j = min(max(alpha_j*jwl_rho0s(jwl_idx)/max(rho, sgm_eps), 0._wp), 1._wp)
1333+
end if
1334+
call s_jwl_mixture_sound_speed_squared(rho, pres, Y_j, alpha_j, jwl_As(jwl_idx), jwl_Bs(jwl_idx), &
1335+
& jwl_R1s(jwl_idx), jwl_R2s(jwl_idx), jwl_omegas(jwl_idx), &
1336+
& jwl_rho0s(jwl_idx), jwl_E0s(jwl_idx), jwl_air_e0s(jwl_idx), &
1337+
& jwl_air_rho0s(jwl_idx), jwl_air_gammas(jwl_idx), c2)
1338+
c = sqrt(max(c2, sgm_eps))
1339+
end block
12971340
else
12981341
if (alt_soundspeed) then ! Wood's mixture sound speed via bulk moduli
12991342
blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1)

‎src/simulation/m_riemann_solver_hll.fpp‎

Lines changed: 20 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ module m_riemann_solver_hll
1212
use m_derived_types
1313
use m_global_parameters
1414
use m_variables_conversion
15+
use m_jwl
1516
use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs, avg_state_roe, &
1617
& avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure
1718
use m_chemistry
@@ -315,8 +316,23 @@ contains
315316
! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
316317
H_R = (E_R + pres_R - pres_mag%R)/rho_R
317318
else
318-
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L
319-
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R
319+
if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq) then
320+
! JWL: reconstruct total energy from p(rho,Y,alpha) via the inverse EOS.
321+
! alpha_rho_L/R(jwl_idx) are the JWL partial densities; alpha_L/R(jwl_idx) are
322+
! the JWL volume fractions, both assembled from primitives in the loops above.
323+
block
324+
real(wp) :: e_L, e_R, Y_jL, Y_jR
325+
Y_jL = alpha_rho_L(jwl_idx)/max(rho_L, sgm_eps)
326+
Y_jR = alpha_rho_R(jwl_idx)/max(rho_R, sgm_eps)
327+
call s_jwl_mix_energy_pr(rho_L, pres_L, Y_jL, alpha_L(jwl_idx), jwl_idx, e_L)
328+
call s_jwl_mix_energy_pr(rho_R, pres_R, Y_jR, alpha_R(jwl_idx), jwl_idx, e_R)
329+
E_L = rho_L*e_L + 5.e-1*rho_L*vel_L_rms + qv_L
330+
E_R = rho_R*e_R + 5.e-1*rho_R*vel_R_rms + qv_R
331+
end block
332+
else
333+
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L
334+
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R
335+
end if
320336
H_L = (E_L + pres_L)/rho_L
321337
H_R = (E_R + pres_R)/rho_R
322338
end if
@@ -356,10 +372,10 @@ contains
356372
@:compute_average_state()
357373

358374
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, &
359-
& qv_L)
375+
& qv_L, alpha_rho_j=alpha_rho_L(jwl_idx))
360376

361377
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, &
362-
& qv_R)
378+
& qv_R, alpha_rho_j=alpha_rho_R(jwl_idx))
363379

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

‎src/simulation/m_riemann_solver_hllc.fpp‎

Lines changed: 61 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ module m_riemann_solver_hllc
1212
use m_derived_types
1313
use m_global_parameters
1414
use m_variables_conversion
15+
use m_jwl
1516
use m_bubbles
1617
use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs, model_eqns_5eq, &
1718
& model_eqns_6eq, model_eqns_4eq, avg_state_roe, avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure
@@ -242,8 +243,22 @@ contains
242243
end do
243244
end if
244245

245-
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L
246-
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R
246+
if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq) then
247+
! JWL: reconstruct total energy from p(rho,Y,alpha) via the inverse EOS,
248+
! then add kinetic energy. alpha_rho_j = qL/R_prim_rsx_vf(*,jwl_idx).
249+
block
250+
real(wp) :: e_L, e_R, Y_jL, Y_jR
251+
Y_jL = qL_prim_rsx_vf(${SF('')}$, jwl_idx)/max(rho_L, sgm_eps)
252+
Y_jR = qR_prim_rsx_vf(${SF(' + 1')}$, jwl_idx)/max(rho_R, sgm_eps)
253+
call s_jwl_mix_energy_pr(rho_L, pres_L, Y_jL, alpha_L(jwl_idx), jwl_idx, e_L)
254+
call s_jwl_mix_energy_pr(rho_R, pres_R, Y_jR, alpha_R(jwl_idx), jwl_idx, e_R)
255+
E_L = rho_L*e_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L
256+
E_R = rho_R*e_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R
257+
end block
258+
else
259+
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L
260+
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R
261+
end if
247262

248263
! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY
249264
if (hypoelasticity) then
@@ -612,8 +627,22 @@ contains
612627
pres_L = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E)
613628
pres_R = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E)
614629

615-
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L
616-
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R
630+
if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq) then
631+
! JWL: reconstruct total energy from p(rho,Y,alpha) via the inverse EOS,
632+
! then add kinetic energy. alpha_rho_j = qL/R_prim_rsx_vf(*,jwl_idx).
633+
block
634+
real(wp) :: e_L, e_R, Y_jL, Y_jR
635+
Y_jL = qL_prim_rsx_vf(${SF('')}$, jwl_idx)/max(rho_L, sgm_eps)
636+
Y_jR = qR_prim_rsx_vf(${SF(' + 1')}$, jwl_idx)/max(rho_R, sgm_eps)
637+
call s_jwl_mix_energy_pr(rho_L, pres_L, Y_jL, alpha_L(jwl_idx), jwl_idx, e_L)
638+
call s_jwl_mix_energy_pr(rho_R, pres_R, Y_jR, alpha_R(jwl_idx), jwl_idx, e_R)
639+
E_L = rho_L*e_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L
640+
E_R = rho_R*e_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R
641+
end block
642+
else
643+
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L
644+
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R
645+
end if
617646

618647
H_L = (E_L + pres_L)/rho_L
619648
H_R = (E_R + pres_R)/rho_R
@@ -871,8 +900,20 @@ contains
871900
pres_L = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E)
872901
pres_R = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E)
873902

874-
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms
875-
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms
903+
if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq) then
904+
block
905+
real(wp) :: e_L, e_R, Y_jL, Y_jR
906+
Y_jL = qL_prim_rsx_vf(${SF('')}$, jwl_idx)/max(rho_L, sgm_eps)
907+
Y_jR = qR_prim_rsx_vf(${SF(' + 1')}$, jwl_idx)/max(rho_R, sgm_eps)
908+
call s_jwl_mix_energy_pr(rho_L, pres_L, Y_jL, alpha_L(jwl_idx), jwl_idx, e_L)
909+
call s_jwl_mix_energy_pr(rho_R, pres_R, Y_jR, alpha_R(jwl_idx), jwl_idx, e_R)
910+
E_L = rho_L*e_L + 5.e-1_wp*rho_L*vel_L_rms
911+
E_R = rho_R*e_R + 5.e-1_wp*rho_R*vel_R_rms
912+
end block
913+
else
914+
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms
915+
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms
916+
end if
876917

877918
H_L = (E_L + pres_L)/rho_L
878919
H_R = (E_R + pres_R)/rho_R
@@ -1313,8 +1354,20 @@ contains
13131354
H_L = (E_L + pres_L)/rho_L
13141355
H_R = (E_R + pres_R)/rho_R
13151356
else
1316-
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L
1317-
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R
1357+
if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq) then
1358+
block
1359+
real(wp) :: e_L, e_R, Y_jL, Y_jR
1360+
Y_jL = qL_prim_rsx_vf(${SF('')}$, jwl_idx)/max(rho_L, sgm_eps)
1361+
Y_jR = qR_prim_rsx_vf(${SF(' + 1')}$, jwl_idx)/max(rho_R, sgm_eps)
1362+
call s_jwl_mix_energy_pr(rho_L, pres_L, Y_jL, alpha_L(jwl_idx), jwl_idx, e_L)
1363+
call s_jwl_mix_energy_pr(rho_R, pres_R, Y_jR, alpha_R(jwl_idx), jwl_idx, e_R)
1364+
E_L = rho_L*e_L + 5.e-1*rho_L*vel_L_rms + qv_L
1365+
E_R = rho_R*e_R + 5.e-1*rho_R*vel_R_rms + qv_R
1366+
end block
1367+
else
1368+
E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L
1369+
E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R
1370+
end if
13181371

13191372
H_L = (E_L + pres_L)/rho_L
13201373
H_R = (E_R + pres_R)/rho_R

0 commit comments

Comments
 (0)