muscl_eps (#1383)

Co-authored-by: Spencer Bryngelson <sbryngelson@gmail.com>

Author: ChrisZYJ

docs/documentation/case.md

@@ -445,6 +445,7 @@ See @ref equations "Equations" for the mathematical models these parameters cont
 | `mp_weno`                  | Logical | Monotonicity preserving WENO |
 | `muscl_order`              | Integer | MUSCL order [1,2] |
 | `muscl_lim`                | Integer | MUSCL Slope Limiter: [1] minmod; [2] monotonized central; [3] Van Albada; [4] Van Leer; [5] SUPERBEE |
+| `muscl_eps`                | Real    | MUSCL limiter slope-product threshold (default: hard-coded thresholds; set to 0 for textbook behavior) |
 | `flux_lim`                 | Integer | Flux limiter for post-process: [1] minmod; [2] MUSCL; [3] OSPRE; [4] SUPERBEE |
 | `int_comp`                 | Logical | THINC Interface Compression |
 | `ic_eps`                   | Real    | Interface compression threshold (default: 1e-4) |
@@ -542,6 +543,10 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
 - `muscl_lim` specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5.
 `muscl_lim = 1`, `2`, `3`, `4`, and `5` correspond to minmod, monotonized central, Van Albada, Van Leer, and SUPERBEE, respectively.
 
+- `muscl_eps` controls the slope-product activation threshold for all MUSCL limiters.
+When not set (default), the threshold is 1e-9 for minmod/MC, and 1e-6 for others.
+Setting `muscl_eps = 0` gives textbook limiter behavior where limiters activate whenever both slopes have the same sign.
+
 - `int_comp` activates interface compression using THINC used in MUSCL Reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
 
 - `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 4.

src/simulation/m_global_parameters.fpp

@@ -148,6 +148,7 @@ module m_global_parameters
     logical :: nv_uvm_pref_gpu  !< Enable explicit gpu memory hints (default FALSE)
     !> @}
 
+    real(wp)           :: muscl_eps                    !< MUSCL limiter slope-product threshold
     real(wp)           :: weno_eps                     !< Binding for the WENO nonlinear weights
     real(wp)           :: teno_CT                      !< Smoothness threshold for TENO
     logical            :: mp_weno                      !< Monotonicity preserving (MP) WENO
@@ -197,6 +198,7 @@ module m_global_parameters
         $:GPU_DECLARE(create='[recon_type, muscl_order, muscl_polyn, muscl_lim]')
     #:endif
 
+    $:GPU_DECLARE(create='[muscl_eps]')
     $:GPU_DECLARE(create='[mpp_lim, model_eqns, mixture_err, alt_soundspeed]')
     $:GPU_DECLARE(create='[avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity]')
     $:GPU_DECLARE(create='[hyperelasticity, hyper_model, elasticity, low_Mach]')
@@ -521,6 +523,7 @@ contains
         model_eqns = dflt_int
         mpp_lim = .false.
         time_stepper = dflt_int
+        muscl_eps = dflt_real
         weno_eps = dflt_real
         teno_CT = dflt_real
         mp_weno = .false.
@@ -856,6 +859,15 @@ contains
             $:GPU_UPDATE(device='[igr, igr_order, igr_iter_solver]')
         #:endif
 
+        ! muscl_eps: use per-limiter defaults when user did not set it
+        if (f_is_default(muscl_eps)) then
+            if (muscl_lim <= 2) then
+                muscl_eps = 1e-9_wp  ! minmod, MC
+            else
+                muscl_eps = 1e-6_wp  ! Van Albada, Van Leer, SUPERBEE
+            end if
+        end if
+
         ! Initialize counts: viscous fluids, surface-tension interfaces, curvature interfaces
         Re_size = 0
         Re_size_max = 0
@@ -1220,6 +1232,7 @@ contains
             $:GPU_UPDATE(device='[num_fluids, num_dims, viscous, num_vels, nb, muscl_lim]')
         #:endif
 
+        $:GPU_UPDATE(device='[muscl_eps]')
         $:GPU_UPDATE(device='[dir_idx, dir_flg, dir_idx_tau]')
 
         $:GPU_UPDATE(device='[relax, relax_model, palpha_eps, ptgalpha_eps]')

src/simulation/m_muscl.fpp

@@ -170,27 +170,26 @@ contains
                                     slope = 0._wp
 
                                     if (muscl_lim == 1) then  ! minmod
-                                        if (slopeL*slopeR > 1e-9_wp) then
+                                        if (slopeL*slopeR > muscl_eps) then
                                             slope = min(abs(slopeL), abs(slopeR))
                                         end if
                                         if (slopeL < 0._wp) slope = -slope
                                     else if (muscl_lim == 2) then  ! MC
-                                        if (slopeL*slopeR > 1e-9_wp) then
+                                        if (slopeL*slopeR > muscl_eps) then
                                             slope = min(2._wp*abs(slopeL), 2._wp*abs(slopeR))
                                             slope = min(slope, 5e-1_wp*(abs(slopeL) + abs(slopeR)))
                                         end if
                                         if (slopeL < 0._wp) slope = -slope
                                     else if (muscl_lim == 3) then  ! Van Albada
-                                        if (abs(slopeL) > 1e-6_wp .and. abs(slopeR) > 1e-6_wp .and. abs(slopeL + slopeR) &
-                                            & > 1e-6_wp .and. slopeL*slopeR > 1e-6_wp) then
+                                        if (slopeL*slopeR > muscl_eps) then
                                             slope = ((slopeL + slopeR)*slopeL*slopeR)/(slopeL**2._wp + slopeR**2._wp)
                                         end if
                                     else if (muscl_lim == 4) then  ! Van Leer
-                                        if (abs(slopeL + slopeR) > 1.e-6_wp .and. slopeL*slopeR > 1.e-6_wp) then
+                                        if (slopeL*slopeR > muscl_eps) then
                                             slope = 2._wp*slopeL*slopeR/(slopeL + slopeR)
                                         end if
                                     else if (muscl_lim == 5) then  ! SUPERBEE
-                                        if (slopeL*slopeR > 1e-6_wp) then
+                                        if (slopeL*slopeR > muscl_eps) then
                                             slope = -1._wp*min(-min(2._wp*abs(slopeL), abs(slopeR)), -min(abs(slopeL), &
                                                                & 2._wp*abs(slopeR)))
                                         end if

src/simulation/m_start_up.fpp

@@ -83,7 +83,7 @@ contains
 
         namelist /user_inputs/ case_dir, run_time_info, m, n, p, dt, &
             t_step_start, t_step_stop, t_step_save, t_step_print, &
-            model_eqns, mpp_lim, time_stepper, weno_eps, &
+            model_eqns, mpp_lim, time_stepper, weno_eps, muscl_eps, &
             rdma_mpi, teno_CT, mp_weno, weno_avg, &
             riemann_solver, low_Mach, wave_speeds, avg_state, &
             bc_x, bc_y, bc_z, &