Add WENO5 Variants (#481)

Co-authored-by: Cameron <71239694+AiredaleDev@users.noreply.github.com>
MFlowCode · Jun 23, 2024 · 1fb1fd3 · 1fb1fd3
1 parent 41d51b7
commit 1fb1fd3
Show file tree

Hide file tree

Showing 66 changed files with 5,149 additions and 59 deletions.
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -360,14 +360,17 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `time_stepper`         | Integer | Runge--Kutta order [1-3] |
 | `adap_dt`              | Logical | Strang splitting scheme with adaptive time stepping |
 | `weno_order`	         | Integer | WENO order [1,3,5] |
-| `weno_eps`	           | Real    | WENO perturbation (avoid division by zero) |
-| `mapped_weno`	         | Logical | WENO with mapping of nonlinear weights |
+| `weno_eps`	         | Real    | WENO perturbation (avoid division by zero) |
+| `mapped_weno`	         | Logical | WENO-M (WENO with mapping of nonlinear weights) |
+| `wenoz`	             | Logical | WENO-Z |
+| `teno`                 | Logical | TENO (Targeted ENO) |
+| `teno_CT`              | Real    | TENO threshold for smoothness detection | 
 | `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*	 |
-| `avg_state`	           | Integer | Averaged state evaluation method: [1] Roe averagen*; [2] Arithmetic mean  |
+| `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	 |
+| `weno_Re_flux`         | Logical | Compute velocity gradient using scaler divergence theorem	 |
 | `weno_avg`          	 | Logical | Arithmetic mean of left and right, WENO-reconstructed, cell-boundary values |
 
 - \* Options that work only with `model_eqns =2`.
@@ -420,7 +423,13 @@ Note that `time_stepper = 3` specifies the total variation diminishing (TVD), th
 - `weno_eps` specifies the lower bound of the WENO nonlinear weights.
 Practically, `weno_eps` $<10^{-6}$ is used.
 
-- `mapped_weno` activates mapping of the nonlinear WENO weights to the more accurate nonlinear weights in order to reinstate the optimal order of accuracy of the reconstruction in the proximity of critical points ([Henrick et al., 2005](references.md#Henrick05)).
+- `mapped_weno` activates the WENO-M scheme in place of the default WENO-JS scheme ([Henrick et al., 2005](references.md#Henrick05)). WENO-M a variant of the WENO scheme that remaps the nonlinear WENO-JS weights by assigning larger weights to non-smooth stencils, reducing dissipation compared to the default WENO-JS scheme, at the expense of higher computational cost. Only one of `mapped_weno`, `wenoz`, and `teno` can be activated.
+
+- `wenoz` activates the WENO-Z scheme in place of the default WENO-JS scheme ([Borges et al., 2008](references.md#Borges08)). WENO-Z is a variant of the WENO scheme that further reduces the dissipation compared to the WENO-M scheme. It has similar computational cost to the WENO-JS scheme.
+
+- `teno` activates the TENO scheme in place of the default WENO-JS scheme ([Fu et al., 2016](references.md#Fu16)). TENO is a variant of the ENO scheme that is the least dissipative, but could be less robust for extreme cases. It uses a threshold to identify smooth and non-smooth stencils, and applies optimal weights to the smooth stencils. Only available for `weno_order = 5`. Requires `teno_CT` to be set.
+
+- `teno_CT` specifies the threshold for the TENO scheme. This dimensionless constant, also known as $C_T$, sets a threshold to identify smooth and non-smooth stencils. Larger values make the scheme more robust but also more dissipative. A recommended value for teno_CT is `1e-6`. When adjusting this parameter, it is recommended to try values like `1e-5` or `1e-7`. 
 
 - `null_weights` activates nullification of the nonlinear WENO weights at the buffer regions outside the domain boundaries when the Riemann extrapolation boundary condition is specified (`bc_[x,y,z]\%beg[end]} = -4`).
 

diff --git a/docs/documentation/references.md b/docs/documentation/references.md
@@ -20,6 +20,10 @@
 
 - <a id="Henrick05">Henrick, A. K., Aslam, T. D., and Powers, J. M. (2005). Mapped weighted essentially nonoscillatory schemes: achieving optimal order near critical points. Journal of Computational Physics, 207(2):542–567.</a>
 
+- <a id="Borges08">Borges, R., Carmona, M., Costa, B., and Don, W. S. (2008). An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of computational physics, 227(6):3191–3211.</a>
+
+- <a id="Fu16">Fu, L., Hu, X. Y., and Adams, N. A. (2016). A family of high-order targeted ENO schemes for compressible-fluid simulations. Journal of Computational Physics, 305:333–359.</a>
+
 - <a id="Johnsen08">Johnsen, E. (2008). Numerical simulations of non-spherical bubble collapse: With applications to shockwave lithotripsy. PhD thesis, California Institute of Technology.</a>
 
 - <a id="Maeda17">Maeda, K. and Colonius, T. (2017). A source term approach for generation of one-way acoustic waves in the euler and navier–stokes equations. Wave Motion, 75:36–49.</a>

diff --git a/examples/1D_shuosher/README.md → examples/1D_shuosher_old/README.md b/examples/1D_shuosher/README.md → examples/1D_shuosher_old/README.md
diff --git a/examples/1D_shuosher/case.py → examples/1D_shuosher_old/case.py b/examples/1D_shuosher/case.py → examples/1D_shuosher_old/case.py
diff --git a/examples/1D_shuosher/initial.png → examples/1D_shuosher_old/initial.png b/examples/1D_shuosher/initial.png → examples/1D_shuosher_old/initial.png
diff --git a/examples/1D_shuosher/result.png → examples/1D_shuosher_old/result.png b/examples/1D_shuosher/result.png → examples/1D_shuosher_old/result.png
diff --git a/examples/1D_shuosher_teno/case.py b/examples/1D_shuosher_teno/case.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+
+import math
+import json
+
+# Numerical setup
+Nx = 1000
+dx = 1./(1.*(Nx+1))
+
+Tend, Nt = 1.8, 2000
+mydt     = Tend/(1.*Nt)
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : 0.,
+    'x_domain%end'                 : 10.,
+    'm'                            : Nx,
+    'n'                            : 0,
+    'p'                            : 0,
+    'dt'                           : mydt,
+    't_step_start'                 : 0,
+    't_step_stop'                  : int(Nt),
+    't_step_save'                  : int(math.ceil(Nt/10.0)),
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'num_patches'                  : 2,
+    'model_eqns'                   : 2,
+    'alt_soundspeed'               : 'F',
+    'num_fluids'                   : 1,
+    'adv_alphan'                   : 'T',
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'F',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1.E-40,
+    'teno'                         : 'T',
+    'teno_CT'                      : 1.E-6,
+    'null_weights'                 : 'F',
+    'mp_weno'                      : 'F',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     : -3,
+    'bc_x%end'                     : -3,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 2,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'rho_wrt'                      :'T',
+    'parallel_io'                  :'T',
+    # ==========================================================================
+
+
+    # Background to cover whole domain with basic line patch
+    # Patch 1 Left (0 < x < 1) ===============================================
+    'patch_icpp(1)%geometry'       : 1,
+    'patch_icpp(1)%x_centroid'     : 0.5,
+    'patch_icpp(1)%length_x'       : 1.,
+    'patch_icpp(1)%vel(1)'         : 2.629,
+    'patch_icpp(1)%pres'           : 10.333,
+    'patch_icpp(1)%alpha_rho(1)'   : 3.857,
+    'patch_icpp(1)%alpha(1)'       : 1.,
+    # ==========================================================================
+
+
+    # One anlytic patch to take care of 1 < x < 10
+    # Patch 2 Analytic =========================================================
+    'patch_icpp(2)%geometry'        : 1,
+    'patch_icpp(2)%x_centroid'      : 5.5,
+    'patch_icpp(2)%length_x'        : 9.,
+    'patch_icpp(2)%vel(1)'          : 0.,
+    'patch_icpp(2)%pres'            : 1.,
+    'patch_icpp(2)%alpha_rho(1)'    : '1 + 0.2*sin(5*x)',
+    'patch_icpp(2)%alpha(1)'        : 1.,
+    # ==========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.E+00/(1.4-1.E+00),
+    'fluid_pp(1)%pi_inf'           : 0.0,
+    # ==========================================================================
+}))
+
+# ==============================================================================
diff --git a/examples/1D_shuosher_wenojs/case.py b/examples/1D_shuosher_wenojs/case.py
@@ -0,0 +1,90 @@
+#!/usr/bin/env python3
+
+import math
+import json
+
+# Numerical setup
+Nx = 1000
+dx = 1./(1.*(Nx+1))
+
+Tend, Nt = 1.8, 2000
+mydt     = Tend/(1.*Nt)
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : 0.,
+    'x_domain%end'                 : 10.,
+    'm'                            : Nx,
+    'n'                            : 0,
+    'p'                            : 0,
+    'dt'                           : mydt,
+    't_step_start'                 : 0,
+    't_step_stop'                  : int(Nt),
+    't_step_save'                  : int(math.ceil(Nt/10.0)),
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'num_patches'                  : 2,
+    'model_eqns'                   : 2,
+    'alt_soundspeed'               : 'F',
+    'num_fluids'                   : 1,
+    'adv_alphan'                   : 'T',
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'F',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1.E-40,
+    'mapped_weno'                  : 'F',
+    'null_weights'                 : 'F',
+    'mp_weno'                      : 'F',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     : -3,
+    'bc_x%end'                     : -3,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 2,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'rho_wrt'                      :'T',
+    'parallel_io'                  :'T',
+    # ==========================================================================
+
+
+    # Background to cover whole domain with basic line patch
+    # Patch 1 Left (0 < x < 1) ===============================================
+    'patch_icpp(1)%geometry'       : 1,
+    'patch_icpp(1)%x_centroid'     : 0.5,
+    'patch_icpp(1)%length_x'       : 1.,
+    'patch_icpp(1)%vel(1)'         : 2.629,
+    'patch_icpp(1)%pres'           : 10.333,
+    'patch_icpp(1)%alpha_rho(1)'   : 3.857,
+    'patch_icpp(1)%alpha(1)'       : 1.,
+    # ==========================================================================
+
+
+    # One anlytic patch to take care of 1 < x < 10
+    # Patch 2 Analytic =========================================================
+    'patch_icpp(2)%geometry'        : 1,
+    'patch_icpp(2)%x_centroid'      : 5.5,
+    'patch_icpp(2)%length_x'        : 9.,
+    'patch_icpp(2)%vel(1)'          : 0.,
+    'patch_icpp(2)%pres'            : 1.,
+    'patch_icpp(2)%alpha_rho(1)'    : '1 + 0.2*sin(5*x)',
+    'patch_icpp(2)%alpha(1)'        : 1.,
+    # ==========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.E+00/(1.4-1.E+00),
+    'fluid_pp(1)%pi_inf'           : 0.0,
+    # ==========================================================================
+}))
+
+# ==============================================================================
diff --git a/examples/1D_shuosher_wenom/case.py b/examples/1D_shuosher_wenom/case.py
@@ -0,0 +1,90 @@
+#!/usr/bin/env python3
+
+import math
+import json
+
+# Numerical setup
+Nx = 1000
+dx = 1./(1.*(Nx+1))
+
+Tend, Nt = 1.8, 2000
+mydt     = Tend/(1.*Nt)
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : 0.,
+    'x_domain%end'                 : 10.,
+    'm'                            : Nx,
+    'n'                            : 0,
+    'p'                            : 0,
+    'dt'                           : mydt,
+    't_step_start'                 : 0,
+    't_step_stop'                  : int(Nt),
+    't_step_save'                  : int(math.ceil(Nt/10.0)),
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'num_patches'                  : 2,
+    'model_eqns'                   : 2,
+    'alt_soundspeed'               : 'F',
+    'num_fluids'                   : 1,
+    'adv_alphan'                   : 'T',
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'F',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1.E-40,
+    'mapped_weno'                  : 'T',
+    'null_weights'                 : 'F',
+    'mp_weno'                      : 'F',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     : -3,
+    'bc_x%end'                     : -3,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 2,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'rho_wrt'                      :'T',
+    'parallel_io'                  :'T',
+    # ==========================================================================
+
+
+    # Background to cover whole domain with basic line patch
+    # Patch 1 Left (0 < x < 1) ===============================================
+    'patch_icpp(1)%geometry'       : 1,
+    'patch_icpp(1)%x_centroid'     : 0.5,
+    'patch_icpp(1)%length_x'       : 1.,
+    'patch_icpp(1)%vel(1)'         : 2.629,
+    'patch_icpp(1)%pres'           : 10.333,
+    'patch_icpp(1)%alpha_rho(1)'   : 3.857,
+    'patch_icpp(1)%alpha(1)'       : 1.,
+    # ==========================================================================
+
+
+    # One anlytic patch to take care of 1 < x < 10
+    # Patch 2 Analytic =========================================================
+    'patch_icpp(2)%geometry'        : 1,
+    'patch_icpp(2)%x_centroid'      : 5.5,
+    'patch_icpp(2)%length_x'        : 9.,
+    'patch_icpp(2)%vel(1)'          : 0.,
+    'patch_icpp(2)%pres'            : 1.,
+    'patch_icpp(2)%alpha_rho(1)'    : '1 + 0.2*sin(5*x)',
+    'patch_icpp(2)%alpha(1)'        : 1.,
+    # ==========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.E+00/(1.4-1.E+00),
+    'fluid_pp(1)%pi_inf'           : 0.0,
+    # ==========================================================================
+}))
+
+# ==============================================================================