!> @file weno_family.f90 !> @brief Shared infrastructure for WENO-family reconstruction schemes. !! !! Provides the parameters and helper subroutines that are common across the !! solver's WENO-family schemes. The original WENO5 helpers are retained for !! WENO5-JS, WENO5-Z, and TENO5. More general candidate-polynomial and !! smoothness-evaluation helpers support WENO-CU6, WENO9-JS, and WENO11-JS !! without duplicating large blocks of algebra. !! !! Exported parameters: !! eps_weno — small stabiliser to avoid division by zero in smoothness weights. !! d0, d1, d2 — linear (optimal) WENO5 weights; sum = 1. !! !! Exported subroutines: !! weno5_substencils — computes 3rd-order substencil reconstructions vp0..vp2. !! weno5_smoothness — computes Jiang-Shu smoothness indicators beta0..beta2. !! compute_candidate_reconstructions — evaluates candidate face polynomials !! from tabulated coefficients. !! compute_candidate_smoothness — evaluates Jiang-Shu smoothness forms !! from tabulated quadratic matrices. !! compute_full_stencil_smoothness — evaluates a quadratic smoothness form !! on the full stencil (used by CU6). !! apply_weno_js_generic — generic Jiang-Shu nonlinear weighting !! for odd-width WENO schemes. !! !! References: !! [1] G.-S. Jiang and C.-W. Shu, "Efficient Implementation of Weighted ENO !! Schemes," J. Comput. Phys., 126:202-228, 1996. module weno_family use precision, only: wp implicit none private !> Small parameter to avoid division by zero in WENO weights. !! See Jiang & Shu (1996), Eq. (22). !! Typical literature range: 1e-6 to 1e-4; order-of-magnitude changes !! have negligible effect in smooth regions but may affect shock sharpness. real(wp), parameter, public :: eps_weno = 1.0e-6_wp !> Linear (optimal) weights for WENO5 5th-order accuracy. !! gamma_0 = 1/10, gamma_1 = 3/5, gamma_2 = 3/10. !! See Jiang & Shu (1996), Eq. (19). real(wp), parameter, public :: d0 = 0.1_wp real(wp), parameter, public :: d1 = 0.6_wp real(wp), parameter, public :: d2 = 0.3_wp !> Linear (optimal) weights for WENO7 7th-order accuracy. !! d7_0=1/35, d7_1=12/35, d7_2=18/35, d7_3=4/35; sum = 1. !! See Balsara & Shu (2000), Eq. (2.16). real(wp), parameter, public :: d7_0 = 1.0_wp / 35.0_wp real(wp), parameter, public :: d7_1 = 12.0_wp / 35.0_wp real(wp), parameter, public :: d7_2 = 18.0_wp / 35.0_wp real(wp), parameter, public :: d7_3 = 4.0_wp / 35.0_wp public :: weno5_substencils, weno5_smoothness public :: compute_candidate_reconstructions, compute_candidate_smoothness public :: compute_full_stencil_smoothness, apply_weno_js_generic contains ! --------------------------------------------------------------------------- !> Compute the three 3rd-order substencil reconstructions for WENO5. !! !! Operates component-wise on a 5-point stencil. The reconstructed values !! are the same for WENO5-JS, WENO5-Z, and TENO5. !! !! Substencil approximations (Jiang & Shu 1996, Eq. 16a-16c): !! vp0 = (1/3)*v1 - (7/6)*v2 + (11/6)*v3 !! vp1 = -(1/6)*v2 + (5/6)*v3 + (1/3)*v4 !! vp2 = (1/3)*v3 + (5/6)*v4 - (1/6)*v5 !! !! @param[in] f_stencil Stencil values at 5 points (neq x 5) !! @param[out] vp0 Substencil 0 reconstruction (neq) !! @param[out] vp1 Substencil 1 reconstruction (neq) !! @param[out] vp2 Substencil 2 reconstruction (neq) ! --------------------------------------------------------------------------- pure subroutine weno5_substencils(f_stencil, vp0, vp1, vp2) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: vp0(:), vp1(:), vp2(:) vp0 = 1.0_wp / 3.0_wp * f_stencil(:, 1) & & - 7.0_wp / 6.0_wp * f_stencil(:, 2) & & + 11.0_wp / 6.0_wp * f_stencil(:, 3) vp1 = -1.0_wp / 6.0_wp * f_stencil(:, 2) & & + 5.0_wp / 6.0_wp * f_stencil(:, 3) & & + 1.0_wp / 3.0_wp * f_stencil(:, 4) vp2 = 1.0_wp / 3.0_wp * f_stencil(:, 3) & & + 5.0_wp / 6.0_wp * f_stencil(:, 4) & & - 1.0_wp / 6.0_wp * f_stencil(:, 5) end subroutine weno5_substencils ! --------------------------------------------------------------------------- !> Compute the Jiang-Shu smoothness indicators for a 5-point WENO5 stencil. !! !! See Jiang & Shu (1996), Eq. (24)-(26). !! !! @param[in] f_stencil Stencil values at 5 points (neq x 5) !! @param[out] beta0 Smoothness indicator for substencil 0 (neq) !! @param[out] beta1 Smoothness indicator for substencil 1 (neq) !! @param[out] beta2 Smoothness indicator for substencil 2 (neq) ! --------------------------------------------------------------------------- pure subroutine weno5_smoothness(f_stencil, beta0, beta1, beta2) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: beta0(:), beta1(:), beta2(:) integer :: i do i = 1, size(f_stencil, 1) beta0(i) = 13.0_wp / 12.0_wp & & * (f_stencil(i, 1) - 2.0_wp * f_stencil(i, 2) + f_stencil(i, 3))**2 & & + 1.0_wp / 4.0_wp & & * (f_stencil(i, 1) - 4.0_wp * f_stencil(i, 2) + 3.0_wp * f_stencil(i, 3))**2 beta1(i) = 13.0_wp / 12.0_wp & & * (f_stencil(i, 2) - 2.0_wp * f_stencil(i, 3) + f_stencil(i, 4))**2 & & + 1.0_wp / 4.0_wp & & * (f_stencil(i, 2) - f_stencil(i, 4))**2 beta2(i) = 13.0_wp / 12.0_wp & & * (f_stencil(i, 3) - 2.0_wp * f_stencil(i, 4) + f_stencil(i, 5))**2 & & + 1.0_wp / 4.0_wp & & * (3.0_wp * f_stencil(i, 3) - 4.0_wp * f_stencil(i, 4) + f_stencil(i, 5))**2 end do end subroutine weno5_smoothness ! --------------------------------------------------------------------------- !> Evaluate candidate reconstruction polynomials from tabulated coefficients. !! !! The candidate coefficients are stored as coeffs(j, k), where j enumerates !! the local substencil column and k enumerates the candidate stencil. The !! candidate stencils are assumed to be contiguous unit shifts within the !! supplied full stencil. !! !! @param[in] f_stencil Full stencil values (neq x stencil_width) !! @param[in] coeffs Candidate polynomial coefficients !! (candidate_width x n_candidate) !! @param[out] q_candidate Candidate face values (neq x n_candidate) ! --------------------------------------------------------------------------- pure subroutine compute_candidate_reconstructions(f_stencil, coeffs, q_candidate) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(in) :: coeffs(:, :) real(wp), intent(out) :: q_candidate(:, :) integer :: j, k q_candidate = 0.0_wp do k = 1, size(coeffs, 2) do j = 1, size(coeffs, 1) q_candidate(:, k) = q_candidate(:, k) + coeffs(j, k) * f_stencil(:, k + j - 1) end do end do end subroutine compute_candidate_reconstructions ! --------------------------------------------------------------------------- !> Evaluate candidate Jiang-Shu smoothness indicators from quadratic forms. !! !! Each candidate smoothness indicator is represented by a symmetric matrix !! beta_mats(:, :, k) acting on the corresponding contiguous substencil. !! !! @param[in] f_stencil Full stencil values (neq x stencil_width) !! @param[in] beta_mats Quadratic-form matrices !! (candidate_width x candidate_width x n_candidate) !! @param[out] beta Smoothness indicators (neq x n_candidate) ! --------------------------------------------------------------------------- pure subroutine compute_candidate_smoothness(f_stencil, beta_mats, beta) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(in) :: beta_mats(:, :, :) real(wp), intent(out) :: beta(:, :) real(wp) :: coeff integer :: a, b, k beta = 0.0_wp do k = 1, size(beta_mats, 3) do a = 1, size(beta_mats, 1) do b = 1, size(beta_mats, 2) coeff = beta_mats(a, b, k) ! The double loop visits both (a,b) and (b,a); halving the ! off-diagonal avoids double-counting. Stored values equal ! the direct polynomial coefficients (not 2x the off-diagonal). if (a /= b) coeff = 0.5_wp * coeff beta(:, k) = beta(:, k) + coeff * f_stencil(:, k + a - 1) * f_stencil(:, k + b - 1) end do end do end do end subroutine compute_candidate_smoothness ! --------------------------------------------------------------------------- !> Evaluate a quadratic smoothness form on the full stencil. !! !! Used by WENO-CU6 for the sixth-order reference smoothness indicator β6. !! !! @param[in] f_stencil Full stencil values (neq x stencil_width) !! @param[in] beta_mat Quadratic-form matrix (stencil_width x stencil_width) !! @param[out] beta Smoothness indicator (neq) ! --------------------------------------------------------------------------- pure subroutine compute_full_stencil_smoothness(f_stencil, beta_mat, beta) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(in) :: beta_mat(:, :) real(wp), intent(out) :: beta(:) real(wp) :: coeff integer :: a, b beta = 0.0_wp do a = 1, size(beta_mat, 1) do b = 1, size(beta_mat, 2) coeff = beta_mat(a, b) ! Off-diagonal halved: double loop visits (a,b) and (b,a) separately. if (a /= b) coeff = 0.5_wp * coeff beta = beta + coeff * f_stencil(:, a) * f_stencil(:, b) end do end do end subroutine compute_full_stencil_smoothness ! --------------------------------------------------------------------------- !> Generic Jiang-Shu WENO reconstruction for odd-width schemes. !! !! Applies the standard Jiang-Shu nonlinear weights !! alpha_k = d_k / (eps + beta_k)^2 !! to the candidate reconstructions supplied by @p coeffs and @p beta_mats. !! !! @param[in] f_stencil Full stencil values (neq x stencil_width) !! @param[in] coeffs Candidate polynomial coefficients !! (candidate_width x n_candidate) !! @param[in] beta_mats Candidate smoothness matrices !! (candidate_width x candidate_width x n_candidate) !! @param[in] linear_weights Optimal linear weights (n_candidate) !! @param[out] f_hat Reconstructed face value (neq) ! --------------------------------------------------------------------------- pure subroutine apply_weno_js_generic(f_stencil, coeffs, beta_mats, linear_weights, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(in) :: coeffs(:, :) real(wp), intent(in) :: beta_mats(:, :, :) real(wp), intent(in) :: linear_weights(:) real(wp), intent(out) :: f_hat(:) real(wp) :: q_candidate(size(f_hat), size(linear_weights)) real(wp) :: beta(size(f_hat), size(linear_weights)) real(wp) :: alpha(size(linear_weights)), w_sum integer :: i, k call compute_candidate_reconstructions(f_stencil, coeffs, q_candidate) call compute_candidate_smoothness(f_stencil, beta_mats, beta) do i = 1, size(f_hat) do k = 1, size(linear_weights) alpha(k) = linear_weights(k) / (eps_weno + beta(i, k))**2 end do w_sum = sum(alpha) f_hat(i) = sum(alpha * q_candidate(i, :)) / w_sum end do end subroutine apply_weno_js_generic end module weno_family