!> @file weno_cu6.f90 !> @brief WENO-CU6 (Hu, Wang & Adams 2010) reconstruction scheme. !! !! Uses four 3-point finite-difference candidate polynomials on a 6-point !! face-centred stencil, together with a 6-point reference smoothness !! indicator. Near strong discontinuities the routine falls back to the !! solver's WENO5-JS branch on the 5-point upwind core to preserve robustness. module weno_cu6 use precision, only: wp use solver_state, only: neq use weno_family, only: eps_weno, d0, d1, d2, compute_candidate_reconstructions, & & compute_candidate_smoothness, compute_full_stencil_smoothness, & & weno5_substencils, weno5_smoothness implicit none private real(wp), parameter :: cu6_c = 20.0_wp integer, parameter :: cu6_q = 4 real(wp), parameter :: fallback_ratio = 1.0_wp real(wp), parameter :: candidate_coeffs(3, 4) = reshape([ & 1.0_wp / 3.0_wp, & -7.0_wp / 6.0_wp, & 11.0_wp / 6.0_wp, & -1.0_wp / 6.0_wp, & 5.0_wp / 6.0_wp, & 1.0_wp / 3.0_wp, & 1.0_wp / 3.0_wp, & 5.0_wp / 6.0_wp, & -1.0_wp / 6.0_wp, & 11.0_wp / 6.0_wp, & -7.0_wp / 6.0_wp, & 1.0_wp / 3.0_wp & ], [3, 4]) real(wp), parameter :: linear_weights(4) = [ & 1.0_wp / 20.0_wp, & 9.0_wp / 20.0_wp, & 9.0_wp / 20.0_wp, & 1.0_wp / 20.0_wp] real(wp), parameter :: beta_mats(3, 3, 4) = reshape([ & 4.0_wp / 3.0_wp, & -19.0_wp / 3.0_wp, & 11.0_wp / 3.0_wp, & -19.0_wp / 3.0_wp, & 25.0_wp / 3.0_wp, & -31.0_wp / 3.0_wp, & 11.0_wp / 3.0_wp, & -31.0_wp / 3.0_wp, & 10.0_wp / 3.0_wp, & 4.0_wp / 3.0_wp, & -13.0_wp / 3.0_wp, & 5.0_wp / 3.0_wp, & -13.0_wp / 3.0_wp, & 13.0_wp / 3.0_wp, & -13.0_wp / 3.0_wp, & 5.0_wp / 3.0_wp, & -13.0_wp / 3.0_wp, & 4.0_wp / 3.0_wp, & 10.0_wp / 3.0_wp, & -31.0_wp / 3.0_wp, & 11.0_wp / 3.0_wp, & -31.0_wp / 3.0_wp, & 25.0_wp / 3.0_wp, & -19.0_wp / 3.0_wp, & 11.0_wp / 3.0_wp, & -19.0_wp / 3.0_wp, & 4.0_wp / 3.0_wp, & 22.0_wp / 3.0_wp, & -73.0_wp / 3.0_wp, & 29.0_wp / 3.0_wp, & -73.0_wp / 3.0_wp, & 61.0_wp / 3.0_wp, & -49.0_wp / 3.0_wp, & 29.0_wp / 3.0_wp, & -49.0_wp / 3.0_wp, & 10.0_wp / 3.0_wp & ], [3, 3, 4]) real(wp), parameter :: beta6_mat(6, 6) = reshape([ & 525910327.0_wp / 232243200.0_wp, & -456216463.0_wp / 23224320.0_wp, & 389975071.0_wp / 11612160.0_wp, & -330534727.0_wp / 11612160.0_wp, & 279429607.0_wp / 23224320.0_wp, & -236379487.0_wp / 116121600.0_wp, & -456216463.0_wp / 23224320.0_wp, & 2146987907.0_wp / 46448640.0_wp, & -1930601747.0_wp / 11612160.0_wp, & 1690889819.0_wp / 11612160.0_wp, & -1463230907.0_wp / 23224320.0_wp, & 251883319.0_wp / 23224320.0_wp, & 389975071.0_wp / 11612160.0_wp, & -1930601747.0_wp / 11612160.0_wp, & 1833221603.0_wp / 11612160.0_wp, & -1679332331.0_wp / 5806080.0_wp, & 1495974539.0_wp / 11612160.0_wp, & -263126407.0_wp / 11612160.0_wp, & -330534727.0_wp / 11612160.0_wp, & 1690889819.0_wp / 11612160.0_wp, & -1679332331.0_wp / 5806080.0_wp, & 1607794163.0_wp / 11612160.0_wp, & -1486026707.0_wp / 11612160.0_wp, & 268747951.0_wp / 11612160.0_wp, & 279429607.0_wp / 23224320.0_wp, & -1463230907.0_wp / 23224320.0_wp, & 1495974539.0_wp / 11612160.0_wp, & -1486026707.0_wp / 11612160.0_wp, & 1432381427.0_wp / 46448640.0_wp, & -268475791.0_wp / 23224320.0_wp, & -236379487.0_wp / 116121600.0_wp, & 251883319.0_wp / 23224320.0_wp, & -263126407.0_wp / 11612160.0_wp, & 268747951.0_wp / 11612160.0_wp, & -268475791.0_wp / 23224320.0_wp, & 263126407.0_wp / 232243200.0_wp & ], [6, 6]) public :: weno_cu6_reconstruct contains pure subroutine weno5_fallback_reconstruct(f_stencil, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) real(wp) :: vp0(neq), vp1(neq), vp2(neq) real(wp) :: beta0(neq), beta1(neq), beta2(neq) real(wp) :: w0(neq), w1(neq), w2(neq), w_sum(neq) integer :: i call weno5_substencils(f_stencil(:, 1:5), vp0, vp1, vp2) call weno5_smoothness(f_stencil(:, 1:5), beta0, beta1, beta2) do i = 1, neq w0(i) = d0 / (eps_weno + beta0(i))**2 w1(i) = d1 / (eps_weno + beta1(i))**2 w2(i) = d2 / (eps_weno + beta2(i))**2 end do w_sum = w0 + w1 + w2 do i = 1, neq f_hat(i) = (vp0(i) * w0(i) + vp1(i) * w1(i) + vp2(i) * w2(i)) / w_sum(i) end do end subroutine weno5_fallback_reconstruct pure subroutine weno_cu6_reconstruct(f_stencil, f_hat) real(wp), intent(in) :: f_stencil(:, :) real(wp), intent(out) :: f_hat(:) real(wp) :: q_candidate(size(f_hat), 4) real(wp) :: beta(size(f_hat), 4) real(wp) :: beta_ref(size(f_hat)) real(wp) :: beta_bar, sensor real(wp) :: tau6, alpha(4), w_sum integer :: i call compute_candidate_reconstructions(f_stencil, candidate_coeffs, q_candidate) call compute_candidate_smoothness(f_stencil, beta_mats, beta) call compute_full_stencil_smoothness(f_stencil, beta6_mat, beta_ref) do i = 1, size(f_hat) beta_bar = sum(linear_weights * beta(i, :)) tau6 = abs(beta_ref(i) - beta_bar) sensor = tau6 / (eps_weno + beta_bar) if (sensor > fallback_ratio) then call weno5_fallback_reconstruct(f_stencil, f_hat) return end if alpha = linear_weights * (cu6_c + tau6 / (eps_weno + beta(i, :)))**cu6_q w_sum = sum(alpha) f_hat(i) = sum(alpha * q_candidate(i, :)) / w_sum end do end subroutine weno_cu6_reconstruct end module weno_cu6