!> @file muscl.f90 !> @brief MUSCL reconstruction with selectable TVD limiters (van Leer 1979). !! !! Provides five MUSCL reconstructors, one per TVD limiter. The dispatcher !! in the reconstruction module selects the appropriate procedure pointer at !! initialisation time based on the limiter name supplied by the user. !! !! Stencil layout (half_stencil=2, stencil_width=3): !! col 1 = f_{i-1}, col 2 = f_i (centre), col 3 = f_{i+1} !! !! Formula (component-wise, left-biased stencil): !! ΔL = s(k,2) - s(k,1) !! ΔR = s(k,3) - s(k,2) !! r = ΔR / ΔL (when |ΔL| > tiny, else 0 → phi = 0) !! f_hat(k) = s(k,2) + 0.5 · phi(r) · ΔL !! !! The right-biased (reversed) stencil for F^- is assembled by !! spatial_discretization in mirror order; the same formula gives the !! correct downwind reconstruction automatically. !! !! References: !! [1] B. van Leer, "Towards the Ultimate Conservative Difference Scheme. V. !! A Second-Order Sequel to Godunov's Method," J. Comput. Phys., !! 32:101-136, 1979. module muscl use precision, only: wp implicit none private public :: muscl_minmod_reconstruct, muscl_superbee_reconstruct public :: muscl_mc_reconstruct, muscl_van_leer_reconstruct public :: muscl_koren_reconstruct !> Abstract interface for a pure scalar TVD limiter function phi(r). abstract interface pure function limiter_iface(r) result(phi) import :: wp real(wp), intent(in) :: r real(wp) :: phi end function limiter_iface end interface contains ! =========================================================================== ! TVD limiter helper functions (private, pure) ! =========================================================================== ! --------------------------------------------------------------------------- !> Minmod TVD limiter. Most dissipative of the standard limiters; guaranteed !! TVD for any CFL within the stability limit. !! !! phi(r) = max(0, min(1, r)) !! !! See Roe (1986), Sweby (1984). ! --------------------------------------------------------------------------- pure function lim_minmod(r) result(phi) real(wp), intent(in) :: r !< ratio of successive differences real(wp) :: phi phi = max(0.0_wp, min(1.0_wp, r)) end function lim_minmod ! --------------------------------------------------------------------------- !> Superbee TVD limiter. Least dissipative TVD limiter; may over-sharpen !! contact discontinuities on smooth problems. !! !! phi(r) = max(0, min(1, 2r), min(2, r)) !! !! See Roe (1985). ! --------------------------------------------------------------------------- pure function lim_superbee(r) result(phi) real(wp), intent(in) :: r real(wp) :: phi phi = max(0.0_wp, min(1.0_wp, 2.0_wp * r), min(2.0_wp, r)) end function lim_superbee ! --------------------------------------------------------------------------- !> Monotonized Central (MC) TVD limiter. Good balance between dissipation !! and resolution; recommended for most shock-capturing applications. !! !! phi(r) = max(0, min((1+r)/2, 2, 2r)) !! !! See van Leer (1977). ! --------------------------------------------------------------------------- pure function lim_mc(r) result(phi) real(wp), intent(in) :: r real(wp) :: phi phi = max(0.0_wp, min(0.5_wp * (1.0_wp + r), 2.0_wp, 2.0_wp * r)) end function lim_mc ! --------------------------------------------------------------------------- !> van Leer smooth TVD limiter. !! !! phi(r) = (r + |r|) / (1 + |r|) = 2·max(r,0) / (1 + |r|) !! !! Smooth, symmetric, and achieves 2nd-order accuracy on smooth data. !! See van Leer (1974). ! --------------------------------------------------------------------------- pure function lim_van_leer(r) result(phi) real(wp), intent(in) :: r real(wp) :: phi phi = (r + abs(r)) / (1.0_wp + abs(r)) end function lim_van_leer ! --------------------------------------------------------------------------- !> Koren TVD limiter. Third-order accurate for smooth data; preserves !! monotonicity near discontinuities. !! !! phi(r) = max(0, min(2r, (1+2r)/3, 2)) !! !! See Koren (1993). ! --------------------------------------------------------------------------- pure function lim_koren(r) result(phi) real(wp), intent(in) :: r real(wp) :: phi phi = max(0.0_wp, min(2.0_wp * r, (1.0_wp + 2.0_wp * r) / 3.0_wp, 2.0_wp)) end function lim_koren ! =========================================================================== ! Core MUSCL kernel and public thin wrappers (one per limiter) ! =========================================================================== ! --------------------------------------------------------------------------- !> Core MUSCL reconstruction kernel (component-wise, left-biased stencil). !! !! Each component k of f_hat is computed as: !! ΔL = s(k,2) - s(k,1), ΔR = s(k,3) - s(k,2) !! r = ΔR / ΔL (when |ΔL| > tiny; else f_hat(k) = s(k,2)) !! f_hat(k) = s(k,2) + 0.5 · lim_fn(r) · ΔL !! !! @param[in] f_stencil Stencil (neq x 3) !! @param[out] f_hat Reconstructed value at the interface !! @param[in] lim_fn TVD limiter function phi(r) ! --------------------------------------------------------------------------- pure subroutine muscl_reconstruct(f_stencil, f_hat, lim_fn) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) procedure(limiter_iface) :: lim_fn integer :: k real(wp) :: dL, dR, r do k = 1, size(f_stencil, 1) dL = f_stencil(k, 2) - f_stencil(k, 1) dR = f_stencil(k, 3) - f_stencil(k, 2) if (abs(dL) < tiny(1.0_wp)) then f_hat(k) = f_stencil(k, 2) else r = dR / dL f_hat(k) = f_stencil(k, 2) + 0.5_wp * lim_fn(r) * dL end if end do end subroutine muscl_reconstruct !> MUSCL reconstruction with the Minmod limiter. !! @param[in] f_stencil Stencil (neq x 3) !! @param[out] f_hat Reconstructed value at the interface pure subroutine muscl_minmod_reconstruct(f_stencil, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) call muscl_reconstruct(f_stencil, f_hat, lim_minmod) end subroutine muscl_minmod_reconstruct !> MUSCL reconstruction with the Superbee limiter. !! @param[in] f_stencil Stencil (neq x 3) !! @param[out] f_hat Reconstructed value at the interface pure subroutine muscl_superbee_reconstruct(f_stencil, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) call muscl_reconstruct(f_stencil, f_hat, lim_superbee) end subroutine muscl_superbee_reconstruct !> MUSCL reconstruction with the Monotonized Central (MC) limiter. !! @param[in] f_stencil Stencil (neq x 3) !! @param[out] f_hat Reconstructed value at the interface pure subroutine muscl_mc_reconstruct(f_stencil, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) call muscl_reconstruct(f_stencil, f_hat, lim_mc) end subroutine muscl_mc_reconstruct !> MUSCL reconstruction with the van Leer smooth limiter. !! @param[in] f_stencil Stencil (neq x 3) !! @param[out] f_hat Reconstructed value at the interface pure subroutine muscl_van_leer_reconstruct(f_stencil, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) call muscl_reconstruct(f_stencil, f_hat, lim_van_leer) end subroutine muscl_van_leer_reconstruct !> MUSCL reconstruction with the Koren limiter. !! @param[in] f_stencil Stencil (neq x 3) !! @param[out] f_hat Reconstructed value at the interface pure subroutine muscl_koren_reconstruct(f_stencil, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) call muscl_reconstruct(f_stencil, f_hat, lim_koren) end subroutine muscl_koren_reconstruct end module muscl