teno5.f90 Source File

!> @file teno5.f90
!> @brief TENO5 (Targeted ENO, 5th-order) reconstruction (Fu, Hu & Adams 2016).
!!
!! Uses the same 5-point stencil and substencil approximations as WENO5-JS
!! (provided by weno_family), but replaces smooth nonlinear weights with a hard
!! threshold that eliminates stencils containing discontinuities before combining
!! the rest using the optimal (linear) weights.  This gives sharper shock
!! resolution than WENO5 while maintaining 5th-order accuracy in smooth regions.
!!
!! Algorithm (Fu, Hu & Adams 2016, Eq. 3.5-3.7):
!!   1. Compute substencil reconstructions vp0..vp2 and smoothness beta0..beta2
!!      (both shared with WENO5-JS via weno_family).
!!   2. Compute global indicator tau5 = |beta0 - beta2| (same as WENO5-Z).
!!   3. Candidate sensitivity: gamma_k = (1 + tau5/(eps+beta_k))^q, q = q_teno.
!!   4. Normalise: chi_k = gamma_k / sum(gamma).
!!   5. Hard cutoff: delta_k = 0 if chi_k < ct_teno, else 1.
!!   6. Weights: alpha_k = delta_k * d_k  (d_k are the WENO5 optimal weights).
!!   7. Fallback (sum(alpha) = 0): use the stencil with the smallest beta_k (ENO).
!!   8. f_hat = sum(alpha_k/sum(alpha) * vp_k).
!!
!! References:
!!   [1] L. Fu, X. Y. Hu, and N. A. Adams, "A family of high-order targeted ENO
!!       schemes for compressible-fluid simulations," J. Comput. Phys.,
!!       305:333-359, 2016.

module teno5
  use precision, only: wp
  use solver_state, only: neq
  use weno_family, only: eps_weno, d0, d1, d2, weno5_substencils, weno5_smoothness
  implicit none
  private

  !> Exponent for the TENO-5 candidate-stencil sensitivity.
  !! q = 6 following Fu, Hu & Adams (2016), Eq. (3.7).
  integer, parameter :: q_teno = 6

  !> Hard-threshold parameter for TENO-5 stencil elimination.
  !! Stencil k is discarded when chi_k < ct_teno.
  !! Fu, Hu & Adams (2016) recommend 1e-5 to 1e-6; smaller values yield
  !! sharper shock resolution but less aggressive oscillation suppression.
  real(wp), parameter :: ct_teno = 1.0e-5_wp

  public :: teno5_reconstruct

contains

  ! ---------------------------------------------------------------------------
  !> Perform TENO5 reconstruction on the supplied stencil.
  !!
  !! @param[in]  f_stencil  Stencil at 5 points (neq x 5)
  !! @param[out] f_hat      Reconstructed value at the interface
  ! ---------------------------------------------------------------------------
  pure subroutine teno5_reconstruct(f_stencil, f_hat)
    real(wp), intent(in) :: f_stencil(:, :)
    real(wp), intent(out) :: f_hat(:)

    real(wp) :: vp0(neq), vp1(neq), vp2(neq)
    real(wp) :: beta0(neq), beta1(neq), beta2(neq)
    real(wp) :: tau5(neq)
    real(wp) :: gamma_0(neq), gamma_1(neq), gamma_2(neq), gamma_sum(neq)  ! candidate sensitivities
    real(wp) :: chi_0(neq), chi_1(neq), chi_2(neq)                        ! normalised sensitivities
    real(wp) :: delta_0(neq), delta_1(neq), delta_2(neq)                  ! binary stencil-active flags
    real(wp) :: alpha_0(neq), alpha_1(neq), alpha_2(neq), alpha_sum(neq)  ! final weights
    integer :: i

    call weno5_substencils(f_stencil, vp0, vp1, vp2)
    call weno5_smoothness(f_stencil, beta0, beta1, beta2)

    ! Global 5th-order smoothness indicator (same as WENO5-Z).
    tau5 = abs(beta0 - beta2)

    ! Candidate stencil sensitivities (Fu et al. 2016, Eq. 3.5).
    do i = 1, neq
      gamma_0(i) = (1.0_wp + tau5(i) / (eps_weno + beta0(i)))**q_teno
      gamma_1(i) = (1.0_wp + tau5(i) / (eps_weno + beta1(i)))**q_teno
      gamma_2(i) = (1.0_wp + tau5(i) / (eps_weno + beta2(i)))**q_teno
    end do
    gamma_sum = gamma_0 + gamma_1 + gamma_2

    ! Normalise sensitivities (Fu et al. 2016, Eq. 3.6).
    chi_0 = gamma_0 / gamma_sum
    chi_1 = gamma_1 / gamma_sum
    chi_2 = gamma_2 / gamma_sum

    ! Hard cutoff: discard stencils below threshold (Fu et al. 2016, Eq. 3.7).
    where (chi_0 < ct_teno)
      delta_0 = 0.0_wp
    elsewhere
      delta_0 = 1.0_wp
    end where
    where (chi_1 < ct_teno)
      delta_1 = 0.0_wp
    elsewhere
      delta_1 = 1.0_wp
    end where
    where (chi_2 < ct_teno)
      delta_2 = 0.0_wp
    elsewhere
      delta_2 = 1.0_wp
    end where

    ! Assemble weights using optimal WENO5 linear weights.
    alpha_0 = delta_0 * d0
    alpha_1 = delta_1 * d1
    alpha_2 = delta_2 * d2
    alpha_sum = alpha_0 + alpha_1 + alpha_2

    ! Combine stencils; fallback to ENO when all stencils eliminated.
    do i = 1, neq
      if (alpha_sum(i) < tiny(1.0_wp)) then
        ! ENO fallback: pick substencil with minimum smoothness indicator.
        if (beta0(i) <= beta1(i) .and. beta0(i) <= beta2(i)) then
          f_hat(i) = vp0(i)
        else if (beta1(i) <= beta2(i)) then
          f_hat(i) = vp1(i)
        else
          f_hat(i) = vp2(i)
        end if
      else
        f_hat(i) = (alpha_0(i) * vp0(i) + alpha_1(i) * vp1(i) &
                &  + alpha_2(i) * vp2(i)) / alpha_sum(i)
      end if
    end do
  end subroutine teno5_reconstruct

end module teno5