Skip to content

Commit f8e3912

Browse files
hyeoksu-leeHyeoksu LeeHyeoksu LeeHyeoksu Lee
authored
low_Mach for HLL Riemann solver and 6-equation model (#833)
Co-authored-by: Hyeoksu Lee <hyeoksu@carpenter02.hsn.ex4000.erdc.hpc.mil> Co-authored-by: Hyeoksu Lee <hyeoksu@carpenter04.hsn.ex4000.erdc.hpc.mil> Co-authored-by: Hyeoksu Lee <hyeoksu@carpenter08.hsn.ex4000.erdc.hpc.mil>
1 parent 72df181 commit f8e3912

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

54 files changed

+3660
-2165
lines changed

‎docs/documentation/case.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -461,7 +461,7 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
461461
`riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively ([Toro, 2013](references.md)).
462462
`riemann_solver = 4` is only for MHD simulations. It resolves 5 of the full seven-wave structure of the MHD equations ([Miyoshi and Kusano, 2005](references.md)).
463463

464-
- `low_Mach` specifies the choice of the low Mach number correction scheme for the HLLC Riemann solver. `low_Mach = 0` is default value and does not apply any correction scheme. `low_Mach = 1` and `2` apply the anti-dissipation pressure correction method ([Chen et al., 2022](references.md)) and the improved velocity reconstruction method ([Thornber et al., 2008](references.md)). This feature requires `riemann_solver = 2` and `model_eqns = 2`.
464+
- `low_Mach` specifies the choice of the low Mach number correction scheme for the HLLC Riemann solver. `low_Mach = 0` is default value and does not apply any correction scheme. `low_Mach = 1` and `2` apply the anti-dissipation pressure correction method ([Chen et al., 2022](references.md)) and the improved velocity reconstruction method ([Thornber et al., 2008](references.md)). This feature requires `model_eqns = 2` or `3`. `low_Mach = 1` works for `riemann_solver = 1` and `2`, but `low_Mach = 2` only works for `riemann_solver = 2`.
465465

466466
- `avg_state` specifies the choice of the method to compute averaged variables at the cell-boundaries from the left and the right states in the Riemann solver by an integer of 1 or 2.
467467
`avg_state = 1` and `2` correspond to Roe- and arithmetic averages, respectively.

‎src/post_process/m_start_up.f90

Lines changed: 8 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -544,31 +544,17 @@ subroutine s_save_data(t_step, varname, pres, c, H)
544544
end if
545545

546546
! Adding the vorticity to the formatted database file
547-
if (p > 0) then
548-
do i = 1, num_vels
549-
if (omega_wrt(i)) then
547+
do i = 1, 3
548+
if (omega_wrt(i)) then
550549

551-
call s_derive_vorticity_component(i, q_prim_vf, q_sf)
550+
call s_derive_vorticity_component(i, q_prim_vf, q_sf)
552551

553-
write (varname, '(A,I0)') 'omega', i
554-
call s_write_variable_to_formatted_database_file(varname, t_step)
555-
556-
varname(:) = ' '
557-
end if
558-
end do
559-
elseif (n > 0) then
560-
do i = 1, num_vels
561-
if (omega_wrt(i)) then
562-
563-
call s_derive_vorticity_component(i, q_prim_vf, q_sf)
564-
565-
write (varname, '(A,I0)') 'omega', i
566-
call s_write_variable_to_formatted_database_file(varname, t_step)
552+
write (varname, '(A,I0)') 'omega', i
553+
call s_write_variable_to_formatted_database_file(varname, t_step)
567554

568-
varname(:) = ' '
569-
end if
570-
end do
571-
end if
555+
varname(:) = ' '
556+
end if
557+
end do
572558

573559
if (ib) then
574560
q_sf = real(ib_markers%sf(-offset_x%beg:m + offset_x%end, -offset_y%beg:n + offset_y%end, -offset_z%beg:p + offset_z%end))

‎src/simulation/include/inline_riemann.fpp

Lines changed: 24 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -78,19 +78,30 @@
7878

7979
#:def compute_low_Mach_correction()
8080

81-
zcoef = min(1._wp, max(vel_L_rms**5e-1_wp/c_L, vel_R_rms**5e-1_wp/c_R))
82-
pcorr = 0._wp
83-
84-
if (low_Mach == 1) then
85-
pcorr = rho_L*rho_R* &
86-
(s_L - vel_L(dir_idx(1)))*(s_R - vel_R(dir_idx(1)))*(vel_R(dir_idx(1)) - vel_L(dir_idx(1)))/ &
87-
(rho_R*(s_R - vel_R(dir_idx(1))) - rho_L*(s_L - vel_L(dir_idx(1))))* &
88-
(zcoef - 1._wp)
89-
else if (low_Mach == 2) then
90-
vel_L_tmp = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_L(dir_idx(1)) - vel_R(dir_idx(1))))
91-
vel_R_tmp = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_R(dir_idx(1)) - vel_L(dir_idx(1))))
92-
vel_L(dir_idx(1)) = vel_L_tmp
93-
vel_R(dir_idx(1)) = vel_R_tmp
81+
if (riemann_solver == 1) then
82+
83+
zcoef = min(1._wp, max(vel_L_rms**5e-1_wp/c_L, vel_R_rms**5e-1_wp/c_R))
84+
pcorr = 0._wp
85+
86+
if (low_Mach == 1) then
87+
pcorr = -(s_P - s_M)*(rho_L + rho_R)/8._wp*(zcoef - 1._wp)
88+
end if
89+
90+
else if (riemann_solver == 2) then
91+
zcoef = min(1._wp, max(vel_L_rms**5e-1_wp/c_L, vel_R_rms**5e-1_wp/c_R))
92+
pcorr = 0._wp
93+
94+
if (low_Mach == 1) then
95+
pcorr = rho_L*rho_R* &
96+
(s_L - vel_L(dir_idx(1)))*(s_R - vel_R(dir_idx(1)))*(vel_R(dir_idx(1)) - vel_L(dir_idx(1)))/ &
97+
(rho_R*(s_R - vel_R(dir_idx(1))) - rho_L*(s_L - vel_L(dir_idx(1))))* &
98+
(zcoef - 1._wp)
99+
else if (low_Mach == 2) then
100+
vel_L_tmp = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_L(dir_idx(1)) - vel_R(dir_idx(1))))
101+
vel_R_tmp = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_R(dir_idx(1)) - vel_L(dir_idx(1))))
102+
vel_L(dir_idx(1)) = vel_L_tmp
103+
vel_R(dir_idx(1)) = vel_R_tmp
104+
end if
94105
end if
95106

96107
#:enddef compute_low_Mach_correction

‎src/simulation/m_checker.fpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -98,8 +98,8 @@ contains
9898
@:PROHIBIT(riemann_solver /= 3 .and. wave_speeds == dflt_int, "wave_speeds must be set if riemann_solver != 3")
9999
@:PROHIBIT(riemann_solver /= 3 .and. avg_state == dflt_int, "avg_state must be set if riemann_solver != 3")
100100
@:PROHIBIT(all(low_Mach /= (/0, 1, 2/)), "low_Mach must be 0, 1 or 2")
101-
@:PROHIBIT(riemann_solver /= 2 .and. low_Mach /= 0, "low_Mach = 1 or 2 requires riemann_solver = 2")
102-
@:PROHIBIT(low_Mach /= 0 .and. model_eqns /= 2, "low_Mach = 1 or 2 requires model_eqns = 2")
101+
@:PROHIBIT(riemann_solver /= 2 .and. low_Mach == 2, "low_Mach = 2 requires riemann_solver = 2")
102+
@:PROHIBIT(low_Mach /= 0 .and. all(model_eqns /= (/2, 3/)), "low_Mach = 1 or 2 requires model_eqns = 2 or 3")
103103
end subroutine s_check_inputs_riemann_solver
104104

105105
!> Checks constraints on geometry and precision

‎src/simulation/m_riemann_solvers.fpp

Lines changed: 46 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -325,8 +325,10 @@ contains
325325

326326
real(wp) :: ptilde_L, ptilde_R
327327
real(wp) :: vel_L_rms, vel_R_rms, vel_avg_rms
328+
real(wp) :: vel_L_tmp, vel_R_tmp
328329
real(wp) :: Ms_L, Ms_R, pres_SL, pres_SR
329330
real(wp) :: alpha_L_sum, alpha_R_sum
331+
real(wp) :: zcoef, pcorr !< low Mach number correction
330332

331333
type(riemann_states) :: c_fast, pres_mag
332334
type(riemann_states_vec3) :: B
@@ -367,7 +369,8 @@ contains
367369
!$acc xi_field_L, xi_field_R, &
368370
!$acc Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, &
369371
!$acc Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, &
370-
!$acc c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm)
372+
!$acc c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, &
373+
!$acc pcorr, zcoef, vel_L_tmp, vel_R_tmp)
371374
do l = is3%beg, is3%end
372375
do k = is2%beg, is2%end
373376
do j = is1%beg, is1%end
@@ -748,6 +751,13 @@ contains
748751
+ (5e-1_wp - sign(5e-1_wp, s_L)) &
749752
*(5e-1_wp + sign(5e-1_wp, s_R))
750753

754+
! Low Mach correction
755+
if (low_Mach == 1) then
756+
@:compute_low_Mach_correction()
757+
else
758+
pcorr = 0._wp
759+
end if
760+
751761
! Mass
752762
if (.not. relativity) then
753763
!$acc loop seq
@@ -852,7 +862,8 @@ contains
852862
+ dir_flg(dir_idx(i))*(pres_L - ptilde_L)) &
853863
+ s_M*s_P*(rho_L*vel_L(dir_idx(i)) &
854864
- rho_R*vel_R(dir_idx(i)))) &
855-
/(s_M - s_P)
865+
/(s_M - s_P) &
866+
+ (s_M/s_L)*(s_P/s_R)*pcorr*(vel_R(dir_idx(i)) - vel_L(dir_idx(i)))
856867
end do
857868
else if (hypoelasticity) then
858869
!$acc loop seq
@@ -882,7 +893,8 @@ contains
882893
+ dir_flg(dir_idx(i))*pres_L) &
883894
+ s_M*s_P*(rho_L*vel_L(dir_idx(i)) &
884895
- rho_R*vel_R(dir_idx(i)))) &
885-
/(s_M - s_P)
896+
/(s_M - s_P) &
897+
+ (s_M/s_L)*(s_P/s_R)*pcorr*(vel_R(dir_idx(i)) - vel_L(dir_idx(i)))
886898
end do
887899
end if
888900

@@ -907,7 +919,8 @@ contains
907919
(s_M*vel_R(dir_idx(1))*(E_R + pres_R - ptilde_R) &
908920
- s_P*vel_L(dir_idx(1))*(E_L + pres_L - ptilde_L) &
909921
+ s_M*s_P*(E_L - E_R)) &
910-
/(s_M - s_P)
922+
/(s_M - s_P) &
923+
+ (s_M/s_L)*(s_P/s_R)*pcorr*(vel_R_rms - vel_L_rms)/2._wp
911924
else if (hypoelasticity) then
912925
!TODO: simplify this so it's not split into 3
913926
if (num_dims == 1) then
@@ -946,7 +959,8 @@ contains
946959
(s_M*vel_R(dir_idx(1))*(E_R + pres_R) &
947960
- s_P*vel_L(dir_idx(1))*(E_L + pres_L) &
948961
+ s_M*s_P*(E_L - E_R)) &
949-
/(s_M - s_P)
962+
/(s_M - s_P) &
963+
+ (s_M/s_L)*(s_P/s_R)*pcorr*(vel_R_rms - vel_L_rms)/2._wp
950964
end if
951965
952966
! Elastic Stresses
@@ -1287,7 +1301,8 @@ contains
12871301
!$acc private(vel_L, vel_R, vel_K_Star, Re_L, Re_R, rho_avg, h_avg, gamma_avg, &
12881302
!$acc s_L, s_R, s_S, vel_avg_rms, alpha_L, alpha_R, Ys_L, Ys_R, Xs_L, Xs_R, &
12891303
!$acc Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, &
1290-
!$acc tau_e_L, tau_e_R, G_L, G_R, flux_ene_e, xi_field_L, xi_field_R)
1304+
!$acc tau_e_L, tau_e_R, G_L, G_R, flux_ene_e, xi_field_L, xi_field_R, pcorr, &
1305+
!$acc zcoef, vel_L_tmp, vel_R_tmp)
12911306
do l = is3%beg, is3%end
12921307
do k = is2%beg, is2%end
12931308
do j = is1%beg, is1%end
@@ -1476,6 +1491,11 @@ contains
14761491
end do
14771492
end if
14781493
1494+
! Low Mach correction
1495+
if (low_Mach == 2) then
1496+
@:compute_low_Mach_correction()
1497+
end if
1498+
14791499
! COMPUTING THE DIRECT WAVE SPEEDS
14801500
if (wave_speeds == 1) then
14811501
if (elasticity) then
@@ -1551,6 +1571,13 @@ contains
15511571
vel_K_Star = vel_L(idx1)*(1_wp - xi_MP) + xi_MP*vel_R(idx1) + &
15521572
xi_MP*xi_PP*(s_S - vel_R(idx1))
15531573
1574+
! Low Mach correction
1575+
if (low_Mach == 1) then
1576+
@:compute_low_Mach_correction()
1577+
else
1578+
pcorr = 0._wp
1579+
end if
1580+
15541581
! COMPUTING FLUXES
15551582
! MASS FLUX.
15561583
!$acc loop seq
@@ -1566,12 +1593,14 @@ contains
15661593
do i = 1, num_dims
15671594
idxi = dir_idx(i)
15681595
flux_rs${XYZ}$_vf(j, k, l, contxe + idxi) = rho_Star*vel_K_Star* &
1569-
(dir_flg(idxi)*vel_K_Star + (1_wp - dir_flg(idxi))*(xi_M*vel_L(idxi) + xi_P*vel_R(idxi))) + dir_flg(idxi)*p_Star
1596+
(dir_flg(idxi)*vel_K_Star + (1_wp - dir_flg(idxi))*(xi_M*vel_L(idxi) + xi_P*vel_R(idxi))) + dir_flg(idxi)*p_Star &
1597+
+ (s_M/s_L)*(s_P/s_R)*dir_flg(idxi)*pcorr
15701598
end do
15711599
15721600
! ENERGY FLUX.
15731601
! f = u*(E-\sigma), q = E, q_star = \xi*E+(s-u)(\rho s_star - \sigma/(s-u))
1574-
flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_star + p_Star)*vel_K_Star
1602+
flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_star + p_Star)*vel_K_Star &
1603+
+ (s_M/s_L)*(s_P/s_R)*pcorr*s_S
15751604
15761605
! ELASTICITY. Elastic shear stress additions for the momentum and energy flux
15771606
if (elasticity) then
@@ -1623,7 +1652,8 @@ contains
16231652
((xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1))* &
16241653
(gammas(i)*p_K_Star + pi_infs(i)) + &
16251654
(xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, i + contxb - 1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + contxb - 1))* &
1626-
qvs(i))*vel_K_Star
1655+
qvs(i))*vel_K_Star &
1656+
+ (s_M/s_L)*(s_P/s_R)*pcorr*s_S*(xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1))
16271657
end do
16281658
16291659
flux_src_rs${XYZ}$_vf(j, k, l, advxb) = vel_src_rs${XYZ}$_vf(j, k, l, idx1)
@@ -1703,6 +1733,7 @@ contains
17031733
do l = is3%beg, is3%end
17041734
do k = is2%beg, is2%end
17051735
do j = is1%beg, is1%end
1736+
17061737
!$acc loop seq
17071738
do i = 1, contxe
17081739
alpha_rho_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, i)
@@ -2200,6 +2231,7 @@ contains
22002231
end do
22012232
end if
22022233
2234+
! Low Mach correction
22032235
if (low_Mach == 2) then
22042236
@:compute_low_Mach_correction()
22052237
end if
@@ -2250,6 +2282,7 @@ contains
22502282
xi_M = (5e-1_wp + sign(5e-1_wp, s_S))
22512283
xi_P = (5e-1_wp - sign(5e-1_wp, s_S))
22522284
2285+
! Low Mach correction
22532286
if (low_Mach == 1) then
22542287
@:compute_low_Mach_correction()
22552288
else
@@ -2667,6 +2700,7 @@ contains
26672700
end do
26682701
end if
26692702
2703+
! Low Mach correction
26702704
if (low_Mach == 2) then
26712705
@:compute_low_Mach_correction()
26722706
end if
@@ -2727,14 +2761,15 @@ contains
27272761
xi_M = (5e-1_wp + sign(5e-1_wp, s_S))
27282762
xi_P = (5e-1_wp - sign(5e-1_wp, s_S))
27292763
2730-
! COMPUTING THE HLLC FLUXES
2731-
! MASS FLUX.
2764+
! Low Mach correction
27322765
if (low_Mach == 1) then
27332766
@:compute_low_Mach_correction()
27342767
else
27352768
pcorr = 0._wp
27362769
end if
27372770
2771+
! COMPUTING THE HLLC FLUXES
2772+
! MASS FLUX.
27382773
!$acc loop seq
27392774
do i = 1, contxe
27402775
flux_rs${XYZ}$_vf(j, k, l, i) = &

0 commit comments

Comments
 (0)