Skip to content

feat: WENO5 sum-of-squares beta on uniform grids (cancellation-free smoothness indicators) #1408

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
sbryngelson wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

feat: WENO5 sum-of-squares beta on uniform grids (cancellation-free smoothness indicators) #1408

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
daf8d8f
feat: WENO5 sum-of-squares beta on uniform grids to eliminate cancell…
sbryngelson May 7, 2026
5b0a606
fix: widen uniform-grid tolerance to sqrt(epsilon(h0)) for single-pre…
sbryngelson May 7, 2026
63e7833
Merge branch 'master' into weno-beta-sos
sbryngelson May 7, 2026
5bb6a34
Merge branch 'master' into weno-beta-sos
sbryngelson May 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 35 additions & 9 deletions src/simulation/m_weno.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,9 @@ module m_weno
integer :: v_size !< Number of WENO-reconstructed cell-average variables
$:GPU_DECLARE(create='[v_size]')

logical :: uniform_grid(3) !< True if grid spacing is uniform in each direction
$:GPU_DECLARE(create='[uniform_grid]')

!> @name Indical bounds in the s1-, s2- and s3-directions
!> @{
type(int_bounds_info) :: is1_weno, is2_weno, is3_weno
Expand Down Expand Up @@ -183,6 +186,7 @@ contains
integer :: i !< Generic loop iterator
real(wp) :: w(1:8) !< Intermediate var for ideal weights: s_cb across overall stencil
real(wp) :: y(1:4) !< Intermediate var for poly & beta: diff(s_cb) across sub-stencil
real(wp) :: h0 !< Reference spacing for uniform-grid detection

! Determine cell count, boundary locations, and BCs for selected WENO direction

Expand Down Expand Up @@ -843,12 +847,24 @@ contains
end if
#:endfor

! Detect whether grid spacing is uniform (enables cancellation-free sum-of-squares beta). Tolerance uses sqrt(epsilon) so it
! works in both double and single precision: ~1.5e-8 relative in double, ~3.5e-4 in single - above FP noise, below real
! stretching.
uniform_grid(weno_dir) = .true.
h0 = (s_cb(s) - s_cb(0))/real(s, wp)
do i = 0, s - 1
if (abs((s_cb(i + 1) - s_cb(i)) - h0) > sqrt(epsilon(h0))*abs(h0)) then
uniform_grid(weno_dir) = .false.
exit
end if
end do

if (weno_dir == 1) then
$:GPU_UPDATE(device='[poly_coef_cbL_x, poly_coef_cbR_x, d_cbL_x, d_cbR_x, beta_coef_x]')
$:GPU_UPDATE(device='[poly_coef_cbL_x, poly_coef_cbR_x, d_cbL_x, d_cbR_x, beta_coef_x, uniform_grid]')
else if (weno_dir == 2) then
$:GPU_UPDATE(device='[poly_coef_cbL_y, poly_coef_cbR_y, d_cbL_y, d_cbR_y, beta_coef_y]')
$:GPU_UPDATE(device='[poly_coef_cbL_y, poly_coef_cbR_y, d_cbL_y, d_cbR_y, beta_coef_y, uniform_grid]')
else
$:GPU_UPDATE(device='[poly_coef_cbL_z, poly_coef_cbR_z, d_cbL_z, d_cbR_z, beta_coef_z]')
$:GPU_UPDATE(device='[poly_coef_cbL_z, poly_coef_cbR_z, d_cbL_z, d_cbR_z, beta_coef_z, uniform_grid]')
end if

! Nullifying WENO coefficients and cell-boundary locations pointers
Expand Down Expand Up @@ -1053,12 +1069,22 @@ contains
poly(2) = v_rs_ws_${XYZ}$ (j, k, l, i) + poly_coef_cbL_${XYZ}$ (j, 2, &
& 0)*dvd(-1) + poly_coef_cbL_${XYZ}$ (j, 2, 1)*dvd(-2)

beta(0) = beta_coef_${XYZ}$ (j, 0, 0)*dvd(1)*dvd(1) + beta_coef_${XYZ}$ (j, 0, &
& 1)*dvd(1)*dvd(0) + beta_coef_${XYZ}$ (j, 0, 2)*dvd(0)*dvd(0) + weno_eps
beta(1) = beta_coef_${XYZ}$ (j, 1, 0)*dvd(0)*dvd(0) + beta_coef_${XYZ}$ (j, 1, &
& 1)*dvd(0)*dvd(-1) + beta_coef_${XYZ}$ (j, 1, 2)*dvd(-1)*dvd(-1) + weno_eps
beta(2) = beta_coef_${XYZ}$ (j, 2, 0)*dvd(-1)*dvd(-1) + beta_coef_${XYZ}$ (j, 2, &
& 1)*dvd(-1)*dvd(-2) + beta_coef_${XYZ}$ (j, 2, 2)*dvd(-2)*dvd(-2) + weno_eps
! Jiang & Shu (1996) sum-of-squares form on uniform grids: all terms non-negative, no
! cancellation. On non-uniform grids, fall back to precomputed coefficients.
if (uniform_grid(${WENO_DIR}$)) then
beta(0) = 13._wp/12._wp*(dvd(1) - dvd(0))**2 + 0.25_wp*(dvd(1) - 3._wp*dvd(0))**2 &
& + weno_eps
beta(1) = 13._wp/12._wp*(dvd(0) - dvd(-1))**2 + 0.25_wp*(dvd(0) + dvd(-1))**2 + weno_eps
beta(2) = 13._wp/12._wp*(dvd(-1) - dvd(-2))**2 + 0.25_wp*(3._wp*dvd(-1) - dvd(-2))**2 &
& + weno_eps
else
beta(0) = beta_coef_${XYZ}$ (j, 0, 0)*dvd(1)*dvd(1) + beta_coef_${XYZ}$ (j, 0, &
& 1)*dvd(1)*dvd(0) + beta_coef_${XYZ}$ (j, 0, 2)*dvd(0)*dvd(0) + weno_eps
beta(1) = beta_coef_${XYZ}$ (j, 1, 0)*dvd(0)*dvd(0) + beta_coef_${XYZ}$ (j, 1, &
& 1)*dvd(0)*dvd(-1) + beta_coef_${XYZ}$ (j, 1, 2)*dvd(-1)*dvd(-1) + weno_eps
beta(2) = beta_coef_${XYZ}$ (j, 2, 0)*dvd(-1)*dvd(-1) + beta_coef_${XYZ}$ (j, 2, &
& 1)*dvd(-1)*dvd(-2) + beta_coef_${XYZ}$ (j, 2, 2)*dvd(-2)*dvd(-2) + weno_eps
end if

if (wenojs) then
do q = 0, weno_num_stencils
Expand Down
Loading