Skip to content

Commit 6e8ded8

Browse files
hyeoksu-leeHyeoksu LeeHyeoksu LeeHyeoksu LeeHyeoksu Lee
authored
Low mach number correction for HLLC Riemann solver (#538)
Co-authored-by: Hyeoksu Lee <[email protected]> Co-authored-by: Hyeoksu Lee <[email protected]> Co-authored-by: Hyeoksu Lee <[email protected]> Co-authored-by: Hyeoksu Lee <[email protected]> Co-authored-by: Hyeoksu Lee <[email protected]> Co-authored-by: Hyeoksu Lee <[email protected]>
1 parent 3d0aff8 commit 6e8ded8

Some content is hidden

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

48 files changed

+3677
-18
lines changed

‎docs/documentation/case.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -368,6 +368,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
368368
| `null_weights` | Logical | Null WENO weights at boundaries |
369369
| `mp_weno` | Logical | Monotonicity preserving WENO |
370370
| `riemann_solver` | Integer | Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact* |
371+
| `low_Mach` | Integer | Low Mach number correction for HLLC Riemann solver: [0] None; [1] Pressure (Chen et al. 2022); [2] Velocity (Thornber et al. 2008) |
371372
| `avg_state` | Integer | Averaged state evaluation method: [1] Roe averagen*; [2] Arithmetic mean |
372373
| `wave_speeds` | Integer | Wave-speed estimation: [1] Direct (Batten et al. 1997); [2] Pressure-velocity* (Toro 1999) |
373374
| `weno_Re_flux` | Logical | Compute velocity gradient using scaler divergence theorem |
@@ -438,6 +439,8 @@ Practically, `weno_eps` $<10^{-6}$ is used.
438439
- `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 3.
439440
`riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively ([Toro, 2013](references.md#Toro13)).
440441

442+
- `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#Chen22)) and the improved velocity reconstruction method ([Thornber et al., 2008](reference.md#Thornber08)). This feature requires `riemann_solver = 2` and `model_eqns = 2`.
443+
441444
- `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.
442445
`avg_state = 1` and `2` correspond to Roe- and arithmetic averages, respectively.
443446

‎docs/documentation/references.md‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,8 @@
1010

1111
- <a id="Bryngelson19">Bryngelson, S. H., Schmidmayer, K., Coralic, V., Meng, J. C., Maeda, K., and Colonius, T. (2019). Mfc: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver. arXiv preprint arXiv:1907.10512.</a>
1212

13+
- <a id="Chen22">Chen, S. S., Li, J. P., Li, Z., Yuan, W., & Gao, Z. H. (2022). Anti-dissipation pressure correction under low Mach numbers for Godunov-type schemes. Journal of Computational Physics, 456, 111027. </a>
14+
1315
- <a id="Childs12">Childs, H., Brugger, E., Whitlock, B., Meredith, J., Ahern, S., Pugmire, D., Biagas, K., Miller, M., Harrison, C., Weber, G. H., Krishnan, H., Fogal, T., Sanderson, A., Garth, C., Bethel, E. W., Camp, D., R¨ubel, O., Durant, M., Favre, J. M., and Navr´atil, P. (2012). VisIt: An End-User Tool For Visualizing and Analyzing Very Large Data. In High Performance Visualization–Enabling Extreme-Scale Scientific Insight, pages 357–372.</a>
1416

1517
- <a id="Coralic15">Coralic, V. (2015). Simulation of shock-induced bubble collapse with application to vascular injury in shockwave lithotripsy. PhD thesis, California Institute of Technology.</a>
@@ -42,6 +44,8 @@
4244

4345
- <a id="Thompson90">Thompson, K. W. (1990). Time-dependent boundary conditions for hyperbolic systems, ii. Journal of computational physics, 89(2):439–461.</a>
4446

47+
- <a id="Thornber08">Thornber, B., Mosedale, A., Drikakis, D., Youngs, D., & Williams, R. J. (2008). An improved reconstruction method for compressible flows with low Mach number features. Journal of computational Physics, 227(10), 4873-4894.</a>
48+
4549
- <a id="Titarev04">Titarev, V. A. and Toro, E. F. (2004). Finite-volume weno schemes for three-dimensional conservation laws. Journal of Computational Physics, 201(1):238–260.</a>
4650

4751
- <a id="Tiwari13">Tiwari, A., Freund, J. B., and Pantano, C. (2013). A diffuse interface model with immiscibility preservation. Journal of computational physics, 252:290–309.</a>

‎examples/2D_GreshoVortex/case.py‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
#!/usr/bin/env python3
2+
3+
import math
4+
import json
5+
6+
gam = 1.4
7+
Ma0 = 1E-3
8+
c0 = math.sqrt(gam)
9+
u0 = Ma0*c0
10+
11+
Nx = 256
12+
Ny = 256
13+
Nz = 0
14+
dx = 1./Nx
15+
16+
cfl = 0.5
17+
T = 1./Ma0
18+
dt = cfl*dx/(u0+c0)
19+
Ntfinal = int(T/dt)
20+
Ntstart = 0
21+
Nfiles = 10
22+
t_save = int(math.ceil((Ntfinal-0)/float(Nfiles)))
23+
Nt = t_save*Nfiles
24+
t_step_start = Ntstart
25+
t_step_stop = int(Nt)
26+
27+
# Configuring case dictionary
28+
print(json.dumps({
29+
# Logistics ================================================================
30+
'run_time_info' : 'T',
31+
# ==========================================================================
32+
33+
# Computational Domain Parameters ==========================================
34+
'x_domain%beg' : 0.,
35+
'x_domain%end' : 1,
36+
'y_domain%beg' : 0.,
37+
'y_domain%end' : 1.,
38+
'm' : Nx,
39+
'n' : Ny,
40+
'p' : Nz,
41+
'dt' : dt,
42+
't_step_start' : t_step_start,
43+
't_step_stop' : t_step_stop,
44+
't_step_save' : t_save,
45+
# ==========================================================================
46+
47+
# Simulation Algorithm Parameters ==========================================
48+
'num_patches' : 1,
49+
'model_eqns' : 2,
50+
'num_fluids' : 1,
51+
'time_stepper' : 3,
52+
'weno_order' : 5,
53+
'weno_eps' : 1.E-16,
54+
# 'mapped_weno' : 'T',
55+
'wenoz' : 'T',
56+
'riemann_solver' : 2,
57+
'low_Mach' : 1,
58+
'wave_speeds' : 1,
59+
'avg_state' : 2,
60+
'bc_x%beg' : -1,
61+
'bc_x%end' : -1,
62+
'bc_y%beg' : -1,
63+
'bc_y%end' : -1,
64+
# ==========================================================================
65+
66+
# Formatted Database Files Structure Parameters ============================
67+
'format' : 1,
68+
'precision' : 2,
69+
'cons_vars_wrt' :'T',
70+
'prim_vars_wrt' :'T',
71+
'parallel_io' :'T',
72+
'fd_order' : 1,
73+
'omega_wrt(3)' :'T',
74+
# ==========================================================================
75+
76+
# Patch 1 ==================================================================
77+
'patch_icpp(1)%geometry' : 7,
78+
'patch_icpp(1)%hcid' : 203,
79+
'patch_icpp(1)%x_centroid' : 0.5,
80+
'patch_icpp(1)%y_centroid' : 0.5,
81+
'patch_icpp(1)%length_x' : 1,
82+
'patch_icpp(1)%length_y' : 1,
83+
'patch_icpp(1)%alpha_rho(1)' : 1.,
84+
'patch_icpp(1)%alpha(1)' : 1.,
85+
'patch_icpp(1)%vel(1)' : u0,
86+
'patch_icpp(1)%vel(2)' : Ma0,
87+
'patch_icpp(1)%pres' : 1.,
88+
# ==========================================================================
89+
90+
# Fluids Physical Parameters ===============================================
91+
'fluid_pp(1)%gamma' : 1.E+00/(gam-1.E+00),
92+
'fluid_pp(1)%pi_inf' : 0.,
93+
# =========================================================================
94+
}))
95+
96+
# ==============================================================================

‎src/simulation/include/inline_riemann.fpp‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,3 +44,21 @@
4444

4545
#:enddef compute_average_state
4646

47+
#:def compute_low_Mach_correction()
48+
49+
zcoef = min(1d0, max(vel_L_rms**5d-1/c_L, vel_R_rms**5d-1/c_R))
50+
pcorr = 0d0
51+
52+
if (low_Mach == 1) then
53+
pcorr = rho_L*rho_R* &
54+
(s_L - vel_L(dir_idx(1)))*(s_R - vel_R(dir_idx(1)))*(vel_R(dir_idx(1)) - vel_L(dir_idx(1)))/ &
55+
(rho_R*(s_R - vel_R(dir_idx(1))) - rho_L*(s_L - vel_L(dir_idx(1))))* &
56+
(zcoef - 1d0)
57+
else if (low_Mach == 2) then
58+
vel_L_tmp = 5d-1*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_L(dir_idx(1)) - vel_R(dir_idx(1))))
59+
vel_R_tmp = 5d-1*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_R(dir_idx(1)) - vel_L(dir_idx(1))))
60+
vel_L(dir_idx(1)) = vel_L_tmp
61+
vel_R(dir_idx(1)) = vel_R_tmp
62+
end if
63+
64+
#:enddef compute_low_Mach_correction

‎src/simulation/m_checker.fpp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,14 @@ contains
130130
elseif (riemann_solver /= 3 .and. avg_state == dflt_int) then
131131
call s_mpi_abort('avg_state must be set if '// &
132132
'riemann_solver != 3. Exiting ...')
133+
elseif (all(low_Mach /= (/0, 1, 2/))) then
134+
call s_mpi_abort('low_Mach must be 0, 1 or 2. Exiting ...')
135+
elseif (riemann_solver /= 2 .and. low_Mach /= 0) then
136+
call s_mpi_abort('low_Mach = 1 or 2 '// &
137+
'requires riemann_solver = 2. Exiting ...')
138+
elseif (low_Mach /= 0 .and. model_eqns /= 2) then
139+
call s_mpi_abort('low_Mach = 1 or 2 requires '// &
140+
'model_eqns = 2. Exiting ...')
133141
end if
134142
end subroutine s_check_inputs_riemann_solver
135143

‎src/simulation/m_global_parameters.fpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -141,6 +141,7 @@ module m_global_parameters
141141
logical :: weno_avg ! Average left/right cell-boundary states
142142
logical :: weno_Re_flux !< WENO reconstruct velocity gradients for viscous stress tensor
143143
integer :: riemann_solver !< Riemann solver algorithm
144+
integer :: low_Mach !< Low Mach number fix to HLLC Riemann solver
144145
integer :: wave_speeds !< Wave speeds estimation method
145146
integer :: avg_state !< Average state evaluation method
146147
logical :: alt_soundspeed !< Alternate mixture sound speed
@@ -166,7 +167,7 @@ module m_global_parameters
166167
!$acc declare create(num_dims, weno_polyn, weno_order, num_fluids, wenojs, mapped_weno, wenoz, teno)
167168
#:endif
168169

169-
!$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity)
170+
!$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach)
170171

171172
logical :: relax !< activate phase change
172173
integer :: relax_model !< Relaxation model
@@ -502,6 +503,7 @@ contains
502503
weno_avg = .false.
503504
weno_Re_flux = .false.
504505
riemann_solver = dflt_int
506+
low_Mach = 0
505507
wave_speeds = dflt_int
506508
avg_state = dflt_int
507509
alt_soundspeed = .false.
@@ -1052,7 +1054,7 @@ contains
10521054
!$acc update device(m, n, p)
10531055

10541056
!$acc update device(alt_soundspeed, acoustic_source, num_source)
1055-
!$acc update device(dt, sys_size, buff_size, pref, rhoref, gamma_idx, pi_inf_idx, E_idx, alf_idx, stress_idx, mpp_lim, bubbles, hypoelasticity, alt_soundspeed, avg_state, num_fluids, model_eqns, num_dims, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT)
1057+
!$acc update device(dt, sys_size, buff_size, pref, rhoref, gamma_idx, pi_inf_idx, E_idx, alf_idx, stress_idx, mpp_lim, bubbles, hypoelasticity, alt_soundspeed, avg_state, num_fluids, model_eqns, num_dims, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT, low_Mach)
10561058

10571059
#:if not MFC_CASE_OPTIMIZATION
10581060
!$acc update device(wenojs, mapped_weno, wenoz, teno)

‎src/simulation/m_mpi_proxy.fpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -188,7 +188,7 @@ contains
188188
189189
#:for VAR in ['t_step_old', 'm', 'n', 'p', 'm_glb', 'n_glb', 'p_glb', &
190190
& 't_step_start','t_step_stop','t_step_save','t_step_print', &
191-
& 'model_eqns','time_stepper', 'riemann_solver', &
191+
& 'model_eqns','time_stepper', 'riemann_solver', 'low_Mach', &
192192
& 'wave_speeds', 'avg_state', 'precision', 'bc_x%beg', 'bc_x%end', &
193193
& 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'fd_order', &
194194
& 'num_probes', 'num_integrals', 'bubble_model', 'thermal', &

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 35 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -878,10 +878,12 @@ contains
878878
real(kind(0d0)) :: R3V2Lbar, R3V2Rbar
879879
880880
real(kind(0d0)) :: vel_L_rms, vel_R_rms, vel_avg_rms
881+
real(kind(0d0)) :: vel_L_tmp, vel_R_tmp
881882
real(kind(0d0)) :: blkmod1, blkmod2
882883
real(kind(0d0)) :: rho_Star, E_Star, p_Star, p_K_Star
883884
real(kind(0d0)) :: pres_SL, pres_SR, Ms_L, Ms_R
884885
real(kind(0d0)) :: start, finish
886+
real(kind(0d0)) :: zcoef, pcorr !< low Mach number correction
885887
integer :: i, j, k, l, q !< Generic loop iterators
886888
integer :: idx1, idxi
887889
@@ -1510,10 +1512,10 @@ contains
15101512
end do
15111513
end do
15121514
end do
1513-
1515+
!$acc end parallel loop
15141516
elseif (model_eqns == 2 .and. bubbles) then
15151517
!$acc parallel loop collapse(3) gang vector default(present) private(R0_L, R0_R, V0_L, V0_R, P0_L, P0_R, pbw_L, pbw_R, vel_L, vel_R, &
1516-
!$acc rho_avg, alpha_L, alpha_R, h_avg, gamma_avg, s_L, s_R, s_S, nbub_L, nbub_R, ptilde_L, ptilde_R, vel_avg_rms, Re_L, Re_R)
1518+
!$acc rho_avg, alpha_L, alpha_R, h_avg, gamma_avg, s_L, s_R, s_S, nbub_L, nbub_R, ptilde_L, ptilde_R, vel_avg_rms, Re_L, Re_R, pcorr, zcoef, vel_L_tmp, vel_R_tmp)
15171519
do l = is3%beg, is3%end
15181520
do k = is2%beg, is2%end
15191521
do j = is1%beg, is1%end
@@ -1764,6 +1766,10 @@ contains
17641766
end do
17651767
end if
17661768
1769+
if (low_Mach == 2) then
1770+
@:compute_low_Mach_correction()
1771+
end if
1772+
17671773
if (wave_speeds == 1) then
17681774
s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R)
17691775
s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L)
@@ -1810,6 +1816,12 @@ contains
18101816
xi_M = (5d-1 + sign(5d-1, s_S))
18111817
xi_P = (5d-1 - sign(5d-1, s_S))
18121818
1819+
if (low_Mach == 1) then
1820+
@:compute_low_Mach_correction()
1821+
else
1822+
pcorr = 0d0
1823+
end if
1824+
18131825
!$acc loop seq
18141826
do i = 1, contxe
18151827
flux_rs${XYZ}$_vf(j, k, l, i) = &
@@ -1843,8 +1855,8 @@ contains
18431855
s_P*(xi_R*(dir_flg(dir_idx(i))*s_S + &
18441856
(1d0 - dir_flg(dir_idx(i)))* &
18451857
vel_R(dir_idx(i))) - vel_R(dir_idx(i)))) + &
1846-
dir_flg(dir_idx(i))*(pres_R - ptilde_R))
1847-
! if (j==0) print*, 'flux_rs_vf', flux_rs_vf(cont_idx%end+dir_idx(i))%sf(j,k,l)
1858+
dir_flg(dir_idx(i))*(pres_R - ptilde_R)) &
1859+
+ (s_M/s_L)*(s_P/s_R)*dir_flg(dir_idx(i))*pcorr
18481860
end do
18491861
18501862
! Energy flux.
@@ -1858,7 +1870,8 @@ contains
18581870
+ xi_P*(vel_R(dir_idx(1))*(E_R + pres_R - ptilde_R) + &
18591871
s_P*(xi_R*(E_R + (s_S - vel_R(dir_idx(1)))* &
18601872
(rho_R*s_S + (pres_R - ptilde_R)/ &
1861-
(s_R - vel_R(dir_idx(1))))) - E_R))
1873+
(s_R - vel_R(dir_idx(1))))) - E_R)) &
1874+
+ (s_M/s_L)*(s_P/s_R)*pcorr*s_S
18621875
18631876
! Volume fraction flux
18641877
@@ -1882,7 +1895,7 @@ contains
18821895
dir_flg(dir_idx(i))* &
18831896
s_P*(xi_R - 1d0))
18841897
1885-
!IF ( (model_eqns == 4) .or. (num_fluids==1) ) vel_src_rs_vf(dir_idx(i))%sf(j,k,l) = 0d0
1898+
!IF ( (model_eqns == 4) .or. (num_fluids==1) ) vel_src_rs_vf(idxi)%sf(j,k,l) = 0d0
18861899
end do
18871900
18881901
flux_src_rs${XYZ}$_vf(j, k, l, advxb) = vel_src_rs${XYZ}$_vf(j, k, l, dir_idx(1))
@@ -1968,7 +1981,7 @@ contains
19681981
!$acc end parallel loop
19691982
else
19701983
!$acc parallel loop collapse(3) gang vector default(present) private(vel_L, vel_R, Re_L, Re_R, &
1971-
!$acc rho_avg, h_avg, gamma_avg, alpha_L, alpha_R, s_L, s_R, s_S, vel_avg_rms) copyin(is1,is2,is3)
1984+
!$acc rho_avg, h_avg, gamma_avg, alpha_L, alpha_R, s_L, s_R, s_S, vel_avg_rms, pcorr, zcoef, vel_L_tmp, vel_R_tmp) copyin(is1,is2,is3)
19721985
do l = is3%beg, is3%end
19731986
do k = is2%beg, is2%end
19741987
do j = is1%beg, is1%end
@@ -2108,6 +2121,10 @@ contains
21082121
end do
21092122
end if
21102123
2124+
if (low_Mach == 2) then
2125+
@:compute_low_Mach_correction()
2126+
end if
2127+
21112128
if (wave_speeds == 1) then
21122129
s_L = min(vel_L(idx1) - c_L, vel_R(idx1) - c_R)
21132130
s_R = max(vel_R(idx1) + c_R, vel_L(idx1) + c_L)
@@ -2155,6 +2172,12 @@ contains
21552172
xi_M = (5d-1 + sign(5d-1, s_S))
21562173
xi_P = (5d-1 - sign(5d-1, s_S))
21572174
2175+
if (low_Mach == 1) then
2176+
@:compute_low_Mach_correction()
2177+
else
2178+
pcorr = 0d0
2179+
end if
2180+
21582181
!$acc loop seq
21592182
do i = 1, contxe
21602183
flux_rs${XYZ}$_vf(j, k, l, i) = &
@@ -2181,8 +2204,8 @@ contains
21812204
s_P*(xi_R*(dir_flg(idxi)*s_S + &
21822205
(1d0 - dir_flg(idxi))* &
21832206
vel_R(idxi)) - vel_R(idxi))) + &
2184-
dir_flg(idxi)*(pres_R))
2185-
! if (j==0) print*, 'flux_rs_vf', flux_rs_vf(cont_idx%end+dir_idx(i))%sf(j,k,l)
2207+
dir_flg(idxi)*(pres_R)) &
2208+
+ (s_M/s_L)*(s_P/s_R)*dir_flg(idxi)*pcorr
21862209
end do
21872210
21882211
! Energy flux.
@@ -2195,7 +2218,8 @@ contains
21952218
+ xi_P*(vel_R(idx1)*(E_R + pres_R) + &
21962219
s_P*(xi_R*(E_R + (s_S - vel_R(idx1))* &
21972220
(rho_R*s_S + pres_R/ &
2198-
(s_R - vel_R(idx1)))) - E_R))
2221+
(s_R - vel_R(idx1)))) - E_R)) &
2222+
+ (s_M/s_L)*(s_P/s_R)*pcorr*s_S
21992223
22002224
! Volume fraction flux
22012225
!$acc loop seq
@@ -2278,6 +2302,7 @@ contains
22782302
end do
22792303
end do
22802304
end do
2305+
!$acc end parallel loop
22812306
end if
22822307
end if
22832308
#:endfor

0 commit comments

Comments
 (0)