feat(WENO): Boundary conditions for non-uniform grids

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 = ((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) * (x2 - x3) * (x2 - x4)) / ((-x1 + x5) * (2 * x2 - x3 - x4) * (-x2 + x5))
+    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 - 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
+
 # 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

test/Components/weno_nonuniform_boundary.jl

@@ -0,0 +1,217 @@
+# Asymmetric node-centered WENO-5 boundary reconstruction (Val{1..5}); Val{3} 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]
+
+# 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)
+    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
+
+# Symbolic 5-point Lagrange derivative weights ℓ_j'(x_target); degree-4 exact reference.
+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],
+    )
+
+    @testset "Symbolic d_k identity vs 5-point Lagrange derivative" begin
+        # Σ 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)
+            for k in 1:5
+                @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-14
+        end
+    end
+
+    @testset "Polynomial exactness of full WENO (degree <= 2)" begin
+        # 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
+            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-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)
+            # Finite output under ε-regularized weights and β.
+            @test isfinite(wfb((x -> x^2).(g), g, Val(T)))
+            # κ ~ 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
+            # 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
+
+        # 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)
+            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
+        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
+        # 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))
+        @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))
+            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))
+            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 build/eval.
+            @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

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