mp5.f90 Source File

!> @file mp5.f90
!> @brief MP5 (Monotonicity-Preserving 5th-order) reconstruction (Suresh & Huynh 1997).
!!
!! Applies a 5th-order unlimited reconstruction and clips it to a
!! monotonicity-preserving bound computed from the local curvature.  Achieves
!! 5th-order accuracy in smooth regions while suppressing oscillations near
!! discontinuities.
!!
!! References:
!!   [1] A. Suresh and H. T. Huynh, "Accurate Monotonicity-Preserving Schemes
!!       with Runge-Kutta Time Stepping," J. Comput. Phys., 136:83-99, 1997.

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

  !> MP5 monotonicity parameter α.  Used in the upper-bound estimate.
  !! α = 4 following Suresh & Huynh (1997), Eq. (2.27).
  real(wp), parameter :: alpha_mp5 = 4.0_wp

  !> MP5 curvature-limiter coefficient β_L.
  !! β_L = 4/3 following Suresh & Huynh (1997), Eq. (2.29).
  real(wp), parameter :: beta_mp5 = 4.0_wp / 3.0_wp

  public :: mp5_reconstruct

contains

  ! ---------------------------------------------------------------------------
  !> Two-argument minmod function.
  !! Returns the value with the smaller magnitude if a and b have the same sign,
  !! else 0.
  ! ---------------------------------------------------------------------------
  pure elemental function minmod2(a, b) result(m)
    real(wp), intent(in) :: a, b
    real(wp) :: m
    if (a * b <= 0.0_wp) then
      m = 0.0_wp
    else if (abs(a) <= abs(b)) then
      m = a
    else
      m = b
    end if
  end function minmod2

  ! ---------------------------------------------------------------------------
  !> Four-argument minmod function.
  !! Returns 0 if the arguments do not all share the same sign; otherwise
  !! returns the value with the smallest magnitude.
  ! ---------------------------------------------------------------------------
  pure elemental function minmod4(a, b, c, d) result(m)
    real(wp), intent(in) :: a, b, c, d
    real(wp) :: m
    integer :: s

    ! Derive an integer sign token via sign() to preserve IEEE-754 signed-zero
    ! semantics (sign(1.0,-0.0) == -1) without real-equality comparisons.
    s = int(sign(1.0_wp, a))
    if (s /= int(sign(1.0_wp, b)) .or. s /= int(sign(1.0_wp, c)) .or. &
        s /= int(sign(1.0_wp, d))) then
      m = 0.0_wp
    else
      m = real(s, wp) * min(abs(a), abs(b), abs(c), abs(d))
    end if
  end function minmod4

  ! ---------------------------------------------------------------------------
  !> Perform MP5 reconstruction on the supplied stencil.
  !!
  !! Applies MP5 reconstruction component-wise to a 5-point stencil.
  !! The unlimited 5th-order approximation is identical to the optimal (linear)
  !! WENO5 reconstruction.  A two-stage monotonicity check and curvature-based
  !! limiter enforce non-oscillatory behaviour near discontinuities while
  !! preserving 5th-order accuracy in smooth regions.
  !!
  !! Stencil column mapping (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}.
  !!
  !! See Suresh & Huynh (1997), Sec. 2.3, Eqs. (2.27)-(2.29).
  !!
  !! @param[in]  f_stencil  Stencil at 5 points (neq x 5)
  !! @param[out] f_hat      Reconstructed value at the interface
  ! ---------------------------------------------------------------------------
  pure subroutine mp5_reconstruct(f_stencil, f_hat)
    real(wp), intent(in) :: f_stencil(:, :)
    real(wp), intent(out) :: f_hat(:)

    integer :: k
    real(wp) :: v1, v2, v3, v4, v5
    real(wp) :: u_unlimited, u_mp
    real(wp) :: d0, d1, d2, d_right_lcm, d_left_lcm
    real(wp) :: u_md, u_lc, u_min, u_max

    do k = 1, neq
      v1 = f_stencil(k, 1)
      v2 = f_stencil(k, 2)
      v3 = f_stencil(k, 3)
      v4 = f_stencil(k, 4)
      v5 = f_stencil(k, 5)

      ! --- Step 1: Unlimited 5th-order approximation ---
      ! Eq. (2.26) of Suresh & Huynh (1997).  Identical to optimal WENO5.
      u_unlimited = (2.0_wp * v1 - 13.0_wp * v2 + 47.0_wp * v3 &
                 & + 27.0_wp * v4 - 3.0_wp * v5) / 60.0_wp

      ! --- Step 2: MP upper bound; accept if u_unlimited already monotone ---
      ! u_MP = v3 + minmod2(v4-v3, alpha*(v3-v2)),  alpha=4.
      u_mp = v3 + minmod2(v4 - v3, alpha_mp5 * (v3 - v2))

      if ((u_unlimited - v3) * (u_unlimited - u_mp) <= 0.0_wp) then
        f_hat(k) = u_unlimited
        cycle
      end if

      ! --- Step 3: Curvature-based limiting (Eqs. 2.28-2.29) ---
      d0 = v3 - 2.0_wp * v2 + v1          ! 2nd difference at i-1
      d1 = v4 - 2.0_wp * v3 + v2          ! 2nd difference at i
      d2 = v5 - 2.0_wp * v4 + v3          ! 2nd difference at i+1

      d_right_lcm = minmod4(4.0_wp * d1 - d2, 4.0_wp * d2 - d1, d1, d2)
      d_left_lcm = minmod4(4.0_wp * d0 - d1, 4.0_wp * d1 - d0, d0, d1)

      u_md = 0.5_wp * (v3 + v4) - 0.5_wp * d_right_lcm
      u_lc = v3 + 0.5_wp * (v3 - v2) + beta_mp5 * d_left_lcm

      u_min = max(min(v3, v4, u_md), min(v3, u_unlimited, u_lc))
      u_max = min(max(v3, v4, u_md), max(v3, u_unlimited, u_lc))

      f_hat(k) = u_unlimited + minmod2(u_min - u_unlimited, u_max - u_unlimited)
    end do
  end subroutine mp5_reconstruct

end module mp5