Low mach number correction for HLLC Riemann solver (MFlowCode#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]>
anandrdbz · Jul 29, 2024 · 6e8ded8 · 6e8ded8
1 parent 3d0aff8
commit 6e8ded8
Show file tree

Hide file tree

Showing 48 changed files with 3,677 additions and 18 deletions.
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -368,6 +368,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `null_weights`         | Logical | Null WENO weights at boundaries |
 | `mp_weno`              | Logical | Monotonicity preserving WENO |
 | `riemann_solver`       | Integer | Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*	 |
+| `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)	 |
 | `avg_state`	         | Integer | Averaged state evaluation method: [1] Roe averagen*; [2] Arithmetic mean  |
 | `wave_speeds`          | Integer | Wave-speed estimation: [1] Direct (Batten et al. 1997); [2] Pressure-velocity* (Toro 1999)	 |
 | `weno_Re_flux`         | Logical | Compute velocity gradient using scaler divergence theorem	 |
@@ -438,6 +439,8 @@ Practically, `weno_eps` $<10^{-6}$ is used.
 - `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 3.
 `riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively ([Toro, 2013](references.md#Toro13)).
 
+- `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`.
+
 - `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.
 `avg_state = 1` and `2` correspond to Roe- and arithmetic averages, respectively.
 

diff --git a/docs/documentation/references.md b/docs/documentation/references.md
@@ -10,6 +10,8 @@
 
 - <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>
 
+- <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>
+
 - <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>
 
 - <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 @@
 
 - <a id="Thompson90">Thompson, K. W. (1990). Time-dependent boundary conditions for hyperbolic systems, ii. Journal of computational physics, 89(2):439–461.</a>
 
+- <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>
+
 - <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>
 
 - <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>

diff --git a/examples/2D_GreshoVortex/case.py b/examples/2D_GreshoVortex/case.py
@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+
+import math
+import json
+
+gam = 1.4
+Ma0 = 1E-3
+c0 = math.sqrt(gam)
+u0 = Ma0*c0
+
+Nx = 256
+Ny = 256
+Nz = 0
+dx = 1./Nx
+
+cfl = 0.5
+T = 1./Ma0
+dt = cfl*dx/(u0+c0)
+Ntfinal = int(T/dt)
+Ntstart = 0
+Nfiles  = 10
+t_save  = int(math.ceil((Ntfinal-0)/float(Nfiles)))
+Nt      = t_save*Nfiles
+t_step_start    = Ntstart
+t_step_stop     = int(Nt)
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================================
+    'run_time_info'                 : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                  :  0.,
+    'x_domain%end'                  :  1,
+    'y_domain%beg'                  :  0.,
+    'y_domain%end'                  :  1.,
+    'm'                             : Nx,
+    'n'                             : Ny,
+    'p'                             : Nz,
+    'dt'                            : dt,
+    't_step_start'                  : t_step_start,
+    't_step_stop'                   : t_step_stop,
+    't_step_save'                   : t_save,
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'num_patches'                   : 1,
+    'model_eqns'                    : 2,
+    'num_fluids'                    : 1,
+    'time_stepper'                  : 3,
+    'weno_order'                    : 5,
+    'weno_eps'                      : 1.E-16,
+    # 'mapped_weno'                   : 'T',
+    'wenoz'                         : 'T',
+    'riemann_solver'                : 2,
+    'low_Mach'                      : 1,
+    'wave_speeds'                   : 1,
+    'avg_state'                     : 2,
+    'bc_x%beg'                      : -1,
+    'bc_x%end'                      : -1,
+    'bc_y%beg'                      : -1,
+    'bc_y%end'                      : -1,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                        : 1,
+    'precision'                     : 2,
+    'cons_vars_wrt'                 :'T',
+    'prim_vars_wrt'                 :'T',
+    'parallel_io'                   :'T',
+    'fd_order'                      : 1,
+    'omega_wrt(3)'                  :'T',
+    # ==========================================================================
+
+    # Patch 1 ==================================================================
+    'patch_icpp(1)%geometry'        : 7,
+    'patch_icpp(1)%hcid'            : 203,
+    'patch_icpp(1)%x_centroid'      : 0.5,
+    'patch_icpp(1)%y_centroid'      : 0.5,
+    'patch_icpp(1)%length_x'        : 1,
+    'patch_icpp(1)%length_y'        : 1,
+    'patch_icpp(1)%alpha_rho(1)'    : 1.,
+    'patch_icpp(1)%alpha(1)'        : 1.,
+    'patch_icpp(1)%vel(1)'          : u0,
+    'patch_icpp(1)%vel(2)'          : Ma0,
+    'patch_icpp(1)%pres'            : 1.,
+    # ==========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.E+00/(gam-1.E+00),
+    'fluid_pp(1)%pi_inf'           : 0.,
+    # =========================================================================
+}))
+
+# ==============================================================================
diff --git a/src/simulation/include/inline_riemann.fpp b/src/simulation/include/inline_riemann.fpp
@@ -44,3 +44,21 @@
 
 #:enddef compute_average_state
 
+#:def compute_low_Mach_correction()
+
+    zcoef = min(1d0, max(vel_L_rms**5d-1/c_L, vel_R_rms**5d-1/c_R))
+    pcorr = 0d0
+
+    if (low_Mach == 1) then
+        pcorr = rho_L*rho_R* &
+                (s_L - vel_L(dir_idx(1)))*(s_R - vel_R(dir_idx(1)))*(vel_R(dir_idx(1)) - vel_L(dir_idx(1)))/ &
+                (rho_R*(s_R - vel_R(dir_idx(1))) - rho_L*(s_L - vel_L(dir_idx(1))))* &
+                (zcoef - 1d0)
+    else if (low_Mach == 2) then
+        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))))
+        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))))
+        vel_L(dir_idx(1)) = vel_L_tmp
+        vel_R(dir_idx(1)) = vel_R_tmp
+    end if
+
+#:enddef compute_low_Mach_correction
diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
@@ -130,6 +130,14 @@ contains
         elseif (riemann_solver /= 3 .and. avg_state == dflt_int) then
             call s_mpi_abort('avg_state must be set if '// &
                              'riemann_solver != 3. Exiting ...')
+        elseif (all(low_Mach /= (/0, 1, 2/))) then
+            call s_mpi_abort('low_Mach must be 0, 1 or 2. Exiting ...')
+        elseif (riemann_solver /= 2 .and. low_Mach /= 0) then
+            call s_mpi_abort('low_Mach = 1 or 2 '// &
+                             'requires riemann_solver = 2. Exiting ...')
+        elseif (low_Mach /= 0 .and. model_eqns /= 2) then
+            call s_mpi_abort('low_Mach = 1 or 2 requires '// &
+                             'model_eqns = 2. Exiting ...')
         end if
     end subroutine s_check_inputs_riemann_solver
 

diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
@@ -141,6 +141,7 @@ module m_global_parameters
     logical :: weno_avg       ! Average left/right cell-boundary states
     logical :: weno_Re_flux   !< WENO reconstruct velocity gradients for viscous stress tensor
     integer :: riemann_solver !< Riemann solver algorithm
+    integer :: low_Mach       !< Low Mach number fix to HLLC Riemann solver
     integer :: wave_speeds    !< Wave speeds estimation method
     integer :: avg_state      !< Average state evaluation method
     logical :: alt_soundspeed !< Alternate mixture sound speed
@@ -166,7 +167,7 @@ module m_global_parameters
         !$acc declare create(num_dims, weno_polyn, weno_order, num_fluids, wenojs, mapped_weno, wenoz, teno)
     #:endif
 
-    !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity)
+    !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach)
 
     logical :: relax          !< activate phase change
     integer :: relax_model    !< Relaxation model
@@ -502,6 +503,7 @@ contains
         weno_avg = .false.
         weno_Re_flux = .false.
         riemann_solver = dflt_int
+        low_Mach = 0
         wave_speeds = dflt_int
         avg_state = dflt_int
         alt_soundspeed = .false.
@@ -1052,7 +1054,7 @@ contains
         !$acc update device(m, n, p)
 
         !$acc update device(alt_soundspeed, acoustic_source, num_source)
-        !$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)
+        !$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)
 
         #:if not MFC_CASE_OPTIMIZATION
             !$acc update device(wenojs, mapped_weno, wenoz, teno)

diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
@@ -188,7 +188,7 @@ contains
 
         #:for VAR in ['t_step_old', 'm', 'n', 'p', 'm_glb', 'n_glb', 'p_glb',  &
             & 't_step_start','t_step_stop','t_step_save','t_step_print',       &
-            & 'model_eqns','time_stepper', 'riemann_solver',                   &
+            & 'model_eqns','time_stepper', 'riemann_solver', 'low_Mach',       &
             & 'wave_speeds', 'avg_state', 'precision', 'bc_x%beg', 'bc_x%end', &
             & 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end',  'fd_order',     &
             & 'num_probes', 'num_integrals', 'bubble_model', 'thermal',        &

diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
@@ -878,10 +878,12 @@ contains
         real(kind(0d0)) :: R3V2Lbar, R3V2Rbar
 
         real(kind(0d0)) :: vel_L_rms, vel_R_rms, vel_avg_rms
+        real(kind(0d0)) :: vel_L_tmp, vel_R_tmp
         real(kind(0d0)) :: blkmod1, blkmod2
         real(kind(0d0)) :: rho_Star, E_Star, p_Star, p_K_Star
         real(kind(0d0)) :: pres_SL, pres_SR, Ms_L, Ms_R
         real(kind(0d0)) :: start, finish
+        real(kind(0d0)) :: zcoef, pcorr !< low Mach number correction
         integer :: i, j, k, l, q !< Generic loop iterators
         integer :: idx1, idxi
 
@@ -1510,10 +1512,10 @@ contains
                             end do
                         end do
                     end do
-
+                    !$acc end parallel loop
                 elseif (model_eqns == 2 .and. bubbles) then
                     !$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, &
-                    !$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)
+                    !$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)
                     do l = is3%beg, is3%end
                         do k = is2%beg, is2%end
                             do j = is1%beg, is1%end
@@ -1764,6 +1766,10 @@ contains
                                     end do
                                 end if
 
+                                if (low_Mach == 2) then
+                                    @:compute_low_Mach_correction()
+                                end if
+
                                 if (wave_speeds == 1) then
                                     s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R)
                                     s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L)
@@ -1810,6 +1816,12 @@ contains
                                 xi_M = (5d-1 + sign(5d-1, s_S))
                                 xi_P = (5d-1 - sign(5d-1, s_S))
 
+                                if (low_Mach == 1) then
+                                    @:compute_low_Mach_correction()
+                                else
+                                    pcorr = 0d0
+                                end if
+
                                 !$acc loop seq
                                 do i = 1, contxe
                                     flux_rs${XYZ}$_vf(j, k, l, i) = &
@@ -1843,8 +1855,8 @@ contains
                                                        s_P*(xi_R*(dir_flg(dir_idx(i))*s_S + &
                                                                   (1d0 - dir_flg(dir_idx(i)))* &
                                                                   vel_R(dir_idx(i))) - vel_R(dir_idx(i)))) + &
-                                                dir_flg(dir_idx(i))*(pres_R - ptilde_R))
-                                    ! if (j==0) print*, 'flux_rs_vf', flux_rs_vf(cont_idx%end+dir_idx(i))%sf(j,k,l)
+                                                dir_flg(dir_idx(i))*(pres_R - ptilde_R)) &
+                                        + (s_M/s_L)*(s_P/s_R)*dir_flg(dir_idx(i))*pcorr
                                 end do
 
                                 ! Energy flux.
@@ -1858,7 +1870,8 @@ contains
                                     + xi_P*(vel_R(dir_idx(1))*(E_R + pres_R - ptilde_R) + &
                                             s_P*(xi_R*(E_R + (s_S - vel_R(dir_idx(1)))* &
                                                        (rho_R*s_S + (pres_R - ptilde_R)/ &
-                                                        (s_R - vel_R(dir_idx(1))))) - E_R))
+                                                        (s_R - vel_R(dir_idx(1))))) - E_R)) &
+                                    + (s_M/s_L)*(s_P/s_R)*pcorr*s_S
 
                                 ! Volume fraction flux
 
@@ -1882,7 +1895,7 @@ contains
                                                 dir_flg(dir_idx(i))* &
                                                 s_P*(xi_R - 1d0))
 
-                                    !IF ( (model_eqns == 4) .or. (num_fluids==1) ) vel_src_rs_vf(dir_idx(i))%sf(j,k,l) = 0d0
+                                    !IF ( (model_eqns == 4) .or. (num_fluids==1) ) vel_src_rs_vf(idxi)%sf(j,k,l) = 0d0
                                 end do
 
                                 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
                     !$acc end parallel loop
                 else
                     !$acc parallel loop collapse(3) gang vector default(present) private(vel_L, vel_R, Re_L, Re_R, &
-                    !$acc rho_avg, h_avg, gamma_avg, alpha_L, alpha_R, s_L, s_R, s_S, vel_avg_rms) copyin(is1,is2,is3)
+                    !$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)
                     do l = is3%beg, is3%end
                         do k = is2%beg, is2%end
                             do j = is1%beg, is1%end
@@ -2108,6 +2121,10 @@ contains
                                     end do
                                 end if
 
+                                if (low_Mach == 2) then
+                                    @:compute_low_Mach_correction()
+                                end if
+
                                 if (wave_speeds == 1) then
                                     s_L = min(vel_L(idx1) - c_L, vel_R(idx1) - c_R)
                                     s_R = max(vel_R(idx1) + c_R, vel_L(idx1) + c_L)
@@ -2155,6 +2172,12 @@ contains
                                 xi_M = (5d-1 + sign(5d-1, s_S))
                                 xi_P = (5d-1 - sign(5d-1, s_S))
 
+                                if (low_Mach == 1) then
+                                    @:compute_low_Mach_correction()
+                                else
+                                    pcorr = 0d0
+                                end if
+
                                 !$acc loop seq
                                 do i = 1, contxe
                                     flux_rs${XYZ}$_vf(j, k, l, i) = &
@@ -2181,8 +2204,8 @@ contains
                                                        s_P*(xi_R*(dir_flg(idxi)*s_S + &
                                                                   (1d0 - dir_flg(idxi))* &
                                                                   vel_R(idxi)) - vel_R(idxi))) + &
-                                                dir_flg(idxi)*(pres_R))
-                                    ! if (j==0) print*, 'flux_rs_vf', flux_rs_vf(cont_idx%end+dir_idx(i))%sf(j,k,l)
+                                                dir_flg(idxi)*(pres_R)) &
+                                        + (s_M/s_L)*(s_P/s_R)*dir_flg(idxi)*pcorr
                                 end do
 
                                 ! Energy flux.
@@ -2195,7 +2218,8 @@ contains
                                     + xi_P*(vel_R(idx1)*(E_R + pres_R) + &
                                             s_P*(xi_R*(E_R + (s_S - vel_R(idx1))* &
                                                        (rho_R*s_S + pres_R/ &
-                                                        (s_R - vel_R(idx1)))) - E_R))
+                                                        (s_R - vel_R(idx1)))) - E_R)) &
+                                    + (s_M/s_L)*(s_P/s_R)*pcorr*s_S
 
                                 ! Volume fraction flux
                                 !$acc loop seq
@@ -2278,6 +2302,7 @@ contains
                             end do
                         end do
                     end do
+                    !$acc end parallel loop
                 end if
             end if
         #:endfor