weno5_z.f90 Source File

!> @file weno5_z.f90
!> @brief WENO5-Z (Borges et al. 2008) reconstruction scheme.
!!
!! 5th-order WENO reconstruction with the Z-weight formula.  Uses the same
!! substencil reconstructions and smoothness indicators as WENO5-JS but
!! replaces the Jiang-Shu weight formula with the Z-weight formula, which
!! achieves full 5th-order accuracy at smooth extrema where WENO5-JS degrades
!! to 3rd order.
!!
!! Weight formula (Borges et al. 2008, Eq. 3.4-3.5):
!!   tau5 = |beta0 - beta2|                         (global indicator)
!!   alpha_k = d_k * (1 + (tau5/(eps+beta_k))^p),  p = p_wenoz = 2
!!   omega_k  = alpha_k / sum(alpha)
!!
!! References:
!!   [1] R. Borges, M. Carmona, B. Costa, and W. S. Don, "An improved weighted
!!       essentially non-oscillatory scheme for hyperbolic conservation laws,"
!!       J. Comput. Phys., 227:3191-3211, 2008.
!!   [2] G.-S. Jiang and C.-W. Shu, "Efficient Implementation of Weighted ENO
!!       Schemes," J. Comput. Phys., 126:202-228, 1996.

module weno5_z
  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 WENO-Z enhancement factor.
  !! p = 2 gives better shock sharpness than p = 1.
  !! See Borges et al. (2008), Eq. (3.5).
  integer, parameter :: p_wenoz = 2

  public :: weno5_z_reconstruct

contains

  ! ---------------------------------------------------------------------------
  !> Perform WENO5-Z 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 weno5_z_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), w0(neq), w1(neq), w2(neq), w_sum(neq)
    integer :: i

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

    ! Global 5th-order smoothness indicator, Borges et al. (2008) Eq. (3.4).
    tau5 = abs(beta0 - beta2)

    ! WENO-Z nonlinear weights, Borges et al. (2008) Eq. (3.5).
    !   alpha_k = d_k * (1 + (tau5 / (eps + beta_k))^p_wenoz)
    do i = 1, neq
      w0(i) = d0 * (1.0_wp + (tau5(i) / (eps_weno + beta0(i)))**p_wenoz)
      w1(i) = d1 * (1.0_wp + (tau5(i) / (eps_weno + beta1(i)))**p_wenoz)
      w2(i) = d2 * (1.0_wp + (tau5(i) / (eps_weno + beta2(i)))**p_wenoz)
    end do
    w_sum = w0 + w1 + w2

    do i = 1, neq
      f_hat(i) = (vp0(i) * w0(i) + vp1(i) * w1(i) + vp2(i) * w2(i)) / w_sum(i)
    end do
  end subroutine weno5_z_reconstruct

end module weno5_z