WENO7

README.md

@@ -132,7 +132,7 @@ They are organized below. Just click the drop-downs!
 
 * Shock and interface capturing schemes
 	* First-order upwinding
- 	* WENO reconstructions of order 3 and 5
+ 	* WENO reconstructions of order 3, 5, and 7
   	* WENO variants: WENO-JS, WENO-M, WENO-Z, TENO
    	* Monotonicity-preserving reconstructions
 	* Reliable handling of high density ratios

docs/documentation/case.md

@@ -123,7 +123,7 @@ For example, $m=n=p=499$ discretizes the domain into $500^3$ cells.
 When the simulation is 2D/axi-symmetric or 1D, it requires that $p=0$ or $p=n=0$, respectively.
 
 - `stretch_[x,y,z]` activates grid stretching in the $[x,y,z]$ directions.
-The grid is gradually stretched such that the domain boundaries are pushed away from the origin along a specified axis.
+The grid is gradually stretched such that the domain boundaries are pushed away from the origin along a specified axis. WENO7 does not support grid stretching.
 
 - `a_[x,y,z]`, `[x,y,z]_a`, and `[x,y,z]_b` are parameters that define the grid stretching function. When grid stretching along the $x$ axis is considered, the stretching function is given as:
 
@@ -347,6 +347,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `weno_eps`	         | Real    | WENO perturbation (avoid division by zero) |
 | `mapped_weno`	         | Logical | WENO-M (WENO with mapping of nonlinear weights) |
 | `wenoz`	             | Logical | WENO-Z |
+| `wenoz_q`              | Real    | WENO-Z power parameter q (only for WENO7) |
 | `teno`                 | Logical | TENO (Targeted ENO) |
 | `teno_CT`              | Real    | TENO threshold for smoothness detection |
 | `null_weights`         | Logical | Null WENO weights at boundaries |
@@ -406,18 +407,20 @@ Note that `time_stepper = 3` specifies the total variation diminishing (TVD), th
 
 - `adap_dt` activates the Strang operator splitting scheme which splits flux and source terms in time marching, and an adaptive time stepping strategy is implemented for the source term. It requires ``bubbles = 'T'``, ``polytropic = 'T'``, ``adv_n = 'T'`` and `time_stepper = 3`.
 
-- `weno_order` specifies the order of WENO scheme that is used for spatial reconstruction of variables by an integer of 1, 3, and 5, that correspond to the 1st, 3rd, and 5th order, respectively.
+- `weno_order` specifies the order of WENO scheme that is used for spatial reconstruction of variables by an integer of 1, 3, 5, and 7, that correspond to the 1st, 3rd, 5th, and 7th order, respectively. WENO7 does not support grid stretching.
 
 - `weno_eps` specifies the lower bound of the WENO nonlinear weights.
-Practically, `weno_eps` $<10^{-6}$ is used.
+It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$ for other WENO variants.
 
 - `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.
+- `wenoz_q` specifies the power parameter `q` used in the WENO-Z scheme. It controls how aggressively the smoothness coefficients scale the weights. A higher value of `wenoz_q` increases the sensitivity to smoothness, improving stability but worsening numerical dissipation. For WENO3 and WENO5, `q=1` is fixed, so `wenoz_q` must not be set. For WENO7, `wenoz_q` can be set to 2, 3, or 4.
 
-- `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`.
+- `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` and `7`. 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` for TENO5. A smaller value can be used for TENO7.
 
 - `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`).
 

examples/1D_shuosher_teno7/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,
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'F',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 7,
+    'weno_eps'                     : 1.E-40,
+    'teno'                         : 'T',
+    'teno_CT'                      : 1.E-9,
+    '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,
+    # ==========================================================================
+}))
+
+# ==============================================================================

examples/1D_shuosher_wenojs/case.py → examples/1D_shuosher_wenojs5/case.py

@@ -37,7 +37,7 @@
     'mixture_err'                  : 'F',
     'time_stepper'                 : 3,
     'weno_order'                   : 5,
-    'weno_eps'                     : 1.E-40,
+    'weno_eps'                     : 1.E-6,
     'mapped_weno'                  : 'F',
     'null_weights'                 : 'F',
     'mp_weno'                      : 'F',

src/common/m_checker_common.fpp

@@ -194,7 +194,7 @@ contains
     !> Checks constraints regarding WENO order.
         !! Called by s_check_inputs_common for all three stages
     subroutine s_check_inputs_weno
-        @:PROHIBIT(all(weno_order /= (/1, 3, 5/)), "weno_order must be 1, 3, or 5")
+        @:PROHIBIT(all(weno_order /= (/1, 3, 5, 7/)), "weno_order must be 1, 3, 5, or 7")
         @:PROHIBIT(m + 1 < weno_order, "m must be at least weno_order - 1")
         @:PROHIBIT(n > 0 .and. n + 1 < weno_order, "n must be at least weno_order - 1")
         @:PROHIBIT(p > 0 .and. p + 1 < weno_order, "p must be at least weno_order - 1")

src/pre_process/m_checker.fpp

@@ -116,6 +116,7 @@ contains
 
         ! Common checks for all directions (stretch_x, stretch_y, and stretch_z)
         #:for X in ['x', 'y', 'z']
+            @:PROHIBIT(stretch_${X}$ .and. weno_order == 7, "weno_order = 7 does not support stretched grids")
             @:PROHIBIT(stretch_${X}$ .and. old_grid, "old_grid and stretch_${X}$ are incompatible")
             @:PROHIBIT(stretch_${X}$ .and. f_is_default(a_${X}$), "a_${X}$ must be set with stretch_${X}$ enabled")
             @:PROHIBIT(stretch_${X}$ .and. f_is_default(${X}$_a), "${X}$_a must be set with stretch_${X}$ enabled")

src/simulation/m_cbc.fpp

@@ -713,7 +713,7 @@ contains
                     ! ==================================================================
 
                     ! PI4 of flux_rs_vf and flux_src_rs_vf at j = 1/2, 3/2 =============
-                elseif (weno_order == 5) then
+                else
                     call s_convert_primitive_to_flux_variables(q_prim_rs${XYZ}$_vf, &
                                                                F_rs${XYZ}$_vf, &
                                                                F_src_rs${XYZ}$_vf, &

src/simulation/m_checker.fpp

@@ -66,12 +66,16 @@ contains
         @:PROHIBIT(weno_order /= 1 .and. f_is_default(weno_eps), &
             "weno_order != 1, but weno_eps is not set. A typical value of weno_eps is 1e-6")
         @:PROHIBIT(weno_eps <= 0d0, "weno_eps must be positive. A typical value of weno_eps is 1e-6")
+        @:PROHIBIT(wenoz .and. weno_order == 7 .and. f_is_default(wenoz_q), &
+            "wenoz is used at 7th order, but wenoz_q is not set. It should be either 2, 3, or 4")
+        @:PROHIBIT(wenoz .and. weno_order == 7 .and. .not. (f_approx_equal(wenoz_q, 2d0) .or. f_approx_equal(wenoz_q, 3d0) .or. f_approx_equal(wenoz_q, 4d0)), &
+            "wenoz_q must be either 2, 3, or 4")
         @:PROHIBIT(teno .and. f_is_default(teno_CT), "teno is used, but teno_CT is not set. A typical value of teno_CT is 1e-6")
         @:PROHIBIT(teno .and. teno_CT <= 0d0, "teno_CT must be positive. A typical value of teno_CT is 1e-6")
         @:PROHIBIT(count([mapped_weno, wenoz, teno]) >= 2, "Only one of mapped_weno, wenoz, or teno can be set to true")
         @:PROHIBIT(weno_order == 1 .and. mapped_weno)
         @:PROHIBIT(weno_order == 1 .and. wenoz)
-        @:PROHIBIT(weno_order /= 5 .and. teno)
+        @:PROHIBIT((weno_order == 1 .or. weno_order == 3) .and. teno)
         @:PROHIBIT(weno_order /= 5 .and. mp_weno)
         @:PROHIBIT(model_eqns == 1 .and. weno_avg)
     end subroutine s_check_inputs_weno

src/simulation/m_global_parameters.fpp

@@ -129,19 +129,23 @@ module m_global_parameters
     #:if MFC_CASE_OPTIMIZATION
         integer, parameter :: weno_polyn = ${weno_polyn}$ !< Degree of the WENO polynomials (polyn)
         integer, parameter :: weno_order = ${weno_order}$ !< Order of the WENO reconstruction
+        integer, parameter :: weno_num_stencils = ${weno_num_stencils}$ !< Number of stencils for WENO reconstruction (only different from weno_polyn for TENO(>5))
         integer, parameter :: num_fluids = ${num_fluids}$ !< number of fluids in the simulation
         logical, parameter :: wenojs = (${wenojs}$ /= 0)            !< WENO-JS (default)
         logical, parameter :: mapped_weno = (${mapped_weno}$ /= 0)  !< WENO-M (WENO with mapping of nonlinear weights)
         logical, parameter :: wenoz = (${wenoz}$ /= 0)              !< WENO-Z
         logical, parameter :: teno = (${teno}$ /= 0)                !< TENO (Targeted ENO)
+        real(kind(0d0)), parameter :: wenoz_q = ${wenoz_q}$         !< Power constant for WENO-Z
     #:else
         integer :: weno_polyn     !< Degree of the WENO polynomials (polyn)
         integer :: weno_order     !< Order of the WENO reconstruction
+        integer :: weno_num_stencils    !< Number of stencils for WENO reconstruction (only different from weno_polyn for TENO(>5))
         integer :: num_fluids     !< number of fluids in the simulation
         logical :: wenojs         !< WENO-JS (default)
         logical :: mapped_weno    !< WENO-M (WENO with mapping of nonlinear weights)
         logical :: wenoz          !< WENO-Z
         logical :: teno           !< TENO (Targeted ENO)
+        real(kind(0d0)) :: wenoz_q  !< Power constant for WENO-Z
     #:endif
 
     real(kind(0d0)) :: weno_eps       !< Binding for the WENO nonlinear weights
@@ -174,7 +178,7 @@ module m_global_parameters
     integer :: cpu_start, cpu_end, cpu_rate
 
     #:if not MFC_CASE_OPTIMIZATION
-        !$acc declare create(num_dims, weno_polyn, weno_order, num_fluids, wenojs, mapped_weno, wenoz, teno)
+        !$acc declare create(num_dims, weno_polyn, weno_order, weno_num_stencils, num_fluids, wenojs, mapped_weno, wenoz, teno, wenoz_q)
     #:endif
 
     !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach)
@@ -551,6 +555,7 @@ contains
             mapped_weno = .false.
             wenoz = .false.
             teno = .false.
+            wenoz_q = dflt_real
         #:endif
 
         chem_params%advection = .false.
@@ -708,7 +713,13 @@ contains
         #:if not MFC_CASE_OPTIMIZATION
             ! Determining the degree of the WENO polynomials
             weno_polyn = (weno_order - 1)/2
+            if (teno) then
+                weno_num_stencils = weno_order - 3
+            else
+                weno_num_stencils = weno_polyn
+            end if
             !$acc update device(weno_polyn)
+            !$acc update device(weno_num_stencils)
             !$acc update device(nb)
             !$acc update device(num_dims, num_fluids)
         #:endif
@@ -1103,6 +1114,7 @@ contains
 
         #:if not MFC_CASE_OPTIMIZATION
             !$acc update device(wenojs, mapped_weno, wenoz, teno)
+            !$acc update device(wenoz_q)
         #:endif
 
         !$acc enter data copyin(nb, R0ref, Ca, Web, Re_inv, weight, R0, V0, bubbles, polytropic, polydisperse, qbmm, R0_type, ptil, bubble_model, thermal, poly_sigma)

src/simulation/m_mpi_proxy.fpp

@@ -230,6 +230,7 @@ contains
             call MPI_BCAST(weno_order, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(nb, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(num_fluids, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(wenoz_q, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
         #:endif
 
         do i = 1, num_fluids_max

src/simulation/m_start_up.fpp

@@ -149,7 +149,7 @@ contains
             rhoref, pref, bubbles, bubble_model, &
             R0ref, chem_params, &
 #:if not MFC_CASE_OPTIMIZATION
-            nb, mapped_weno, wenoz, teno, weno_order, num_fluids, &
+            nb, mapped_weno, wenoz, teno, wenoz_q, weno_order, num_fluids, &
 #:endif
             Ca, Web, Re_inv, &
             acoustic_source, acoustic, num_source, &