positivity_limiter.f90 Source File

!> @file positivity_limiter.f90
!> @brief Zhang-Shu positivity-preserving limiter for face-reconstructed states.
!!
!! Applied on the FDS path (AUSM+, HLL, HLLC, Roe) after reconstruction.
!! Scales each reconstructed face state toward the adjacent cell average so that
!! density and pressure remain non-negative, preventing NaN propagation on
!! challenging flows.
!!
!! The limiter is a no-op when the reconstructed state is already physical, so
!! it does not degrade accuracy in smooth regions.
!!
!! Usage: set `state%use_positivity_lim = .true.` before calling compute_resid(state).
!!
!! Reference:
!!   X. Zhang & C.-W. Shu, "On maximum-principle-satisfying high order schemes
!!   for scalar conservation laws," J. Comput. Phys. 229:3091-3120, 2010.

module positivity_limiter
  use precision, only: wp
  implicit none
  private

  public :: limit_positivity

contains

  ! ---------------------------------------------------------------------------
  !> Scale a reconstructed face state toward the cell average to ensure
  !! density ρ ≥ ε and pressure p ≥ ε (Zhang-Shu 2010).
  !!
  !! Algorithm (applied sequentially):
  !!   1. If ρ_face < ε:
  !!        θ = (ρ_avg - ε) / (ρ_avg - ρ_face),  clamped to [0, 1]
  !!        q_face ← q_avg + θ * (q_face - q_avg)
  !!   2. Recompute pressure from (updated) q_face.
  !!      If p_face < ε:
  !!        θ = (p_avg - ε) / (p_avg - p_face),  clamped to [0, 1]
  !!        q_face ← q_avg + θ * (q_face - q_avg)
  !!
  !! In smooth regions both steps are no-ops (θ = 1) and the limiter does not
  !! alter the reconstructed state.
  !!
  !! @param[inout] q_face  Reconstructed conserved state at the face (neq=3).
  !!                        Modified in-place.
  !! @param[in]    q_avg   Cell-average conserved state adjacent to this face.
  !! @param[in]    gam     Ratio of specific heats (used for pressure calculation).
  ! ---------------------------------------------------------------------------
  subroutine limit_positivity(q_face, q_avg, gam)
    real(wp), intent(inout) :: q_face(:)
    real(wp), intent(in) :: q_avg(:)
    real(wp), intent(in) :: gam

    !> Small floor for density and pressure.  Below this value a state is
    !! considered unphysical and the limiter activates.
    !! Chosen as ≈100×machine-ε (~2.2e-16) so θ is well-conditioned while
    !! remaining negligible for any physically realistic flow.
    real(wp), parameter :: eps = 1.0e-13_wp

    real(wp) :: rho_avg, rho_face, p_avg, p_face, u_avg, u_face, theta

    rho_avg = q_avg(1)
    rho_face = q_face(1)

    ! ---- Step 1: enforce ρ ≥ ε ----
    if (rho_face < eps) then
      ! rho_avg must be positive (cell average should be physical).
      ! theta blends face toward avg until rho = eps at the face.
      theta = (rho_avg - eps) / (rho_avg - rho_face)
      theta = min(1.0_wp, max(0.0_wp, theta))
      q_face = q_avg + theta * (q_face - q_avg)
    end if

    ! ---- Step 2: enforce p ≥ ε ----
    ! Recompute pressure from (possibly updated) q_face.
    rho_face = q_face(1)
    u_face = q_face(2) / rho_face
    p_face = (q_face(3) - 0.5_wp * rho_face * u_face**2) * (gam - 1.0_wp)

    if (p_face < eps) then
      u_avg = q_avg(2) / rho_avg
      p_avg = (q_avg(3) - 0.5_wp * rho_avg * u_avg**2) * (gam - 1.0_wp)

      theta = (p_avg - eps) / (p_avg - p_face)
      theta = min(1.0_wp, max(0.0_wp, theta))
      q_face = q_avg + theta * (q_face - q_avg)
    end if

  end subroutine limit_positivity

end module positivity_limiter