Skip to content

Commit 6cfcb5e

Browse files
ChrisZYJChrisZYJ
committed
Checker & elastic energy macro
1 parent 7a28684 commit 6cfcb5e

6 files changed

Lines changed: 30 additions & 64 deletions

File tree

‎src/simulation/include/inline_riemann.fpp‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -91,3 +91,21 @@
9191
end if
9292
end if
9393
#:enddef compute_low_Mach_correction
94+
95+
#:def compute_hypo_elastic_energy(EL, ER, shear_cond)
96+
if (.not. hypo_energy_guard .or. ((G_L > verysmall) .and. (G_R > verysmall))) then
97+
${EL}$ = ${EL}$ + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
98+
${ER}$ = ${ER}$ + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
99+
if (${shear_cond}$) then
100+
${EL}$ = ${EL}$ + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
101+
${ER}$ = ${ER}$ + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
102+
end if
103+
end if
104+
#:enddef compute_hypo_elastic_energy
105+
106+
#:def compute_elastic_wave_speeds_lr()
107+
s_L = min(vel_L(dir_idx(1)) - sqrt(max(verysmall, c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)), &
108+
& vel_R(dir_idx(1)) - sqrt(max(verysmall, c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)))
109+
s_R = max(vel_R(dir_idx(1)) + sqrt(max(verysmall, c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)), &
110+
& vel_L(dir_idx(1)) + sqrt(max(verysmall, c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)))
111+
#:enddef compute_elastic_wave_speeds_lr

‎src/simulation/m_checker.fpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -138,6 +138,7 @@ contains
138138
& "HLLD hypoelasticity requires at least 2D (n must be > 0)")
139139
@:PROHIBIT(hypoelasticity .and. riemann_solver == 4 .and. num_fluids /= 2, &
140140
& "HLLD hypoelasticity currently requires exactly 2 fluid components")
141+
@:PROHIBIT(riemann_solver == 4 .and. mhd .and. hypoelasticity, "HLLD does not support combined MHD and hypoelasticity")
141142
@:PROHIBIT(hypoelasticity .and. (riemann_solver == 1 .or. riemann_solver == 2) .and. num_fluids > 2, &
142143
& "HLL/HLLC hypoelasticity supports at most 2 fluid components (blkmod uses fluid_pp(1:2) only)")
143144
@:PROHIBIT(hypoelasticity .and. riemann_solver == 4 .and. viscous, &

‎src/simulation/m_riemann_solver_hll.fpp‎

Lines changed: 2 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -344,16 +344,7 @@ contains
344344
do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
345345
tau_e_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%stress%beg - 1 + i)
346346
tau_e_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%stress%beg - 1 + i)
347-
! Elastic energy (guard skips when G near zero)
348-
if (.not. hypo_energy_guard .or. ((G_L > verysmall) .and. (G_R > verysmall))) then
349-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
350-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
351-
! Double for shear stresses
352-
if (any(eqn_idx%stress%beg - 1 + i == shear_indices)) then
353-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
354-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
355-
end if
356-
end if
347+
@:compute_hypo_elastic_energy(E_L, E_R, any(eqn_idx%stress%beg - 1 + i == shear_indices))
357348
end do
358349
end if
359350

@@ -394,14 +385,7 @@ contains
394385
s_R = max(vel_R(dir_idx(1)) + c_fast%R, vel_L(dir_idx(1)) + c_fast%L)
395386
else if (hypoelasticity) then
396387
! Elastic wave speed, Rodriguez et al. JCP (2019)
397-
s_L = min(vel_L(dir_idx(1)) - sqrt(max(verysmall, &
398-
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)), &
399-
& vel_R(dir_idx(1)) - sqrt(max(verysmall, &
400-
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)))
401-
s_R = max(vel_R(dir_idx(1)) + sqrt(max(verysmall, &
402-
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)), &
403-
& vel_L(dir_idx(1)) + sqrt(max(verysmall, &
404-
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)))
388+
@:compute_elastic_wave_speeds_lr()
405389
else if (hyperelasticity) then
406390
s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L), &
407391
& vel_R(dir_idx(1)) - sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R))

‎src/simulation/m_riemann_solver_hllc.fpp‎

Lines changed: 4 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -305,16 +305,7 @@ contains
305305
end do
306306
$:GPU_LOOP(parallelism='[seq]')
307307
do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
308-
! Elastic energy (guard skips when G near zero)
309-
if (.not. hypo_energy_guard .or. ((G_L > verysmall) .and. (G_R > verysmall))) then
310-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
311-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
312-
! Additional terms in 2D and 3D
313-
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
314-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
315-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
316-
end if
317-
end if
308+
@:compute_hypo_elastic_energy(E_L, E_R, any(eqn_idx%stress%beg - 1 + i == shear_indices))
318309
end do
319310
end if
320311

@@ -376,14 +367,7 @@ contains
376367
if (wave_speeds == wave_speeds_direct) then
377368
if (elasticity) then
378369
! Elastic wave speed, Rodriguez et al. JCP (2019)
379-
s_L = min(vel_L(dir_idx(1)) - sqrt(max(verysmall, &
380-
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)), &
381-
& vel_R(dir_idx(1)) - sqrt(max(verysmall, &
382-
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)))
383-
s_R = max(vel_R(dir_idx(1)) + sqrt(max(verysmall, &
384-
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)), &
385-
& vel_L(dir_idx(1)) + sqrt(max(verysmall, &
386-
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)))
370+
@:compute_elastic_wave_speeds_lr()
387371
s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + tau_e_L(dir_idx_tau(1)) &
388372
& + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1)) &
389373
& *(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R &
@@ -1413,16 +1397,7 @@ contains
14131397
end do
14141398
$:GPU_LOOP(parallelism='[seq]')
14151399
do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
1416-
! Elastic energy (guard skips when G near zero)
1417-
if (.not. hypo_energy_guard .or. ((G_L > verysmall) .and. (G_R > verysmall))) then
1418-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
1419-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
1420-
! Additional terms in 2D and 3D
1421-
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
1422-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
1423-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
1424-
end if
1425-
end if
1400+
@:compute_hypo_elastic_energy(E_L, E_R, any(eqn_idx%stress%beg - 1 + i == shear_indices))
14261401
end do
14271402
end if
14281403

@@ -1487,14 +1462,7 @@ contains
14871462
if (wave_speeds == wave_speeds_direct) then
14881463
if (elasticity) then
14891464
! Elastic wave speed, Rodriguez et al. JCP (2019)
1490-
s_L = min(vel_L(dir_idx(1)) - sqrt(max(verysmall, &
1491-
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)), &
1492-
& vel_R(dir_idx(1)) - sqrt(max(verysmall, &
1493-
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)))
1494-
s_R = max(vel_R(dir_idx(1)) + sqrt(max(verysmall, &
1495-
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)), &
1496-
& vel_L(dir_idx(1)) + sqrt(max(verysmall, &
1497-
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)))
1465+
@:compute_elastic_wave_speeds_lr()
14981466
s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + tau_e_L(dir_idx_tau(1)) &
14991467
& + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1)) &
15001468
& *(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R &

‎src/simulation/m_riemann_solver_hypo_hlld.fpp‎

Lines changed: 4 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
!> @brief Dual-pass HLLD approximate Riemann solver for hypoelastic flows, with non-conservative interface-velocity coupling
66
#:include 'case.fpp'
77
#:include 'macros.fpp'
8+
#:include 'inline_riemann.fpp'
89

910
module m_riemann_solver_hypo_hlld
1011

@@ -281,16 +282,9 @@ contains
281282

282283
$:GPU_LOOP(parallelism='[seq]')
283284
do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
284-
! Elastic energy (guard skips when G near zero)
285-
if (.not. hypo_energy_guard .or. ((G_L > verysmall) .and. (G_R > verysmall))) then
286-
E%L = E%L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
287-
E%R = E%R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
288-
! Shear terms doubled: 2D/2D-axisym i==2 only; 3D i==2,4,5
289-
if ((n > 0 .and. p == 0 .and. i == 2) .or. (p > 0 .and. (i == 2 .or. i == 4 .or. i == 5))) then
290-
E%L = E%L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
291-
E%R = E%R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
292-
end if
293-
end if
285+
@:compute_hypo_elastic_energy(E%L, E%R, &
286+
& (n > 0 .and. p == 0 .and. i == 2) .or. (p > 0 .and. (i == 2 &
287+
& .or. i == 4 .or. i == 5)))
294288
end do
295289

296290
H%L = (E%L + pres%L)/rho%L

‎toolchain/mfc/case_validator.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1074,6 +1074,7 @@ def check_mhd_simulation(self):
10741074
self.prohibit(mhd and wave_speeds is not None and wave_speeds == 2, "MHD requires wave_speeds = 1")
10751075
hypoelasticity = self.get("hypoelasticity", "F") == "T"
10761076
self.prohibit(riemann_solver == 4 and not mhd and not hypoelasticity, "HLLD (riemann_solver = 4) requires MHD or hypoelasticity")
1077+
self.prohibit(riemann_solver == 4 and mhd and hypoelasticity, "HLLD does not support combined MHD and hypoelasticity")
10771078
self.prohibit(riemann_solver == 4 and relativity, "HLLD is not available for RMHD (relativity)")
10781079
self.prohibit(hyper_cleaning and not mhd, "Hyperbolic cleaning requires mhd to be enabled")
10791080
self.prohibit(hyper_cleaning and n is not None and n == 0, "Hyperbolic cleaning is not supported for 1D simulations")

0 commit comments

Comments
 (0)