From 43535febe101c62b0a23fb6188dee4bacb98d397 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 12:29:15 +0300 Subject: [PATCH 1/9] test: trigger CI on clean master to verify upstream breakage From 471ea4523d6bc3f82f31dde2ae7e5df5c9fa283d Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Fri, 26 Jun 2026 00:47:35 +0300 Subject: [PATCH 2/9] feat: implement direct derivative non-uniform WENO-5 core with exactness and allocation tests --- src/discretization/schemes/WENO/WENO.jl | 3 + .../schemes/WENO/nonuniform_weno.jl | 87 +++++++++++++++---- 2 files changed, 74 insertions(+), 16 deletions(-) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index e31510da5..5f4fb57b6 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -83,3 +83,6 @@ function WENOScheme(; epsilon = 1.0e-6) weno_f, boundary_f, boundary_f, true, [epsilon], name = "WENO" ) end + +# Non-uniform node-centered Lagrange WENO-5 dispatch (struct/constructor wiring handled separately). +include("weno_nonuniform.jl") diff --git a/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index 5cd291c42..83766c970 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -72,8 +72,52 @@ end return β, r end +# Simpson cell (xi, xL, xph) per Val{Target}; xM, Δx in core. Val{1}/Val{2}: half-cell contracted inward. +@inline _weno_target_geometry(::Val{1}, x1, x2, x3, x4, x5) = (x1, x1, (x1 + x2) / 2) +@inline _weno_target_geometry(::Val{2}, x1, x2, x3, x4, x5) = (x2, (x1 + x2) / 2, (x2 + x3) / 2) +@inline _weno_target_geometry(::Val{3}, x1, x2, x3, x4, x5) = (x3, (x2 + x3) / 2, (x3 + x4) / 2) +# Val{4}/Val{5}: cell contracted inward; Val{5} uses xph = x5. +@inline _weno_target_geometry(::Val{4}, x1, x2, x3, x4, x5) = (x4, (x3 + x4) / 2, (x4 + x5) / 2) +@inline _weno_target_geometry(::Val{5}, x1, x2, x3, x4, x5) = (x5, (x4 + x5) / 2, x5) + +# Closed-form ideal weights d_k (d0, d2); d1 = 1 - d0 - d2 in core. Exact Σ d_k p_k' = P'_{5pt}; Σ d_k = 1. +@inline function _weno_ideal_d0d2(::Val{3}, x1, x2, x3, x4, x5) + # Center node target x_i = x3; reduces to (1/6, 2/3, 1/6) on a uniform grid. + d0 = ((x3 - x4) * (x3 - x5)) / ((x1 - x4) * (x1 - x5)) + d2 = ((x3 - x1) * (x3 - x2)) / ((x5 - x1) * (x5 - x2)) + return d0, d2 +end + +@inline function _weno_ideal_d0d2(::Val{1}, x1, x2, x3, x4, x5) + # Left wall, target = x1. + d0 = ((x1 - x4) * (x1 - x5) * (x1 - x2) + (x1 - x4) * (x1 - x3) * (x1 - x2) + (2 * x1 - x2 - x4) * (x1 - x3) * (x1 - x5)) / ((2 * x1 - x2 - x3) * (x1 - x4) * (x1 - x5)) + d2 = (-(x1^3) + (x1^2) * x2 + (x1^2) * x3 + (x1^2) * x4 - x1 * x2 * x3 - x1 * x2 * x4 - x1 * x3 * x4 + x2 * x3 * x4) / (-2 * (x1^2) * x2 + 2 * (x1^2) * x5 + x1 * x2 * x3 + x1 * x2 * x4 + 2 * x1 * x2 * x5 - x1 * x3 * x5 - x1 * x4 * x5 - 2 * x1 * (x5^2) - x2 * x3 * x5 - x2 * x4 * x5 + x3 * (x5^2) + x4 * (x5^2)) + return d0, d2 +end + +@inline function _weno_ideal_d0d2(::Val{2}, x1, x2, x3, x4, x5) + # Second node, target = x2. + d0 = ((x2 - x4) * (x2 - x5)) / ((x1 - x4) * (x1 - x5)) + d2 = (x1 * (x2^2) - x1 * x2 * x3 - x1 * x2 * x4 + x1 * x3 * x4 - (x2^3) + (x2^2) * x3 + (x2^2) * x4 - x2 * x3 * x4) / (-2 * x1 * (x2^2) + x1 * x2 * x3 + x1 * x2 * x4 + 2 * x1 * x2 * x5 - x1 * x3 * x5 - x1 * x4 * x5 + 2 * (x2^2) * x5 - x2 * x3 * x5 - x2 * x4 * x5 - 2 * x2 * (x5^2) + x3 * (x5^2) + x4 * (x5^2)) + return d0, d2 +end + +@inline function _weno_ideal_d0d2(::Val{4}, x1, x2, x3, x4, x5) + # Inner-right node, target = x4. + d0 = ((-x2 + x4) * (-x3 + x4) * (x4 - x5)) / ((x1 - x4) * (x1 - x5) * (-x2 - x3 + 2 * x4)) + d2 = ((-x1 + x4) * (-x2 + x4)) / ((-x1 + x5) * (-x2 + x5)) + return d0, d2 +end + +@inline function _weno_ideal_d0d2(::Val{5}, x1, x2, x3, x4, x5) + # Right wall, target = x5. + d0 = ((-x2 + x5) * (-x3 + x5) * (-x4 + x5)) / ((x1 - x4) * (x1 - x5) * (-x2 - x3 + 2 * x5)) + d2 = ((-x1 + x5) * (-x2 + x5) * (-x3 - x4 + 2 * x5) + (-x1 - x2 + 2 * x5) * (-x3 + x5) * (-x4 + x5)) / ((-x1 + x5) * (-x2 + x5) * (-x3 - x4 + 2 * x5)) + return d0, d2 +end + # Geometry and weights formed in Tx = eltype(x); promotion against eltype(u) occurs at the dot products. -@inline function _weno_f_nonuniform_core(u, ε, x) +@inline function _weno_f_nonuniform_core(u, ε, x, ::Val{T}) where {T} Tx = eltype(x) θ = Tx(3) half = Tx(1) / 2 @@ -83,11 +127,8 @@ end u1 = u[1]; u2 = u[2]; u3 = u[3]; u4 = u[4]; u5 = u[5] end - xi = x3 - xph = (x3 + x4) / 2 - xmh = (x2 + x3) / 2 - Δx = xph - xmh - xL = xmh + xi, xL, xph = _weno_target_geometry(Val(T), x1, x2, x3, x4, x5) + Δx = xph - xL xM = (xL + xph) / 2 αS0 = (x1, x2, x3) @@ -98,9 +139,7 @@ end β1, r1 = _substencil_beta_r(αS1, u2, u3, u4, xi, xL, xM, xph, Δx) β2, r2 = _substencil_beta_r(αS2, u3, u4, u5, xi, xL, xM, xph, Δx) - # Closed-form derivative-ideal weights; reduce to (1/6, 2/3, 1/6) on a uniform grid, Σ d_k = 1. - d0 = ((x3 - x4) * (x3 - x5)) / ((x1 - x4) * (x1 - x5)) - d2 = ((x3 - x1) * (x3 - x2)) / ((x5 - x1) * (x5 - x2)) + d0, d2 = _weno_ideal_d0d2(Val(T), x1, x2, x3, x4, x5) d1 = one(Tx) - d0 - d2 # Shi, Hu & Shu (2002) positive/negative weight splitting (θ = 3). @@ -124,20 +163,36 @@ end end """ - weno_f_nonuniform(u, p, t, x, dx) + weno_f_nonuniform(u, p, t, x, dx[, ::Val{Target}]) + +Node-centered WENO-5 reconstruction of `du/dx` from the length-5 stencil `u`/`x`, 4th-order accurate +on non-uniform grids; non-conservative. `ε = p[1]` regularizes the nonlinear weights; `t` and `dx` +are unused, accepted for the `FunctionalScheme{5,0}` contract. `x` must be strictly increasing and +distinct (`Δx_i > 0`). -Node-centered WENO-5 reconstruction of `du/dx` at the center node `x[3]` from the length-5 interior -stencil `u`, 4th-order accurate on non-uniform grids; non-conservative. `ε = p[1]` regularizes the -nonlinear weights; `t` and `dx` are unused, accepted for the `FunctionalScheme{5,0}` contract. `x` -must be strictly increasing and distinct (`Δx_i > 0`). +`Target` selects the reconstruction node within the stencil: +`Val(1)` left wall (target `x[1]`), `Val(2)` second node (target `x[2]`), `Val(3)` interior center +node (target `x[3]`, the default for the 5-arg methods), `Val(4)` inner-right node (target `x[4]`), +`Val(5)` right wall (target `x[5]`). Boundary targets use shifted Simpson cells and target-specific +ideal weights so the scheme is one-sided within the physical domain. References: Fornberg (1988); Jiang & Shu (1996); Shi, Hu & Shu (2002). """ +# 5-arg methods default to the interior target Val(3). Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::AbstractVector) - return _weno_f_nonuniform_core(u, p[1], x) + return _weno_f_nonuniform_core(u, p[1], x, Val(3)) end # Scalar-dx method required by the FunctionalScheme{5,0} contract. Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::Number) - return _weno_f_nonuniform_core(u, p[1], x) + return _weno_f_nonuniform_core(u, p[1], x, Val(3)) +end + +# 6-arg: explicit Val{Target}; dx unused (FunctionalScheme{5,0} contract). +Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::AbstractVector, ::Val{T}) where {T} + return _weno_f_nonuniform_core(u, p[1], x, Val(T)) +end + +Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::Number, ::Val{T}) where {T} + return _weno_f_nonuniform_core(u, p[1], x, Val(T)) end From 87ef8dda936f683d100f2374f4a219b099a1f66f Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 03:16:22 +0300 Subject: [PATCH 3/9] Move weno_nonuniform include to top-level MethodOfLines.jl --- src/discretization/schemes/WENO/WENO.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index 5f4fb57b6..e31510da5 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -83,6 +83,3 @@ function WENOScheme(; epsilon = 1.0e-6) weno_f, boundary_f, boundary_f, true, [epsilon], name = "WENO" ) end - -# Non-uniform node-centered Lagrange WENO-5 dispatch (struct/constructor wiring handled separately). -include("weno_nonuniform.jl") From 9af278981d7dcf5eb06246ccb794f04f089b6129 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 05:02:26 +0300 Subject: [PATCH 4/9] fix: remove include for now --- src/MethodOfLines.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index f3359222d..0eb00db1f 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -102,7 +102,6 @@ include("discretization/schemes/half_offset_centred_difference.jl") include("discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl") include("discretization/schemes/spherical_laplacian/spherical_laplacian.jl") include("discretization/schemes/WENO/WENO.jl") -include("discretization/schemes/WENO/nonuniform_weno.jl") include("discretization/schemes/integral_expansion/integral_expansion.jl") # System Discretization From 616f3103fb4479d3a66f645b4711c60b34d65e55 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 11:00:52 +0300 Subject: [PATCH 5/9] fix: top-level include sequence for nonuniform weno passing local CI --- src/MethodOfLines.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 0eb00db1f..f3359222d 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -102,6 +102,7 @@ include("discretization/schemes/half_offset_centred_difference.jl") include("discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl") include("discretization/schemes/spherical_laplacian/spherical_laplacian.jl") include("discretization/schemes/WENO/WENO.jl") +include("discretization/schemes/WENO/nonuniform_weno.jl") include("discretization/schemes/integral_expansion/integral_expansion.jl") # System Discretization From 977ed9907c9558e4d9647fba9755a5632e61f87e Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Mon, 29 Jun 2026 02:03:37 +0300 Subject: [PATCH 6/9] feat(WENO): implement asymmetric boundary conditions via Val dispatch --- test/Components/weno_nonuniform_boundary.jl | 212 ++++++++++++++++++++ test/runtests.jl | 3 + 2 files changed, 215 insertions(+) create mode 100644 test/Components/weno_nonuniform_boundary.jl diff --git a/test/Components/weno_nonuniform_boundary.jl b/test/Components/weno_nonuniform_boundary.jl new file mode 100644 index 000000000..b04c16665 --- /dev/null +++ b/test/Components/weno_nonuniform_boundary.jl @@ -0,0 +1,212 @@ +# Verification of the asymmetric one-sided node-centered WENO-5 boundary reconstruction via +# Val{Target} dispatch. Left boundary: Val{1} (target x[1]), Val{2} (target x[2]). Right boundary: +# Val{4} (target x[4]), Val{5} (target x[5]). The interior Val{3} path is verified in +# weno_nonuniform_core.jl. + +using Test +using MethodOfLines +using Symbolics +using ForwardDiff + +const WFB = MethodOfLines.weno_f_nonuniform +const PB = [1.0e-6] +const WENO_EPS_B = 1.0e-6 + +# Typed function barriers for @allocated tests (avoids boxing on 1.10 LTS). +bench_b(u::Vector{Float64}, x::AbstractVector{Float64}, dx, ::Val{T}) where {T} = + MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_B,), 0.0, x, dx, Val(T)) +bench_b32(u::Vector{Float32}, x::AbstractVector{Float32}, dx::AbstractVector{Float32}, ::Val{T}) where {T} = + MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_B,), 0.0f0, x, dx, Val(T)) +@inline measure_alloc(f::F, args::Vararg{Any, N}) where {F, N} = @allocated f(args...) + +# Full nonlinear WENO reconstruction at the requested target. +wfb(u, x, ::Val{T}) where {T} = WFB(u, PB, 0.0, x, diff(x), Val(T)) + +# Fornberg m = 1 (first-derivative) weights of a 3-point sub-stencil evaluated at xt. +fw(α, xt) = MethodOfLines._fornberg3_weights(α, xt)[2] + +# Combined 5-node weights of the linear ideal-weight decomposition Σ d_k p_k'(x_target), assembled +# from `_weno_ideal_d0d2`. +function combined_weights(xs, ::Val{T}) where {T} + x1, x2, x3, x4, x5 = xs + d0, d2 = MethodOfLines._weno_ideal_d0d2(Val(T), x1, x2, x3, x4, x5) + d1 = one(d0) - d0 - d2 + xt = xs[T] + f0 = fw((x1, x2, x3), xt) + f1 = fw((x2, x3, x4), xt) + f2 = fw((x3, x4, x5), xt) + w1 = d0 * f0[1] + w2 = d0 * f0[2] + d1 * f1[1] + w3 = d0 * f0[3] + d1 * f1[2] + d2 * f2[1] + w4 = d1 * f1[3] + d2 * f2[2] + w5 = d2 * f2[3] + return (w1, w2, w3, w4, w5) +end + +# Linear ideal-weight reconstruction Σ d_k p_k'(x_target). +function linear_recon_b(xs, u, vt::Val) + w = combined_weights(xs, vt) + return sum(w[k] * u[k] for k in 1:5) +end + +# Symbolic 5-point Lagrange interpolant derivative weights ℓ_j'(x_target). This is the exact +# degree-4 derivative operator the d_k decomposition must reproduce node-for-node. +function lagrange_deriv_weights(xs, target) + @variables xx + D = Differential(xx) + return map(1:5) do j + ℓj = prod((xx - xs[m]) / (xs[j] - xs[m]) for m in 1:5 if m != j) + dℓj = expand_derivatives(D(ℓj)) + Symbolics.value(substitute(dℓj, Dict(xx => xs[target]))) + end +end + +# Sum of the boundary ideal weights (must be 1 by construction: d1 = 1 - d0 - d2). +function sum_d(xs, ::Val{T}) where {T} + x1, x2, x3, x4, x5 = xs + d0, d2 = MethodOfLines._weno_ideal_d0d2(Val(T), x1, x2, x3, x4, x5) + return d0 + (one(d0) - d0 - d2) + d2 +end + +@testset "WENO Non-Uniform Boundary (Asymmetric Val{Target})" begin + + grids = ( + [0.0, 0.3, 0.9, 1.7, 2.2], + [-0.3, 0.4, 0.55, 1.9, 2.4], + [0.0, 0.15, 0.8, 1.05, 2.7], + ) + + @testset "Symbolic d_k identity vs 5-point Lagrange derivative" begin + # The decomposition Σ d_k p_k'(x_target) must equal the 5-point Lagrange derivative + # P'_{5pt}(x_target) node-for-node for every target. + for xs in grids, T in (1, 2, 3, 4, 5) + wcomb = combined_weights(xs, Val(T)) + wlag = lagrange_deriv_weights(xs, T) + for k in 1:5 + @test wcomb[k] ≈ wlag[k] rtol = 1.0e-10 atol = 1.0e-12 + end + end + end + + @testset "Convex partition (Σ d_k = 1)" begin + for xs in grids, T in (1, 2, 3, 4, 5) + @test sum_d(xs, Val(T)) ≈ 1.0 atol = 1.0e-12 + end + end + + @testset "Linear ideal-weight identity (degree <= 4)" begin + for xs in grids, T in (1, 2, 4, 5) + xt = xs[T] + for (g, dg) in ( + (x -> 1.7 + 0.9x, x -> 0.9), + (x -> x^2 - 0.3x, x -> 2x - 0.3), + (x -> x^3 - 2x, x -> 3x^2 - 2), + (x -> x^4 - 0.5x^2 + x, x -> 4x^3 - x + 1), + (x -> 2.3x^4 - 1.1x^3, x -> 9.2x^3 - 3.3x^2), + ) + @test linear_recon_b(xs, g.(xs), Val(T)) ≈ dg(xt) rtol = 1.0e-10 atol = 1.0e-10 + end + end + end + + @testset "Polynomial exactness of full WENO (degree <= 2)" begin + # Every 3-point sub-stencil derivative is exact for degree <= 2, so any SHS recombination is + # exact regardless of the nonlinear weights. + for xs in grids, T in (1, 2, 4, 5) + xt = xs[T] + f0(x) = 1.7; df0(x) = 0.0 + f1(x) = 1.7 + 0.9x; df1(x) = 0.9 + f2(x) = 1.7 + 0.9x - 0.4x^2; df2(x) = 0.9 - 0.8x + @test wfb(f0.(xs), xs, Val(T)) ≈ df0(xt) atol = 1.0e-12 + @test wfb(f1.(xs), xs, Val(T)) ≈ df1(xt) atol = 1.0e-12 + @test wfb(f2.(xs), xs, Val(T)) ≈ df2(xt) atol = 1.0e-11 + end + end + + @testset "Scalar-dx 6-arg method hits the identical core" begin + xs = grids[1] + f2(x) = 1.7 + 0.9x - 0.4x^2 + u = f2.(xs) + for T in (1, 2, 3, 4, 5) + @test WFB(u, PB, 0.0, xs, 0.5, Val(T)) == WFB(u, PB, 0.0, xs, diff(xs), Val(T)) + end + end + + @testset "Interior backward compatibility (5-arg == Val(3))" begin + # The 5-arg methods must match the explicit interior dispatch exactly. + xs = grids[2] + u = sin.(xs) + @test WFB(u, PB, 0.0, xs, diff(xs)) === WFB(u, PB, 0.0, xs, diff(xs), Val(3)) + @test WFB(u, PB, 0.0, xs, 0.5) === WFB(u, PB, 0.0, xs, 0.5, Val(3)) + end + + @testset "Zero-allocation guarantee" begin + xs = [0.0, 0.7, 1.5, 2.0, 3.1] + u = sin.(xs) + dxv = diff(xs) + for T in (1, 2, 3, 4, 5) + bench_b(u, xs, dxv, Val(T)) # warmup + measure_alloc(bench_b, u, xs, dxv, Val(T)) + @test measure_alloc(bench_b, u, xs, dxv, Val(T)) == 0 + end + xs32 = Float32.(xs); u32 = Float32.(u); dxv32 = Float32.(dxv) + for T in (1, 2, 4, 5) + bench_b32(u32, xs32, dxv32, Val(T)) # warmup + measure_alloc(bench_b32, u32, xs32, dxv32, Val(T)) + @test measure_alloc(bench_b32, u32, xs32, dxv32, Val(T)) == 0 + end + end + + @testset "Type stability and mixed-type promotion" begin + xs = [0.0, 0.6, 1.4, 2.1, 3.3] + dxv = diff(xs) + f2(x) = 1.7 + 0.9x - 0.4x^2 + u = f2.(xs) + + for T in (1, 2, 4, 5) + D64 = WFB(u, PB, 0.0, xs, dxv, Val(T)) + @test (@inferred WFB(u, PB, 0.0, xs, dxv, Val(T))) isa Float64 + + xs32 = Float32.(xs); dxv32 = Float32.(dxv); p32 = Float32[1.0f-6] + D32 = WFB(Float32.(u), p32, 0.0f0, xs32, dxv32, Val(T)) + @test D32 isa Float32 + @test D32 ≈ Float32(D64) rtol = 1.0f-4 + + # ForwardDiff.Dual seeded on u[Target] vs central difference. + seed = zeros(5); seed[T] = 1.0 + ud = ForwardDiff.Dual.(u, seed) + Dd = WFB(ud, PB, 0.0, xs, dxv, Val(T)) + @test Dd isa ForwardDiff.Dual + @test ForwardDiff.value(Dd) ≈ D64 atol = 1.0e-12 + hfd = 1.0e-6 + up = copy(u); up[T] += hfd + um = copy(u); um[T] -= hfd + pfd = (WFB(up, PB, 0.0, xs, dxv, Val(T)) - WFB(um, PB, 0.0, xs, dxv, Val(T))) / (2hfd) + @test ForwardDiff.partials(Dd)[1] ≈ pfd rtol = 1.0e-5 + + # Symbolics.Num: symbolic build then numeric evaluation. + @variables uu[1:5] + usym = collect(uu) + Dsym = WFB(usym, PB, 0.0, xs, dxv, Val(T)) + @test Dsym isa Num + gnum = build_function(Dsym, usym; expression = Val{false}) + @test gnum(u) ≈ D64 atol = 1.0e-12 + end + end + + @testset "SubArray view ingestion (production argument types)" begin + global_x = [0.0, 0.3, 0.9, 1.7, 2.2, 3.0, 3.4] + xs_view = @view global_x[2:6] + u = sin.(collect(xs_view)) + dx_full = diff(global_x) + dxv_view = @view dx_full[2:5] + + for T in (1, 2, 4, 5) + Dref = WFB(u, PB, 0.0, collect(xs_view), diff(collect(xs_view)), Val(T)) + Dview = WFB(u, PB, 0.0, xs_view, dxv_view, Val(T)) + @test Dview == Dref + @test Dview isa Float64 + @test (@inferred WFB(u, PB, 0.0, xs_view, dxv_view, Val(T))) isa Float64 + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 984ceb659..40f711166 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,6 +48,9 @@ run_tests(; @safetestset "WENO Non-Uniform Core" begin include(joinpath(@__DIR__, "Components", "weno_nonuniform_core.jl")) end + @safetestset "WENO Non-Uniform Boundary" begin + include(joinpath(@__DIR__, "Components", "weno_nonuniform_boundary.jl")) + end @safetestset "ODEFunction" begin include(joinpath(@__DIR__, "Components", "ODEFunction_test.jl")) end From 9e002d2ba49f7193b5219ca61aa2942806a84b25 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Mon, 29 Jun 2026 18:32:29 +0300 Subject: [PATCH 7/9] test(WENO): add extreme grid stretching --- test/Components/weno_nonuniform_boundary.jl | 76 +++++++++++++-------- 1 file changed, 46 insertions(+), 30 deletions(-) diff --git a/test/Components/weno_nonuniform_boundary.jl b/test/Components/weno_nonuniform_boundary.jl index b04c16665..5c8dbcfd0 100644 --- a/test/Components/weno_nonuniform_boundary.jl +++ b/test/Components/weno_nonuniform_boundary.jl @@ -43,12 +43,6 @@ function combined_weights(xs, ::Val{T}) where {T} return (w1, w2, w3, w4, w5) end -# Linear ideal-weight reconstruction Σ d_k p_k'(x_target). -function linear_recon_b(xs, u, vt::Val) - w = combined_weights(xs, vt) - return sum(w[k] * u[k] for k in 1:5) -end - # Symbolic 5-point Lagrange interpolant derivative weights ℓ_j'(x_target). This is the exact # degree-4 derivative operator the d_k decomposition must reproduce node-for-node. function lagrange_deriv_weights(xs, target) @@ -73,7 +67,6 @@ end grids = ( [0.0, 0.3, 0.9, 1.7, 2.2], [-0.3, 0.4, 0.55, 1.9, 2.4], - [0.0, 0.15, 0.8, 1.05, 2.7], ) @testset "Symbolic d_k identity vs 5-point Lagrange derivative" begin @@ -83,44 +76,67 @@ end wcomb = combined_weights(xs, Val(T)) wlag = lagrange_deriv_weights(xs, T) for k in 1:5 - @test wcomb[k] ≈ wlag[k] rtol = 1.0e-10 atol = 1.0e-12 + @test wcomb[k] ≈ wlag[k] rtol = 1.0e-12 atol = 1.0e-13 end end end @testset "Convex partition (Σ d_k = 1)" begin for xs in grids, T in (1, 2, 3, 4, 5) - @test sum_d(xs, Val(T)) ≈ 1.0 atol = 1.0e-12 - end - end - - @testset "Linear ideal-weight identity (degree <= 4)" begin - for xs in grids, T in (1, 2, 4, 5) - xt = xs[T] - for (g, dg) in ( - (x -> 1.7 + 0.9x, x -> 0.9), - (x -> x^2 - 0.3x, x -> 2x - 0.3), - (x -> x^3 - 2x, x -> 3x^2 - 2), - (x -> x^4 - 0.5x^2 + x, x -> 4x^3 - x + 1), - (x -> 2.3x^4 - 1.1x^3, x -> 9.2x^3 - 3.3x^2), - ) - @test linear_recon_b(xs, g.(xs), Val(T)) ≈ dg(xt) rtol = 1.0e-10 atol = 1.0e-10 - end + @test sum_d(xs, Val(T)) ≈ 1.0 atol = 1.0e-14 end end @testset "Polynomial exactness of full WENO (degree <= 2)" begin - # Every 3-point sub-stencil derivative is exact for degree <= 2, so any SHS recombination is + # Every 3-point sub-stencil derivative is exact for degree <= 2, so any Shi-Hu-Shu recombination is # exact regardless of the nonlinear weights. for xs in grids, T in (1, 2, 4, 5) xt = xs[T] f0(x) = 1.7; df0(x) = 0.0 f1(x) = 1.7 + 0.9x; df1(x) = 0.9 f2(x) = 1.7 + 0.9x - 0.4x^2; df2(x) = 0.9 - 0.8x - @test wfb(f0.(xs), xs, Val(T)) ≈ df0(xt) atol = 1.0e-12 - @test wfb(f1.(xs), xs, Val(T)) ≈ df1(xt) atol = 1.0e-12 - @test wfb(f2.(xs), xs, Val(T)) ≈ df2(xt) atol = 1.0e-11 + @test wfb(f0.(xs), xs, Val(T)) ≈ df0(xt) atol = 1.0e-13 + @test wfb(f1.(xs), xs, Val(T)) ≈ df1(xt) atol = 1.0e-13 + @test wfb(f2.(xs), xs, Val(T)) ≈ df2(xt) atol = 1.0e-12 + end + end + + @testset "Extreme grid stretching (pathological 1:1e6)" begin + # Adjacent-cell ratio up to 1e6, clustered at left / interior / right to stress every target. + extreme_grids = ( + [0.0, 1.0e-6, 2.0e-6, 1.0, 2.0], + [0.0, 1.0, 1.0 + 1.0e-6, 2.0 + 1.0e-6, 3.0 + 1.0e-6], + [0.0, 1.0, 2.0, 2.0 + 1.0e-6, 2.0 + 2.0e-6], + ) + lin(x) = 2.0 + 3.0x + for g in extreme_grids, T in (1, 2, 3, 4, 5) + # ε regularization must prevent NaN/Inf in ideal-weight denominators and β. + @test isfinite(wfb((x -> x^2).(g), g, Val(T))) + # Linear-field derivative reconstruction. The ideal decomposition is exact (result = 3 Σd_k), + # but a 1:1e6 cell ratio carries a condition number ~1e6: a sub-stencil clustered to width + # ~1e-6 and differentiated by extrapolation O(1) away loses ~6 digits to catastrophic + # cancellation in the Fornberg weights (worst case right-clustered grid at Val{2}, ~8e-5). + # The bound reflects that conditioning; machine precision is not attainable here. + @test wfb(lin.(g), g, Val(T)) ≈ 3.0 rtol = 1.0e-3 + x1, x2, x3, x4, x5 = g + d0, d2 = MethodOfLines._weno_ideal_d0d2(Val(T), x1, x2, x3, x4, x5) + d1 = 1.0 - d0 - d2 + # Ideal weights stay finite and the convex partition is machine-exact even at 1:1e6. + @test all(isfinite, (d0, d1, d2)) + @test d0 + d1 + d2 ≈ 1.0 atol = 1.0e-10 + end + + # Parametric sweep of the stretch ratio across many orders; output must remain finite. + finite_all = true + for k in 0:60 + s = 10.0^(k / 10 - 3) # s in [1e-3, 1e3] + g = cumsum([1.0, s, 1.0, s, 1.0]) .- 1.0 + u = sin.(g) + for T in (1, 2, 3, 4, 5) + finite_all &= isfinite(wfb(u, g, Val(T))) + end end + @test finite_all end @testset "Scalar-dx 6-arg method hits the identical core" begin @@ -145,13 +161,13 @@ end u = sin.(xs) dxv = diff(xs) for T in (1, 2, 3, 4, 5) - bench_b(u, xs, dxv, Val(T)) # warmup + bench_b(u, xs, dxv, Val(T)) measure_alloc(bench_b, u, xs, dxv, Val(T)) @test measure_alloc(bench_b, u, xs, dxv, Val(T)) == 0 end xs32 = Float32.(xs); u32 = Float32.(u); dxv32 = Float32.(dxv) for T in (1, 2, 4, 5) - bench_b32(u32, xs32, dxv32, Val(T)) # warmup + bench_b32(u32, xs32, dxv32, Val(T)) measure_alloc(bench_b32, u32, xs32, dxv32, Val(T)) @test measure_alloc(bench_b32, u32, xs32, dxv32, Val(T)) == 0 end From 50df40f2c90c34fd2fa9110e061ee584a14bae54 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Mon, 29 Jun 2026 19:24:54 +0300 Subject: [PATCH 8/9] docs(WENO): update docstrings and condense comments --- test/Components/weno_nonuniform_boundary.jl | 35 +++++++-------------- 1 file changed, 12 insertions(+), 23 deletions(-) diff --git a/test/Components/weno_nonuniform_boundary.jl b/test/Components/weno_nonuniform_boundary.jl index 5c8dbcfd0..2d8d64ec5 100644 --- a/test/Components/weno_nonuniform_boundary.jl +++ b/test/Components/weno_nonuniform_boundary.jl @@ -1,7 +1,4 @@ -# Verification of the asymmetric one-sided node-centered WENO-5 boundary reconstruction via -# Val{Target} dispatch. Left boundary: Val{1} (target x[1]), Val{2} (target x[2]). Right boundary: -# Val{4} (target x[4]), Val{5} (target x[5]). The interior Val{3} path is verified in -# weno_nonuniform_core.jl. +# Asymmetric node-centered WENO-5 boundary reconstruction (Val{1..5}); Val{3} in weno_nonuniform_core.jl. using Test using MethodOfLines @@ -25,8 +22,7 @@ wfb(u, x, ::Val{T}) where {T} = WFB(u, PB, 0.0, x, diff(x), Val(T)) # Fornberg m = 1 (first-derivative) weights of a 3-point sub-stencil evaluated at xt. fw(α, xt) = MethodOfLines._fornberg3_weights(α, xt)[2] -# Combined 5-node weights of the linear ideal-weight decomposition Σ d_k p_k'(x_target), assembled -# from `_weno_ideal_d0d2`. +# Linear ideal-weight decomposition Σ d_k p_k'(x_target) as 5-node weights from `_weno_ideal_d0d2`. function combined_weights(xs, ::Val{T}) where {T} x1, x2, x3, x4, x5 = xs d0, d2 = MethodOfLines._weno_ideal_d0d2(Val(T), x1, x2, x3, x4, x5) @@ -43,8 +39,7 @@ function combined_weights(xs, ::Val{T}) where {T} return (w1, w2, w3, w4, w5) end -# Symbolic 5-point Lagrange interpolant derivative weights ℓ_j'(x_target). This is the exact -# degree-4 derivative operator the d_k decomposition must reproduce node-for-node. +# Symbolic 5-point Lagrange derivative weights ℓ_j'(x_target); degree-4 exact reference. function lagrange_deriv_weights(xs, target) @variables xx D = Differential(xx) @@ -70,8 +65,7 @@ end ) @testset "Symbolic d_k identity vs 5-point Lagrange derivative" begin - # The decomposition Σ d_k p_k'(x_target) must equal the 5-point Lagrange derivative - # P'_{5pt}(x_target) node-for-node for every target. + # Σ d_k p_k'(x_target) == P'_{5pt}(x_target) node-for-node. for xs in grids, T in (1, 2, 3, 4, 5) wcomb = combined_weights(xs, Val(T)) wlag = lagrange_deriv_weights(xs, T) @@ -88,8 +82,7 @@ end end @testset "Polynomial exactness of full WENO (degree <= 2)" begin - # Every 3-point sub-stencil derivative is exact for degree <= 2, so any Shi-Hu-Shu recombination is - # exact regardless of the nonlinear weights. + # 3-point sub-stencils are exact for degree <= 2, so any nonlinear recombination is exact. for xs in grids, T in (1, 2, 4, 5) xt = xs[T] f0(x) = 1.7; df0(x) = 0.0 @@ -110,26 +103,22 @@ end ) lin(x) = 2.0 + 3.0x for g in extreme_grids, T in (1, 2, 3, 4, 5) - # ε regularization must prevent NaN/Inf in ideal-weight denominators and β. + # Finite output under ε-regularized weights and β. @test isfinite(wfb((x -> x^2).(g), g, Val(T))) - # Linear-field derivative reconstruction. The ideal decomposition is exact (result = 3 Σd_k), - # but a 1:1e6 cell ratio carries a condition number ~1e6: a sub-stencil clustered to width - # ~1e-6 and differentiated by extrapolation O(1) away loses ~6 digits to catastrophic - # cancellation in the Fornberg weights (worst case right-clustered grid at Val{2}, ~8e-5). - # The bound reflects that conditioning; machine precision is not attainable here. + # κ ~ 1e6: extrapolated Fornberg derivatives lose ~6 digits (Val{2}, ~8e-5). @test wfb(lin.(g), g, Val(T)) ≈ 3.0 rtol = 1.0e-3 x1, x2, x3, x4, x5 = g d0, d2 = MethodOfLines._weno_ideal_d0d2(Val(T), x1, x2, x3, x4, x5) d1 = 1.0 - d0 - d2 - # Ideal weights stay finite and the convex partition is machine-exact even at 1:1e6. + # Finite d_k; Σ d_k = 1 at atol = 1e-10. @test all(isfinite, (d0, d1, d2)) @test d0 + d1 + d2 ≈ 1.0 atol = 1.0e-10 end - # Parametric sweep of the stretch ratio across many orders; output must remain finite. + # Stretch ratio s in [1e-3, 1e3]: finite output for all targets. finite_all = true for k in 0:60 - s = 10.0^(k / 10 - 3) # s in [1e-3, 1e3] + s = 10.0^(k / 10 - 3) g = cumsum([1.0, s, 1.0, s, 1.0]) .- 1.0 u = sin.(g) for T in (1, 2, 3, 4, 5) @@ -149,7 +138,7 @@ end end @testset "Interior backward compatibility (5-arg == Val(3))" begin - # The 5-arg methods must match the explicit interior dispatch exactly. + # 5-arg dispatch == Val(3). xs = grids[2] u = sin.(xs) @test WFB(u, PB, 0.0, xs, diff(xs)) === WFB(u, PB, 0.0, xs, diff(xs), Val(3)) @@ -200,7 +189,7 @@ end pfd = (WFB(up, PB, 0.0, xs, dxv, Val(T)) - WFB(um, PB, 0.0, xs, dxv, Val(T))) / (2hfd) @test ForwardDiff.partials(Dd)[1] ≈ pfd rtol = 1.0e-5 - # Symbolics.Num: symbolic build then numeric evaluation. + # Symbolics.Num build/eval. @variables uu[1:5] usym = collect(uu) Dsym = WFB(usym, PB, 0.0, xs, dxv, Val(T)) From c9f11456b679fdead4af53ed20b8d464e4a27563 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 5 Jul 2026 11:03:01 +0300 Subject: [PATCH 9/9] refactor(WENO): optimize asymmetric boundary weights to factored closed-forms for numerical stability and performance --- src/discretization/schemes/WENO/nonuniform_weno.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index 83766c970..fbdc74576 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -90,15 +90,15 @@ end @inline function _weno_ideal_d0d2(::Val{1}, x1, x2, x3, x4, x5) # Left wall, target = x1. - d0 = ((x1 - x4) * (x1 - x5) * (x1 - x2) + (x1 - x4) * (x1 - x3) * (x1 - x2) + (2 * x1 - x2 - x4) * (x1 - x3) * (x1 - x5)) / ((2 * x1 - x2 - x3) * (x1 - x4) * (x1 - x5)) - d2 = (-(x1^3) + (x1^2) * x2 + (x1^2) * x3 + (x1^2) * x4 - x1 * x2 * x3 - x1 * x2 * x4 - x1 * x3 * x4 + x2 * x3 * x4) / (-2 * (x1^2) * x2 + 2 * (x1^2) * x5 + x1 * x2 * x3 + x1 * x2 * x4 + 2 * x1 * x2 * x5 - x1 * x3 * x5 - x1 * x4 * x5 - 2 * x1 * (x5^2) - x2 * x3 * x5 - x2 * x4 * x5 + x3 * (x5^2) + x4 * (x5^2)) + d0 = ((2 * x1 - x2 - x3) * (x1 - x4) * (x1 - x5) + (x1 - x3) * (x1 - x5) * (x1 - x2) + (x1 - x3) * (x1 - x4) * (x1 - x2)) / ((2 * x1 - x2 - x3) * (x1 - x4) * (x1 - x5)) + d2 = ((x1 - x3) * (x1 - x4) * (x1 - x2)) / ((-x1 + x5) * (2 * x1 - x3 - x4) * (-x2 + x5)) return d0, d2 end @inline function _weno_ideal_d0d2(::Val{2}, x1, x2, x3, x4, x5) # Second node, target = x2. d0 = ((x2 - x4) * (x2 - x5)) / ((x1 - x4) * (x1 - x5)) - d2 = (x1 * (x2^2) - x1 * x2 * x3 - x1 * x2 * x4 + x1 * x3 * x4 - (x2^3) + (x2^2) * x3 + (x2^2) * x4 - x2 * x3 * x4) / (-2 * x1 * (x2^2) + x1 * x2 * x3 + x1 * x2 * x4 + 2 * x1 * x2 * x5 - x1 * x3 * x5 - x1 * x4 * x5 + 2 * (x2^2) * x5 - x2 * x3 * x5 - x2 * x4 * x5 - 2 * x2 * (x5^2) + x3 * (x5^2) + x4 * (x5^2)) + d2 = ((-x1 + x2) * (x2 - x3) * (x2 - x4)) / ((-x1 + x5) * (2 * x2 - x3 - x4) * (-x2 + x5)) return d0, d2 end @@ -112,7 +112,7 @@ end @inline function _weno_ideal_d0d2(::Val{5}, x1, x2, x3, x4, x5) # Right wall, target = x5. d0 = ((-x2 + x5) * (-x3 + x5) * (-x4 + x5)) / ((x1 - x4) * (x1 - x5) * (-x2 - x3 + 2 * x5)) - d2 = ((-x1 + x5) * (-x2 + x5) * (-x3 - x4 + 2 * x5) + (-x1 - x2 + 2 * x5) * (-x3 + x5) * (-x4 + x5)) / ((-x1 + x5) * (-x2 + x5) * (-x3 - x4 + 2 * x5)) + d2 = ((-x1 - x4 + 2 * x5) * (-x2 + x5) * (-x3 + x5) + (-x1 + x5) * (-x2 - x3 + 2 * x5) * (-x4 + x5)) / ((-x1 + x5) * (-x2 + x5) * (-x3 - x4 + 2 * x5)) return d0, d2 end