!> @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 weno_family, only: eps_weno, compute_candidate_reconstructions, & & compute_candidate_smoothness, compute_full_stencil_smoothness use weno5_js, only: weno5_js_reconstruct 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 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) ! Deliberate face-level (all-or-nothing) fallback: the loop scans the ! components, but the FIRST component whose discontinuity sensor exceeds ! fallback_ratio reconstructs the ENTIRE face with WENO5 and returns — it ! does NOT fall back only that one component. With fallback_ratio = 1.0, ! a face switches to WENO5 once tau6 exceeds the mean candidate smoothness ! (sensor > 1), keeping the whole face's reconstruction consistent across ! its components rather than mixing CU6 and WENO5 per component. 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 ! Face-level fallback to WENO5-JS on the 5-point upwind core. call weno5_js_reconstruct(f_stencil(:, 1:5), 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