eno3.f90 Source File


Source Code

!> @file eno3.f90
!> @brief Classical ENO3 (Harten et al. 1987) reconstruction scheme.
!!
!! Selects the smoothest 3-cell substencil per equation component using
!! Newton divided differences.
!!
!! References:
!!   [1] A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy, "Uniformly
!!       high-order accurate essentially non-oscillatory schemes, III,"
!!       J. Comput. Phys., 71:231-303, 1987.

module eno3
  use precision, only: wp
  use solver_state, only: neq
  implicit none
  private

  public :: eno3_reconstruct

contains

  ! ---------------------------------------------------------------------------
  !> Perform classical ENO3 reconstruction on the supplied stencil.
  !!
  !! Selects the smoothest 3-cell substencil per equation component using
  !! Newton divided differences (Harten et al., 1987).
  !!
  !! Stencil column mapping (5-point, half_stencil=3):
  !!   col 1 = v_{i-2},  col 2 = v_{i-1},  col 3 = v_i (centre),
  !!   col 4 = v_{i+1},  col 5 = v_{i+2}.
  !!
  !! Two-step stencil selection per component:
  !!   Step 1 — compare |d^1[i-1]| vs |d^1[i]|  (1st-order divided differences):
  !!            d^1[j] = v_{j+1} - v_j
  !!   Step 2 — compare |d^2[s-1]| vs |d^2[s]|  (2nd-order divided differences):
  !!            d^2[j] = (v_{j+2} - 2*v_{j+1} + v_j) / 2
  !!
  !! Selected substencil starts at column s ∈ {1, 2, 3} and uses the same
  !! reconstruction formulas as the three WENO5 substencils:
  !!   S0 (s=1): f̂ =  (1/3)*v1 - (7/6)*v2 + (11/6)*v3
  !!   S1 (s=2): f̂ = -(1/6)*v1 + (5/6)*v2  + (1/3)*v3
  !!   S2 (s=3): f̂ =  (1/3)*v1 + (5/6)*v2  - (1/6)*v3
  !!
  !! @param[in]  f_stencil  Stencil at 5 points (neq x 5)
  !! @param[out] f_hat      Reconstructed value at the interface
  ! ---------------------------------------------------------------------------
  pure subroutine eno3_reconstruct(f_stencil, f_hat)
    real(wp), intent(in) :: f_stencil(:, :)
    real(wp), intent(out) :: f_hat(:)

    real(wp) :: d1_l, d1_r, d2_l, d2_r
    real(wp) :: v1, v2, v3
    integer :: ieq, s

    do ieq = 1, neq
      ! --- Step 1: 1st-order divided differences at i-1 and i ---
      !   d^1[i-1] = f_i     - f_{i-1} = col3 - col2
      !   d^1[i]   = f_{i+1} - f_i     = col4 - col3
      d1_l = f_stencil(ieq, 3) - f_stencil(ieq, 2)
      d1_r = f_stencil(ieq, 4) - f_stencil(ieq, 3)

      if (abs(d1_l) <= abs(d1_r)) then
        s = 2   ! 2-cell {i-1, i} = {col2, col3}
      else
        s = 3   ! 2-cell {i, i+1} = {col3, col4}
      end if

      ! --- Step 2: 2nd-order divided differences for left/right extension ---
      !   d^2[s-1] = (col(s+1) - 2*col(s) + col(s-1)) / 2  [left extension]
      !   d^2[s]   = (col(s+2) - 2*col(s+1) + col(s))  / 2  [right extension]
      d2_l = (f_stencil(ieq, s + 1) - 2.0_wp * f_stencil(ieq, s) &
          &   + f_stencil(ieq, s - 1)) * 0.5_wp
      d2_r = (f_stencil(ieq, s + 2) - 2.0_wp * f_stencil(ieq, s + 1) &
          &   + f_stencil(ieq, s)) * 0.5_wp

      if (abs(d2_l) <= abs(d2_r)) s = s - 1   ! extend left

      ! --- Apply selected substencil reconstruction formula ---
      v1 = f_stencil(ieq, s)
      v2 = f_stencil(ieq, s + 1)
      v3 = f_stencil(ieq, s + 2)

      select case (s)
      case (1)   ! S0 = {i-2, i-1, i}
        f_hat(ieq) = 1.0_wp / 3.0_wp * v1 - 7.0_wp / 6.0_wp * v2 + 11.0_wp / 6.0_wp * v3
      case (2)   ! S1 = {i-1, i, i+1}
        f_hat(ieq) = -1.0_wp / 6.0_wp * v1 + 5.0_wp / 6.0_wp * v2 + 1.0_wp / 3.0_wp * v3
      case (3)   ! S2 = {i, i+1, i+2}
        f_hat(ieq) = 1.0_wp / 3.0_wp * v1 + 5.0_wp / 6.0_wp * v2 - 1.0_wp / 6.0_wp * v3
      end select
    end do
  end subroutine eno3_reconstruct

end module eno3