From e11e45d21eafe22b0adaf81ec5cef254ebb6233c Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Fri, 26 Jun 2026 00:47:35 +0300 Subject: [PATCH 1/9] feat: implement direct derivative non-uniform WENO-5 core with exactness and allocation tests --- Project.toml | 4 +- src/MethodOfLines.jl | 1 + src/discretization/schemes/WENO/WENO.jl | 3 + .../schemes/WENO/nonuniform_weno.jl | 167 ++++++++++++- test/Components/weno_nonuniform_core.jl | 222 ++++++++++++++++++ test/runtests.jl | 3 + 6 files changed, 395 insertions(+), 5 deletions(-) create mode 100644 test/Components/weno_nonuniform_core.jl diff --git a/Project.toml b/Project.toml index 83e54808e..b813c6465 100644 --- a/Project.toml +++ b/Project.toml @@ -26,6 +26,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Combinatorics = "1" DiffEqBase = "7.3" DomainSets = "0.7, 0.8" +ForwardDiff = "0.10, 1" IfElse = "0.1" Interpolations = "0.14, 0.15, 0.16" Latexify = "0.16" @@ -51,6 +52,7 @@ Test = "1" julia = "1.10" [extras] +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6" OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" @@ -62,4 +64,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "NonlinearSolve", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqSSPRK", "SafeTestsets", "SciMLTesting", "StableRNGs"] +test = ["Test", "ForwardDiff", "NonlinearSolve", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqSSPRK", "SafeTestsets", "SciMLTesting", "StableRNGs"] 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 diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index 4b37dfbfb..e0aa0d364 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -85,3 +85,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 461d4dcd2..77f41b954 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -1,6 +1,165 @@ -@noinline function weno_f_nonuniform(u, p, t, x, dx::AbstractVector) - throw(ArgumentError("WENO on non-uniform grids is not yet implemented.")) +# Node-centered Lagrange WENO-5 on non-uniform grids. The three 3-point sub-stencil derivative +# estimates p_k'(x_i) are combined with closed-form derivative-ideal weights d_k under a +# Shi-Hu-Shu positive/negative split. References: Fornberg (1988); Jiang & Shu (1996); +# Shi, Hu & Shu (2002). + +# Fornberg (1988) finite-difference weights for a 3-point stencil, returning the m = 0, 1, 2 +# derivative operators. Unrolled over SVectors for stack allocation; index map ν = 1,2,3 -> 0,1,2. +@inline function _fornberg3_weights(α::NTuple{3, T}, xt::T) where {T} + α0, α1, α2 = α + + # n = 0 base case. + n1m0 = one(T); n1m1 = zero(T); n1m2 = zero(T) + c1 = one(T) + + # n = 1: introduce α1. + c2 = α1 - α0 + r1 = c1 / c2 + tA = α1 - xt + a1m0 = (tA * n1m0) / c2 + a1m1 = (tA * n1m1 - n1m0) / c2 + a1m2 = (tA * n1m2 - 2 * n1m1) / c2 + sB = α0 - xt + a2m0 = r1 * (-(sB) * n1m0) + a2m1 = r1 * (n1m0 - sB * n1m1) + a2m2 = r1 * (2 * n1m1 - sB * n1m2) + c1 = c2 + + n1m0 = a1m0; n1m1 = a1m1; n1m2 = a1m2 + n2m0 = a2m0; n2m1 = a2m1; n2m2 = a2m2 + + # n = 2: introduce α2. + c2 = (α2 - α0) * (α2 - α1) + r2 = c1 / c2 + c3a = α2 - α0 + c3b = α2 - α1 + tA2 = α2 - xt + b1m0 = (tA2 * n1m0) / c3a + b1m1 = (tA2 * n1m1 - n1m0) / c3a + b1m2 = (tA2 * n1m2 - 2 * n1m1) / c3a + b2m0 = (tA2 * n2m0) / c3b + b2m1 = (tA2 * n2m1 - n2m0) / c3b + b2m2 = (tA2 * n2m2 - 2 * n2m1) / c3b + sB2 = α1 - xt + b3m0 = r2 * (-(sB2) * n2m0) + b3m1 = r2 * (n2m0 - sB2 * n2m1) + b3m2 = r2 * (2 * n2m1 - sB2 * n2m2) + + m0 = SVector{3, T}(b1m0, b2m0, b3m0) # values + m1 = SVector{3, T}(b1m1, b2m1, b3m1) # first derivative + m2 = SVector{3, T}(b1m2, b2m2, b3m2) # second derivative + return m0, m1, m2 +end + +@inline _dot3(w::SVector{3}, a, b, c) = w[1] * a + w[2] * b + w[3] * c + +# Smoothness indicator β_k (Jiang & Shu 1996) via Simpson's rule over cell i, exact for the +# degree-2 sub-stencil integrand, together with the sub-stencil derivative estimate r_k = p_k'(x_i). +@inline function _substencil_beta_r( + α::NTuple{3, T}, ua, ub, uc, xi, xL, xM, xph, Δx + ) where {T} + _, m1i, _ = _fornberg3_weights(α, xi) + _, m1L, _ = _fornberg3_weights(α, xL) + _, m1M, m2M = _fornberg3_weights(α, xM) + _, m1R, _ = _fornberg3_weights(α, xph) + + r = _dot3(m1i, ua, ub, uc) # p_k'(x_i) + pL = _dot3(m1L, ua, ub, uc) # p_k'(x_{i-1/2}) + pM = _dot3(m1M, ua, ub, uc) # p_k'(x_M) + pR = _dot3(m1R, ua, ub, uc) # p_k'(x_{i+1/2}) + pp = _dot3(m2M, ua, ub, uc) # p_k'' (constant on a 3-point stencil) + + I1 = (Δx / 6) * (pL^2 + 4 * pM^2 + pR^2) # ∫ (p')^2 dx (Simpson, exact) + I2 = Δx * pp^2 # ∫ (p'')^2 dx + val = Δx * I1 + Δx^3 * I2 + β = IfElse.ifelse(val < zero(val), zero(val), val) + return β, r +end + +# Zero-allocation core. Geometry and weights are formed in Tx = eltype(x); promotion against +# eltype(u) (Symbolics.Num, ForwardDiff.Dual, Float32) occurs at the dot products. +@inline function _weno_f_nonuniform_core(u, p, x) + Tx = eltype(x) + ε = p[1] + θ = Tx(3) + half = Tx(1) / 2 + + @inbounds begin + x1 = Tx(x[1]); x2 = Tx(x[2]); x3 = Tx(x[3]); x4 = Tx(x[4]); x5 = Tx(x[5]) + u1 = u[1]; u2 = u[2]; u3 = u[3]; u4 = u[4]; u5 = u[5] + end + + # Reconstruction target is the center node x_i = x3. The cell [x_{i-1/2}, x_{i+1/2}] and its + # Simpson nodes are used only by the smoothness indicators. + xi = x3 + xph = (x3 + x4) / 2 # x_{i+1/2} + xmh = (x2 + x3) / 2 # x_{i-1/2} + Δx = xph - xmh # cell width Δx_i + xL = xmh + xM = (xL + xph) / 2 + + αS0 = (x1, x2, x3) + αS1 = (x2, x3, x4) + αS2 = (x3, x4, x5) + + β0, r0 = _substencil_beta_r(αS0, u1, u2, u3, xi, xL, xM, xph, Δx) + β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: the unique d_k with Σ_k d_k p_k'(x_i) = P'_5(x_i), + # reducing to (1/6, 2/3, 1/6) on a uniform grid. Σ_k d_k = 1. + d0 = ((x3 - x4) * (x3 - x5)) / ((x1 - x4) * (x1 - x5)) + d2 = ((x3 - x1) * (x3 - x2)) / ((x5 - x1) * (x5 - x2)) + d1 = one(Tx) - d0 - d2 + + # Shi-Hu-Shu (2002) positive/negative split (θ = 3). In this node-centered direct-derivative + # topology the d_k remain non-negative, so the split is presently a defensive mechanism and a + # placeholder for future Taylor-series flux expansions that may introduce negative weights. The + # positive and negative branches are kept separate and recombined as σp·Rp − σm·Rm. + dp0 = half * (d0 + θ * abs(d0)); dp1 = half * (d1 + θ * abs(d1)); dp2 = half * (d2 + θ * abs(d2)) + dm0 = dp0 - d0; dm1 = dp1 - d1; dm2 = dp2 - d2 + σp = dp0 + dp1 + dp2 + σm = dm0 + dm1 + dm2 + + ap0 = (dp0 / σp) / (ε + β0)^2; ap1 = (dp1 / σp) / (ε + β1)^2; ap2 = (dp2 / σp) / (ε + β2)^2 + sp = ap0 + ap1 + ap2 + ωp0 = ap0 / sp; ωp1 = ap1 / sp; ωp2 = ap2 / sp + + am0 = (dm0 / σm) / (ε + β0)^2; am1 = (dm1 / σm) / (ε + β1)^2; am2 = (dm2 / σm) / (ε + β2)^2 + sm = am0 + am1 + am2 + ωm0 = am0 / sm; ωm1 = am1 / sm; ωm2 = am2 / sm + + Rp = ωp0 * r0 + ωp1 * r1 + ωp2 * r2 + Rm = ωm0 * r0 + ωm1 * r1 + ωm2 * r2 + + return σp * Rp - σm * Rm end -Base.@propagate_inbounds @inline weno_f_nonuniform(u, p, t, x, dx::Real) = - weno_f_uniform(u, p, t, x, dx) +""" + weno_f_nonuniform(u, p, t, x, dx) + +Node-centered Lagrange WENO-5 reconstruction of a first spatial derivative on a non-uniform grid. + +Given the length-5 interior stencil values `u` and the corresponding absolute node coordinates +`x`, this implementation reconstructs the direct spatial derivative at the center node `x[3]` to +achieve 4th-order accuracy on non-uniform grids. It is a non-conservative formulation. The +smoothness parameter `ε = p[1]` regularizes the nonlinear weights; `t` and `dx` are accepted to +satisfy the `FunctionalScheme{5,0}` contract but are unused, as all geometry is taken from `x`. + +The coordinates `x` must be strictly monotonically increasing and distinct (`Δx_i > 0`): the +formulation divides by node differences and is undefined on degenerate grids. + +References: Fornberg (1988); Jiang & Shu (1996); Shi, Hu & Shu (2002). +""" +@inline function weno_f_nonuniform(u, p, t, x, dx::AbstractVector) + @boundscheck (length(u) == 5 && length(x) == 5) || + throw(ArgumentError("weno_f_nonuniform requires a length-5 interior stencil (got length(u)=$(length(u)), length(x)=$(length(x))).")) + return _weno_f_nonuniform_core(u, p, x) +end + +# Scalar dx method: uniform-stepsize fallback required by the FunctionalScheme contract. +@inline function weno_f_nonuniform(u, p, t, x, dx::Number) + @boundscheck (length(u) == 5 && length(x) == 5) || + throw(ArgumentError("weno_f_nonuniform requires a length-5 interior stencil (got length(u)=$(length(u)), length(x)=$(length(x))).")) + return _weno_f_nonuniform_core(u, p, x) +end diff --git a/test/Components/weno_nonuniform_core.jl b/test/Components/weno_nonuniform_core.jl new file mode 100644 index 000000000..e5b4a28d0 --- /dev/null +++ b/test/Components/weno_nonuniform_core.jl @@ -0,0 +1,222 @@ +# Verification of the node-centered direct-derivative WENO-5 core (`weno_f_nonuniform`, +# FunctionalScheme{5,0}, returning du/dx at the center node x[3]). +# +# Properties under test: the nonlinear scheme is exact through degree 2 (ω_k = d_k there) while the +# linear ideal-weight decomposition is exact through degree 4; the derivative-ideal weights d_k form +# a convex partition (d_k >= 0, Σ d_k = 1) on every monotone grid, so the Shi-Hu-Shu negative branch +# is inert and is verified through convexity rather than a negative-weight path. + +using Test +using MethodOfLines +using Symbolics +using ForwardDiff + +const WF = MethodOfLines.weno_f_nonuniform +const P = [1.0e-6] + +# Public-API helpers (vector / scalar dx). `dx` is unused by the kernel but required by the contract. +wf(u, x) = WF(u, P, 0.0, x, diff(x)) +wf_scalar(u, x) = WF(u, P, 0.0, x, 0.5) + +# Closed-form derivative-ideal weights (mirror of the core; the oracle in the linear-identity set is +# the analytic derivative, so a faulty formula fails against it). +function ideal_weights(x) + x1, x2, x3, x4, x5 = x + d0 = ((x3 - x4) * (x3 - x5)) / ((x1 - x4) * (x1 - x5)) + d2 = ((x3 - x1) * (x3 - x2)) / ((x5 - x1) * (x5 - x2)) + d1 = 1 - d0 - d2 + return (d0, d1, d2) +end + +# 3-point sub-stencil first derivative at xt via the shipped Fornberg engine (m = 1 weights). +function sub_deriv(α, ua, ub, uc, xt) + w1 = MethodOfLines._fornberg3_weights(α, xt)[2] + return w1[1] * ua + w1[2] * ub + w1[3] * uc +end + +# Linear (ideal-weight) reconstruction Σ d_k p_k'(x_i) assembled from the shipped pieces. +function linear_recon(x, u) + d0, d1, d2 = ideal_weights(x) + r0 = sub_deriv((x[1], x[2], x[3]), u[1], u[2], u[3], x[3]) + r1 = sub_deriv((x[2], x[3], x[4]), u[2], u[3], u[4], x[3]) + r2 = sub_deriv((x[3], x[4], x[5]), u[3], u[4], u[5], x[3]) + return d0 * r0 + d1 * r1 + d2 * r2 +end + +@testset "WENO Non-Uniform Core (Direct Derivative Reconstruction)" begin + + @testset "Polynomial exactness (degree <= 2)" begin + xs = [0.0, 0.6, 1.4, 2.1, 3.3] # ratios 0.6:0.8:0.7:1.2, center xs[3] = 1.4 + xc = xs[3] + + 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 wf(f0.(xs), xs) ≈ df0(xc) atol = 1e-14 + @test wf(f1.(xs), xs) ≈ df1(xc) atol = 1e-13 + @test wf(f2.(xs), xs) ≈ df2(xc) atol = 1e-12 + + # Scalar-dx fallback hits the identical core: bit-for-bit equal. + @test wf_scalar(f2.(xs), xs) == wf(f2.(xs), xs) + + # The nonlinear scheme is not exact for degree 3 at finite h (ω_k != d_k). + f3(x) = x^3 + @test !isapprox(wf(f3.(xs), xs), 3xc^2; atol = 1e-6) + end + + @testset "Linear ideal-weight identity (degree <= 4)" begin + for xs in ([0.0, 0.6, 1.4, 2.1, 3.3], [-0.3, 0.4, 0.55, 1.9, 2.4]) + xc = xs[3] + # Convex partition: Σ d_k = 1 and d_k in [0, 1] (the SHS negative branch is inert). + d0, d1, d2 = ideal_weights(xs) + @test d0 + d1 + d2 ≈ 1.0 atol = 1e-14 + @test all(0 .<= (d0, d1, d2) .<= 1) + + for (g, dg) in ( + (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(xs, g.(xs)) ≈ dg(xc) rtol = 1e-11 atol = 1e-11 + end + end + end + + @testset "Order of convergence (MMS)" begin + f(x) = sin(1.3x) + 0.5x + df(x) = 1.3cos(1.3x) + 0.5 + + # Frozen offsets: inner nodes perturbed ~10-20% off uniform, ends fixed. The relative + # geometry is constant under refinement (self-similar), so the asymptotic order is isolated. + o = [-2.0, -1.13, 0.08, 0.91, 2.0] + xc = 1.0 + hs = [0.2, 0.1, 0.05, 0.025, 0.0125] + + errs = map(hs) do h + xs = xc .+ h .* o + abs(wf(f.(xs), xs) - df(xs[3])) # evaluated at the actual center node xs[3] + end + orders = [log2(errs[k] / errs[k+1]) for k in 1:length(hs)-1] + + # Asymptotic order -> 4 (this FD scheme is 4th order on non-uniform grids). + @test orders[end] > 3.85 + @test orders[end-1] > 3.7 + # Least-squares slope over the finest region (h <= 0.05) confirms ~4th order. + lh = log.(hs[3:end]); le = log.(errs[3:end]); n = length(lh) + slope = (n * sum(lh .* le) - sum(lh) * sum(le)) / (n * sum(lh .^ 2) - sum(lh)^2) + @test slope > 3.8 + # A 1st/2nd/3rd-order regression would drag every order well below this. + @test all(>(3.0), orders) + end + + @testset "Extreme grid stretching (1:1e6)" begin + grids = ( + [0.0, 1.0, 1.0 + 1e-6, 2.0 + 1e-6, 3.0 + 1e-6], # 1e-6 center cell + [0.0, 1e-6, 2e-6, 1.0, 2.0], # extreme clustering near the center + ) + for g in grids + # Finiteness under severe ill-conditioning. + @test isfinite(wf((x -> x^2).(g), g)) + + # A linear field's derivative is exact regardless of the grid (Σ d_k = 1, every r_k + # equals the slope), so the SHS recombination must not corrupt it. + lin(x) = 2.0 + 3.0x + @test wf(lin.(g), g) ≈ 3.0 rtol = 1e-6 + + # Stability (not accuracy) target for a quadratic on a brutally stretched stencil. + @test abs(wf((x -> x^2).(g), g) - 2 * g[3]) < 1.0 + + # d_k stay a convex partition even at 1:1e6, so d1 < 0 is unreachable. + d0, d1, d2 = ideal_weights(g) + @test all(isfinite, (d0, d1, d2)) + @test d1 >= 0 + @test d0 + d1 + d2 ≈ 1.0 atol = 1e-10 + end + + # A batch of deterministic extreme-ratio grids must never produce NaN/Inf. + finite_all = true + for k in 0:60 + s = 10.0 ^ (k / 10 - 3) # cell-ratio scale spanning 1e-3 .. 1e3 + g = cumsum([1.0, s, 1.0, s, 1.0]) .- 1.0 + finite_all &= isfinite(wf((x -> sin(x)).(g), g)) + end + @test finite_all + end + + @testset "Zero-allocation guarantee" begin + xs = [0.0, 0.7, 1.5, 2.0, 3.1] + u = sin.(xs) + dxv = diff(xs) + + bench(uu, x, dx) = WF(uu, P, 0.0, x, dx) + bench(u, xs, dxv) # warmup (compile) + bench(u, xs, 0.5) + + @test @allocated(bench(u, xs, dxv)) == 0 + @test @allocated(bench(u, xs, 0.5)) == 0 + + xs32 = Float32.(xs); u32 = Float32.(u); dxv32 = Float32.(dxv); p32 = Float32[1.0f-6] + bench32(uu, x, dx) = WF(uu, p32, 0.0f0, x, dx) + bench32(u32, xs32, dxv32) + @test @allocated(bench32(u32, xs32, dxv32)) == 0 + end + + @testset "Type stability" 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) + D64 = WF(u, P, 0.0, xs, dxv) + + # Static inference: concrete Float64 return. + @test (@inferred WF(u, P, 0.0, xs, dxv)) isa Float64 + + # Float32 end-to-end. + xs32 = Float32.(xs); dxv32 = Float32.(dxv); p32 = Float32[1.0f-6] + D32 = WF(Float32.(u), p32, 0.0f0, xs32, dxv32) + @test D32 isa Float32 + @test D32 ≈ Float32(D64) rtol = 1.0f-4 + + # ForwardDiff.Dual: differentiate the output w.r.t. u[3]; compare to a central difference. + seed = [0.0, 0.0, 1.0, 0.0, 0.0] + ud = ForwardDiff.Dual.(u, seed) + Dd = WF(ud, P, 0.0, xs, dxv) + @test Dd isa ForwardDiff.Dual + @test ForwardDiff.value(Dd) ≈ D64 atol = 1e-12 + hfd = 1e-6 + up = copy(u); up[3] += hfd + um = copy(u); um[3] -= hfd + pfd = (WF(up, P, 0.0, xs, dxv) - WF(um, P, 0.0, xs, dxv)) / (2hfd) + @test ForwardDiff.partials(Dd)[1] ≈ pfd rtol = 1e-5 + + # Symbolics.Num: symbolic build then numeric evaluation must equal the direct numeric call. + @variables uu[1:5] + usym = collect(uu) + Dsym = WF(usym, P, 0.0, xs, dxv) + @test Dsym isa Num + gnum = build_function(Dsym, usym; expression = Val{false}) + @test gnum(u) ≈ D64 atol = 1e-12 + end + + @testset "SubArray view ingestion (production argument types)" begin + # The discretizer passes x as a `@view s.grid[x][itap]` and the non-uniform dx as a view of + # the spacing vector. Verify that compiler views preserve the result, inference, and zero + # allocations, rather than only plain `Vector`s. + global_x = [0.0, 0.3, 0.9, 1.7, 2.2, 3.0, 3.4] # length-7 backing grid + xs_view = @view global_x[2:6] # length-5 SubArray (the interior stencil) + u = sin.(collect(xs_view)) # discvars indexing materializes a Vector + dx_full = diff(global_x) + dxv_view = @view dx_full[2:5] # length-4 SubArray of spacings + + Dref = WF(u, P, 0.0, collect(xs_view), diff(collect(xs_view))) + Dview = WF(u, P, 0.0, xs_view, dxv_view) + @test Dview == Dref + @test Dview isa Float64 + @test (@inferred WF(u, P, 0.0, xs_view, dxv_view)) isa Float64 + + benchv(uu, x, dx) = WF(uu, P, 0.0, x, dx) + benchv(u, xs_view, dxv_view) # warmup + @test @allocated(benchv(u, xs_view, dxv_view)) == 0 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index a578097c9..984ceb659 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,6 +45,9 @@ run_tests(; @safetestset "WENO dispatch" begin include(joinpath(@__DIR__, "Components", "weno_dispatch.jl")) end + @safetestset "WENO Non-Uniform Core" begin + include(joinpath(@__DIR__, "Components", "weno_nonuniform_core.jl")) + end @safetestset "ODEFunction" begin include(joinpath(@__DIR__, "Components", "ODEFunction_test.jl")) end From 927eb3c41fa194a053d3707de89bf10bf082c478 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Fri, 26 Jun 2026 14:04:48 +0300 Subject: [PATCH 2/9] Fix 16-byte allocation on Julia LTS for non-uniform WENO --- .../schemes/WENO/nonuniform_weno.jl | 20 ++++++------ test/Components/weno_nonuniform_core.jl | 32 ++++++++++++------- 2 files changed, 29 insertions(+), 23 deletions(-) diff --git a/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index 77f41b954..0411e4758 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -53,6 +53,9 @@ end @inline _dot3(w::SVector{3}, a, b, c) = w[1] * a + w[2] * b + w[3] * c +@inline _beta_nonneg(val::T) where {T<:AbstractFloat} = max(val, zero(val)) +@inline _beta_nonneg(val) = IfElse.ifelse(val < zero(val), zero(val), val) + # Smoothness indicator β_k (Jiang & Shu 1996) via Simpson's rule over cell i, exact for the # degree-2 sub-stencil integrand, together with the sub-stencil derivative estimate r_k = p_k'(x_i). @inline function _substencil_beta_r( @@ -72,15 +75,14 @@ end I1 = (Δx / 6) * (pL^2 + 4 * pM^2 + pR^2) # ∫ (p')^2 dx (Simpson, exact) I2 = Δx * pp^2 # ∫ (p'')^2 dx val = Δx * I1 + Δx^3 * I2 - β = IfElse.ifelse(val < zero(val), zero(val), val) + β = _beta_nonneg(val) return β, r end # Zero-allocation core. Geometry and weights are formed in Tx = eltype(x); promotion against # eltype(u) (Symbolics.Num, ForwardDiff.Dual, Float32) occurs at the dot products. -@inline function _weno_f_nonuniform_core(u, p, x) +@inline function _weno_f_nonuniform_core(u, ε, x) Tx = eltype(x) - ε = p[1] θ = Tx(3) half = Tx(1) / 2 @@ -151,15 +153,11 @@ formulation divides by node differences and is undefined on degenerate grids. References: Fornberg (1988); Jiang & Shu (1996); Shi, Hu & Shu (2002). """ -@inline function weno_f_nonuniform(u, p, t, x, dx::AbstractVector) - @boundscheck (length(u) == 5 && length(x) == 5) || - throw(ArgumentError("weno_f_nonuniform requires a length-5 interior stencil (got length(u)=$(length(u)), length(x)=$(length(x))).")) - return _weno_f_nonuniform_core(u, p, x) +Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::AbstractVector) + return _weno_f_nonuniform_core(u, p[1], x) end # Scalar dx method: uniform-stepsize fallback required by the FunctionalScheme contract. -@inline function weno_f_nonuniform(u, p, t, x, dx::Number) - @boundscheck (length(u) == 5 && length(x) == 5) || - throw(ArgumentError("weno_f_nonuniform requires a length-5 interior stencil (got length(u)=$(length(u)), length(x)=$(length(x))).")) - return _weno_f_nonuniform_core(u, p, x) +Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::Number) + return _weno_f_nonuniform_core(u, p[1], x) end diff --git a/test/Components/weno_nonuniform_core.jl b/test/Components/weno_nonuniform_core.jl index e5b4a28d0..48fd49d98 100644 --- a/test/Components/weno_nonuniform_core.jl +++ b/test/Components/weno_nonuniform_core.jl @@ -13,6 +13,16 @@ using ForwardDiff const WF = MethodOfLines.weno_f_nonuniform const P = [1.0e-6] +const WENO_EPS_F64 = 1.0e-6 +const WENO_EPS_F32 = 1.0f-6 + +# Typed function barriers for @allocated tests (avoid closure boxing on LTS). +bench_weno_f64(u::Vector{Float64}, x::AbstractVector{Float64}, dx) = + MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, x, dx) +bench_weno_f32(u::Vector{Float32}, x::AbstractVector{Float32}, dx::AbstractVector{Float32}) = + MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F32,), 0.0f0, x, dx) +bench_weno_sub(u::Vector{Float64}, x::SubArray{Float64,1}, dx::SubArray{Float64,1}) = + MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, x, dx) # Public-API helpers (vector / scalar dx). `dx` is unused by the kernel but required by the contract. wf(u, x) = WF(u, P, 0.0, x, diff(x)) @@ -149,17 +159,16 @@ end u = sin.(xs) dxv = diff(xs) - bench(uu, x, dx) = WF(uu, P, 0.0, x, dx) - bench(u, xs, dxv) # warmup (compile) - bench(u, xs, 0.5) + bench_weno_f64(u, xs, dxv) # warmup (compile) + bench_weno_f64(u, xs, 0.5) - @test @allocated(bench(u, xs, dxv)) == 0 - @test @allocated(bench(u, xs, 0.5)) == 0 + @test @allocated(bench_weno_f64(u, xs, dxv)) == 0 + @test @allocated(bench_weno_f64(u, xs, 0.5)) == 0 + @test @allocated(MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, xs, dxv)) == 0 - xs32 = Float32.(xs); u32 = Float32.(u); dxv32 = Float32.(dxv); p32 = Float32[1.0f-6] - bench32(uu, x, dx) = WF(uu, p32, 0.0f0, x, dx) - bench32(u32, xs32, dxv32) - @test @allocated(bench32(u32, xs32, dxv32)) == 0 + xs32 = Float32.(xs); u32 = Float32.(u); dxv32 = Float32.(dxv) + bench_weno_f32(u32, xs32, dxv32) + @test @allocated(bench_weno_f32(u32, xs32, dxv32)) == 0 end @testset "Type stability" begin @@ -215,8 +224,7 @@ end @test Dview isa Float64 @test (@inferred WF(u, P, 0.0, xs_view, dxv_view)) isa Float64 - benchv(uu, x, dx) = WF(uu, P, 0.0, x, dx) - benchv(u, xs_view, dxv_view) # warmup - @test @allocated(benchv(u, xs_view, dxv_view)) == 0 + bench_weno_sub(u, xs_view, dxv_view) # warmup + @test @allocated(bench_weno_sub(u, xs_view, dxv_view)) == 0 end end From 8adf592bed562fe565cd073c6ae23a0669cbdd5c Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Fri, 26 Jun 2026 14:58:27 +0300 Subject: [PATCH 3/9] test: fix allocation artifacts and apply strict Runic formatting to WENO core --- .../schemes/WENO/nonuniform_weno.jl | 2 +- test/Components/weno_nonuniform_core.jl | 67 +++++++++++-------- 2 files changed, 39 insertions(+), 30 deletions(-) diff --git a/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index 0411e4758..bc215b262 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -53,7 +53,7 @@ end @inline _dot3(w::SVector{3}, a, b, c) = w[1] * a + w[2] * b + w[3] * c -@inline _beta_nonneg(val::T) where {T<:AbstractFloat} = max(val, zero(val)) +@inline _beta_nonneg(val::T) where {T <: AbstractFloat} = max(val, zero(val)) @inline _beta_nonneg(val) = IfElse.ifelse(val < zero(val), zero(val), val) # Smoothness indicator β_k (Jiang & Shu 1996) via Simpson's rule over cell i, exact for the diff --git a/test/Components/weno_nonuniform_core.jl b/test/Components/weno_nonuniform_core.jl index 48fd49d98..c9bf99aa3 100644 --- a/test/Components/weno_nonuniform_core.jl +++ b/test/Components/weno_nonuniform_core.jl @@ -21,7 +21,7 @@ bench_weno_f64(u::Vector{Float64}, x::AbstractVector{Float64}, dx) = MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, x, dx) bench_weno_f32(u::Vector{Float32}, x::AbstractVector{Float32}, dx::AbstractVector{Float32}) = MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F32,), 0.0f0, x, dx) -bench_weno_sub(u::Vector{Float64}, x::SubArray{Float64,1}, dx::SubArray{Float64,1}) = +bench_weno_sub(u::Vector{Float64}, x::SubArray{Float64, 1}, dx::SubArray{Float64, 1}) = MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, x, dx) # Public-API helpers (vector / scalar dx). `dx` is unused by the kernel but required by the contract. @@ -63,32 +63,32 @@ end 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 wf(f0.(xs), xs) ≈ df0(xc) atol = 1e-14 - @test wf(f1.(xs), xs) ≈ df1(xc) atol = 1e-13 - @test wf(f2.(xs), xs) ≈ df2(xc) atol = 1e-12 + @test wf(f0.(xs), xs) ≈ df0(xc) atol = 1.0e-14 + @test wf(f1.(xs), xs) ≈ df1(xc) atol = 1.0e-13 + @test wf(f2.(xs), xs) ≈ df2(xc) atol = 1.0e-12 # Scalar-dx fallback hits the identical core: bit-for-bit equal. @test wf_scalar(f2.(xs), xs) == wf(f2.(xs), xs) # The nonlinear scheme is not exact for degree 3 at finite h (ω_k != d_k). f3(x) = x^3 - @test !isapprox(wf(f3.(xs), xs), 3xc^2; atol = 1e-6) + @test !isapprox(wf(f3.(xs), xs), 3xc^2; atol = 1.0e-6) end @testset "Linear ideal-weight identity (degree <= 4)" begin for xs in ([0.0, 0.6, 1.4, 2.1, 3.3], [-0.3, 0.4, 0.55, 1.9, 2.4]) xc = xs[3] - # Convex partition: Σ d_k = 1 and d_k in [0, 1] (the SHS negative branch is inert). + # Convex partition: Σ d_k = 1 and d_k in [0, 1] (the Shi-Hu-Shu negative branch is inert). d0, d1, d2 = ideal_weights(xs) - @test d0 + d1 + d2 ≈ 1.0 atol = 1e-14 + @test d0 + d1 + d2 ≈ 1.0 atol = 1.0e-14 @test all(0 .<= (d0, d1, d2) .<= 1) for (g, dg) in ( - (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), + (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(xs, g.(xs)) ≈ dg(xc) rtol = 1e-11 atol = 1e-11 + @test linear_recon(xs, g.(xs)) ≈ dg(xc) rtol = 1.0e-11 atol = 1.0e-11 end end end @@ -107,11 +107,11 @@ end xs = xc .+ h .* o abs(wf(f.(xs), xs) - df(xs[3])) # evaluated at the actual center node xs[3] end - orders = [log2(errs[k] / errs[k+1]) for k in 1:length(hs)-1] + orders = [log2(errs[k] / errs[k + 1]) for k in 1:(length(hs) - 1)] # Asymptotic order -> 4 (this FD scheme is 4th order on non-uniform grids). @test orders[end] > 3.85 - @test orders[end-1] > 3.7 + @test orders[end - 1] > 3.7 # Least-squares slope over the finest region (h <= 0.05) confirms ~4th order. lh = log.(hs[3:end]); le = log.(errs[3:end]); n = length(lh) slope = (n * sum(lh .* le) - sum(lh) * sum(le)) / (n * sum(lh .^ 2) - sum(lh)^2) @@ -122,17 +122,17 @@ end @testset "Extreme grid stretching (1:1e6)" begin grids = ( - [0.0, 1.0, 1.0 + 1e-6, 2.0 + 1e-6, 3.0 + 1e-6], # 1e-6 center cell - [0.0, 1e-6, 2e-6, 1.0, 2.0], # extreme clustering near the center + [0.0, 1.0, 1.0 + 1.0e-6, 2.0 + 1.0e-6, 3.0 + 1.0e-6], # 1e-6 center cell + [0.0, 1.0e-6, 2.0e-6, 1.0, 2.0], # extreme clustering near the center ) for g in grids # Finiteness under severe ill-conditioning. @test isfinite(wf((x -> x^2).(g), g)) # A linear field's derivative is exact regardless of the grid (Σ d_k = 1, every r_k - # equals the slope), so the SHS recombination must not corrupt it. + # equals the slope), so the Shi-Hu-Shu recombination must not corrupt it. lin(x) = 2.0 + 3.0x - @test wf(lin.(g), g) ≈ 3.0 rtol = 1e-6 + @test wf(lin.(g), g) ≈ 3.0 rtol = 1.0e-6 # Stability (not accuracy) target for a quadratic on a brutally stretched stencil. @test abs(wf((x -> x^2).(g), g) - 2 * g[3]) < 1.0 @@ -141,13 +141,13 @@ end d0, d1, d2 = ideal_weights(g) @test all(isfinite, (d0, d1, d2)) @test d1 >= 0 - @test d0 + d1 + d2 ≈ 1.0 atol = 1e-10 + @test d0 + d1 + d2 ≈ 1.0 atol = 1.0e-10 end # A batch of deterministic extreme-ratio grids must never produce NaN/Inf. finite_all = true for k in 0:60 - s = 10.0 ^ (k / 10 - 3) # cell-ratio scale spanning 1e-3 .. 1e3 + s = 10.0^(k / 10 - 3) # cell-ratio scale spanning 1e-3 .. 1e3 g = cumsum([1.0, s, 1.0, s, 1.0]) .- 1.0 finite_all &= isfinite(wf((x -> sin(x)).(g), g)) end @@ -159,16 +159,23 @@ end u = sin.(xs) dxv = diff(xs) - bench_weno_f64(u, xs, dxv) # warmup (compile) + # Warmup (compile) every form exercised under @allocated, including the direct invocation. + bench_weno_f64(u, xs, dxv) bench_weno_f64(u, xs, 0.5) + MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, xs, dxv) - @test @allocated(bench_weno_f64(u, xs, dxv)) == 0 - @test @allocated(bench_weno_f64(u, xs, 0.5)) == 0 - @test @allocated(MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, xs, dxv)) == 0 + # `let` blocks isolate the locals from the @testset scope so LTS does not box captured globals. + let u_loc = u, xs_loc = xs, dxv_loc = dxv + @test @allocated(bench_weno_f64(u_loc, xs_loc, dxv_loc)) == 0 + @test @allocated(bench_weno_f64(u_loc, xs_loc, 0.5)) == 0 + @test @allocated(MethodOfLines.weno_f_nonuniform(u_loc, (WENO_EPS_F64,), 0.0, xs_loc, dxv_loc)) == 0 + end xs32 = Float32.(xs); u32 = Float32.(u); dxv32 = Float32.(dxv) bench_weno_f32(u32, xs32, dxv32) - @test @allocated(bench_weno_f32(u32, xs32, dxv32)) == 0 + let u_loc = u32, xs_loc = xs32, dxv_loc = dxv32 + @test @allocated(bench_weno_f32(u_loc, xs_loc, dxv_loc)) == 0 + end end @testset "Type stability" begin @@ -192,12 +199,12 @@ end ud = ForwardDiff.Dual.(u, seed) Dd = WF(ud, P, 0.0, xs, dxv) @test Dd isa ForwardDiff.Dual - @test ForwardDiff.value(Dd) ≈ D64 atol = 1e-12 - hfd = 1e-6 + @test ForwardDiff.value(Dd) ≈ D64 atol = 1.0e-12 + hfd = 1.0e-6 up = copy(u); up[3] += hfd um = copy(u); um[3] -= hfd pfd = (WF(up, P, 0.0, xs, dxv) - WF(um, P, 0.0, xs, dxv)) / (2hfd) - @test ForwardDiff.partials(Dd)[1] ≈ pfd rtol = 1e-5 + @test ForwardDiff.partials(Dd)[1] ≈ pfd rtol = 1.0e-5 # Symbolics.Num: symbolic build then numeric evaluation must equal the direct numeric call. @variables uu[1:5] @@ -205,7 +212,7 @@ end Dsym = WF(usym, P, 0.0, xs, dxv) @test Dsym isa Num gnum = build_function(Dsym, usym; expression = Val{false}) - @test gnum(u) ≈ D64 atol = 1e-12 + @test gnum(u) ≈ D64 atol = 1.0e-12 end @testset "SubArray view ingestion (production argument types)" begin @@ -225,6 +232,8 @@ end @test (@inferred WF(u, P, 0.0, xs_view, dxv_view)) isa Float64 bench_weno_sub(u, xs_view, dxv_view) # warmup - @test @allocated(bench_weno_sub(u, xs_view, dxv_view)) == 0 + let u_loc = u, xs_loc = xs_view, dxv_loc = dxv_view + @test @allocated(bench_weno_sub(u_loc, xs_loc, dxv_loc)) == 0 + end end end From 81717cb281aa5a65f716288179d3db7d2655ff6f Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Fri, 26 Jun 2026 16:42:27 +0300 Subject: [PATCH 4/9] test: introduce specialized function barrier to strictly measure zero-allocation on LTS --- test/Components/weno_nonuniform_core.jl | 28 +++++++++++++------------ 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/test/Components/weno_nonuniform_core.jl b/test/Components/weno_nonuniform_core.jl index c9bf99aa3..a9985d23d 100644 --- a/test/Components/weno_nonuniform_core.jl +++ b/test/Components/weno_nonuniform_core.jl @@ -24,6 +24,12 @@ bench_weno_f32(u::Vector{Float32}, x::AbstractVector{Float32}, dx::AbstractVecto bench_weno_sub(u::Vector{Float64}, x::SubArray{Float64, 1}, dx::SubArray{Float64, 1}) = MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, x, dx) +# Allocation must be measured from inside a fully specialized function barrier. On Julia 1.10 LTS, +# evaluating `@allocated f(...)` directly in (or under) the `@testset` scope boxes the Float64 result +# (16 bytes); `begin ...; nothing end` and `let` wrappers do not prevent this. Parametrizing on the +# callable `F` forces specialization so the call and its return value stay unboxed -> a true 0. +@inline measure_alloc(f::F, args::Vararg{Any, N}) where {F, N} = @allocated f(args...) + # Public-API helpers (vector / scalar dx). `dx` is unused by the kernel but required by the contract. wf(u, x) = WF(u, P, 0.0, x, diff(x)) wf_scalar(u, x) = WF(u, P, 0.0, x, 0.5) @@ -159,23 +165,20 @@ end u = sin.(xs) dxv = diff(xs) - # Warmup (compile) every form exercised under @allocated, including the direct invocation. + # Warmup (compile) every form measured below, including the direct invocation. bench_weno_f64(u, xs, dxv) bench_weno_f64(u, xs, 0.5) MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, xs, dxv) + measure_alloc(bench_weno_f64, u, xs, dxv) - # `let` blocks isolate the locals from the @testset scope so LTS does not box captured globals. - let u_loc = u, xs_loc = xs, dxv_loc = dxv - @test @allocated(bench_weno_f64(u_loc, xs_loc, dxv_loc)) == 0 - @test @allocated(bench_weno_f64(u_loc, xs_loc, 0.5)) == 0 - @test @allocated(MethodOfLines.weno_f_nonuniform(u_loc, (WENO_EPS_F64,), 0.0, xs_loc, dxv_loc)) == 0 - end + @test measure_alloc(bench_weno_f64, u, xs, dxv) == 0 + @test measure_alloc(bench_weno_f64, u, xs, 0.5) == 0 + @test measure_alloc(MethodOfLines.weno_f_nonuniform, u, (WENO_EPS_F64,), 0.0, xs, dxv) == 0 xs32 = Float32.(xs); u32 = Float32.(u); dxv32 = Float32.(dxv) bench_weno_f32(u32, xs32, dxv32) - let u_loc = u32, xs_loc = xs32, dxv_loc = dxv32 - @test @allocated(bench_weno_f32(u_loc, xs_loc, dxv_loc)) == 0 - end + measure_alloc(bench_weno_f32, u32, xs32, dxv32) + @test measure_alloc(bench_weno_f32, u32, xs32, dxv32) == 0 end @testset "Type stability" begin @@ -232,8 +235,7 @@ end @test (@inferred WF(u, P, 0.0, xs_view, dxv_view)) isa Float64 bench_weno_sub(u, xs_view, dxv_view) # warmup - let u_loc = u, xs_loc = xs_view, dxv_loc = dxv_view - @test @allocated(bench_weno_sub(u_loc, xs_loc, dxv_loc)) == 0 - end + measure_alloc(bench_weno_sub, u, xs_view, dxv_view) + @test measure_alloc(bench_weno_sub, u, xs_view, dxv_view) == 0 end end From c1ca4051357a4a37596226fd9c2b73aeba1e5e80 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Fri, 26 Jun 2026 22:47:18 +0300 Subject: [PATCH 5/9] docs: streamline comments and enforce strict literature citations in WENO core --- .../schemes/WENO/nonuniform_weno.jl | 70 ++++++----------- test/Components/weno_nonuniform_core.jl | 78 +++++++------------ 2 files changed, 54 insertions(+), 94 deletions(-) diff --git a/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index bc215b262..5cd291c42 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -1,18 +1,13 @@ -# Node-centered Lagrange WENO-5 on non-uniform grids. The three 3-point sub-stencil derivative -# estimates p_k'(x_i) are combined with closed-form derivative-ideal weights d_k under a -# Shi-Hu-Shu positive/negative split. References: Fornberg (1988); Jiang & Shu (1996); -# Shi, Hu & Shu (2002). +# Node-centered Lagrange WENO-5 first-derivative reconstruction on non-uniform grids. +# References: Fornberg (1988); Jiang & Shu (1996); Shi, Hu & Shu (2002). -# Fornberg (1988) finite-difference weights for a 3-point stencil, returning the m = 0, 1, 2 -# derivative operators. Unrolled over SVectors for stack allocation; index map ν = 1,2,3 -> 0,1,2. +# Fornberg (1988) finite-difference weights for a 3-point stencil; returns the m = 0, 1, 2 operators. @inline function _fornberg3_weights(α::NTuple{3, T}, xt::T) where {T} α0, α1, α2 = α - # n = 0 base case. n1m0 = one(T); n1m1 = zero(T); n1m2 = zero(T) c1 = one(T) - # n = 1: introduce α1. c2 = α1 - α0 r1 = c1 / c2 tA = α1 - xt @@ -28,7 +23,6 @@ n1m0 = a1m0; n1m1 = a1m1; n1m2 = a1m2 n2m0 = a2m0; n2m1 = a2m1; n2m2 = a2m2 - # n = 2: introduce α2. c2 = (α2 - α0) * (α2 - α1) r2 = c1 / c2 c3a = α2 - α0 @@ -45,9 +39,9 @@ b3m1 = r2 * (n2m0 - sB2 * n2m1) b3m2 = r2 * (2 * n2m1 - sB2 * n2m2) - m0 = SVector{3, T}(b1m0, b2m0, b3m0) # values - m1 = SVector{3, T}(b1m1, b2m1, b3m1) # first derivative - m2 = SVector{3, T}(b1m2, b2m2, b3m2) # second derivative + m0 = SVector{3, T}(b1m0, b2m0, b3m0) + m1 = SVector{3, T}(b1m1, b2m1, b3m1) + m2 = SVector{3, T}(b1m2, b2m2, b3m2) return m0, m1, m2 end @@ -56,8 +50,7 @@ end @inline _beta_nonneg(val::T) where {T <: AbstractFloat} = max(val, zero(val)) @inline _beta_nonneg(val) = IfElse.ifelse(val < zero(val), zero(val), val) -# Smoothness indicator β_k (Jiang & Shu 1996) via Simpson's rule over cell i, exact for the -# degree-2 sub-stencil integrand, together with the sub-stencil derivative estimate r_k = p_k'(x_i). +# Smoothness indicator β_k (Jiang & Shu 1996), Simpson quadrature over cell i. @inline function _substencil_beta_r( α::NTuple{3, T}, ua, ub, uc, xi, xL, xM, xph, Δx ) where {T} @@ -66,21 +59,20 @@ end _, m1M, m2M = _fornberg3_weights(α, xM) _, m1R, _ = _fornberg3_weights(α, xph) - r = _dot3(m1i, ua, ub, uc) # p_k'(x_i) - pL = _dot3(m1L, ua, ub, uc) # p_k'(x_{i-1/2}) - pM = _dot3(m1M, ua, ub, uc) # p_k'(x_M) - pR = _dot3(m1R, ua, ub, uc) # p_k'(x_{i+1/2}) - pp = _dot3(m2M, ua, ub, uc) # p_k'' (constant on a 3-point stencil) + r = _dot3(m1i, ua, ub, uc) + pL = _dot3(m1L, ua, ub, uc) + pM = _dot3(m1M, ua, ub, uc) + pR = _dot3(m1R, ua, ub, uc) + pp = _dot3(m2M, ua, ub, uc) - I1 = (Δx / 6) * (pL^2 + 4 * pM^2 + pR^2) # ∫ (p')^2 dx (Simpson, exact) - I2 = Δx * pp^2 # ∫ (p'')^2 dx + I1 = (Δx / 6) * (pL^2 + 4 * pM^2 + pR^2) + I2 = Δx * pp^2 val = Δx * I1 + Δx^3 * I2 β = _beta_nonneg(val) return β, r end -# Zero-allocation core. Geometry and weights are formed in Tx = eltype(x); promotion against -# eltype(u) (Symbolics.Num, ForwardDiff.Dual, Float32) occurs at the dot products. +# 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) Tx = eltype(x) θ = Tx(3) @@ -91,12 +83,10 @@ end u1 = u[1]; u2 = u[2]; u3 = u[3]; u4 = u[4]; u5 = u[5] end - # Reconstruction target is the center node x_i = x3. The cell [x_{i-1/2}, x_{i+1/2}] and its - # Simpson nodes are used only by the smoothness indicators. xi = x3 - xph = (x3 + x4) / 2 # x_{i+1/2} - xmh = (x2 + x3) / 2 # x_{i-1/2} - Δx = xph - xmh # cell width Δx_i + xph = (x3 + x4) / 2 + xmh = (x2 + x3) / 2 + Δx = xph - xmh xL = xmh xM = (xL + xph) / 2 @@ -108,16 +98,12 @@ 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: the unique d_k with Σ_k d_k p_k'(x_i) = P'_5(x_i), - # reducing to (1/6, 2/3, 1/6) on a uniform grid. Σ_k d_k = 1. + # 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)) d1 = one(Tx) - d0 - d2 - # Shi-Hu-Shu (2002) positive/negative split (θ = 3). In this node-centered direct-derivative - # topology the d_k remain non-negative, so the split is presently a defensive mechanism and a - # placeholder for future Taylor-series flux expansions that may introduce negative weights. The - # positive and negative branches are kept separate and recombined as σp·Rp − σm·Rm. + # Shi, Hu & Shu (2002) positive/negative weight splitting (θ = 3). dp0 = half * (d0 + θ * abs(d0)); dp1 = half * (d1 + θ * abs(d1)); dp2 = half * (d2 + θ * abs(d2)) dm0 = dp0 - d0; dm1 = dp1 - d1; dm2 = dp2 - d2 σp = dp0 + dp1 + dp2 @@ -140,16 +126,10 @@ end """ weno_f_nonuniform(u, p, t, x, dx) -Node-centered Lagrange WENO-5 reconstruction of a first spatial derivative on a non-uniform grid. - -Given the length-5 interior stencil values `u` and the corresponding absolute node coordinates -`x`, this implementation reconstructs the direct spatial derivative at the center node `x[3]` to -achieve 4th-order accuracy on non-uniform grids. It is a non-conservative formulation. The -smoothness parameter `ε = p[1]` regularizes the nonlinear weights; `t` and `dx` are accepted to -satisfy the `FunctionalScheme{5,0}` contract but are unused, as all geometry is taken from `x`. - -The coordinates `x` must be strictly monotonically increasing and distinct (`Δx_i > 0`): the -formulation divides by node differences and is undefined on degenerate grids. +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`). References: Fornberg (1988); Jiang & Shu (1996); Shi, Hu & Shu (2002). """ @@ -157,7 +137,7 @@ Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::Abst return _weno_f_nonuniform_core(u, p[1], x) end -# Scalar dx method: uniform-stepsize fallback required by the FunctionalScheme contract. +# 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) end diff --git a/test/Components/weno_nonuniform_core.jl b/test/Components/weno_nonuniform_core.jl index a9985d23d..accaa6bb1 100644 --- a/test/Components/weno_nonuniform_core.jl +++ b/test/Components/weno_nonuniform_core.jl @@ -1,10 +1,5 @@ -# Verification of the node-centered direct-derivative WENO-5 core (`weno_f_nonuniform`, -# FunctionalScheme{5,0}, returning du/dx at the center node x[3]). -# -# Properties under test: the nonlinear scheme is exact through degree 2 (ω_k = d_k there) while the -# linear ideal-weight decomposition is exact through degree 4; the derivative-ideal weights d_k form -# a convex partition (d_k >= 0, Σ d_k = 1) on every monotone grid, so the Shi-Hu-Shu negative branch -# is inert and is verified through convexity rather than a negative-weight path. +# Verification of the node-centered WENO-5 core: exact through degree 2 (nonlinear) and degree 4 +# (linear ideal-weight decomposition). using Test using MethodOfLines @@ -16,7 +11,7 @@ const P = [1.0e-6] const WENO_EPS_F64 = 1.0e-6 const WENO_EPS_F32 = 1.0f-6 -# Typed function barriers for @allocated tests (avoid closure boxing on LTS). +# Typed function barriers for @allocated tests. bench_weno_f64(u::Vector{Float64}, x::AbstractVector{Float64}, dx) = MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, x, dx) bench_weno_f32(u::Vector{Float32}, x::AbstractVector{Float32}, dx::AbstractVector{Float32}) = @@ -24,18 +19,14 @@ bench_weno_f32(u::Vector{Float32}, x::AbstractVector{Float32}, dx::AbstractVecto bench_weno_sub(u::Vector{Float64}, x::SubArray{Float64, 1}, dx::SubArray{Float64, 1}) = MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, x, dx) -# Allocation must be measured from inside a fully specialized function barrier. On Julia 1.10 LTS, -# evaluating `@allocated f(...)` directly in (or under) the `@testset` scope boxes the Float64 result -# (16 bytes); `begin ...; nothing end` and `let` wrappers do not prevent this. Parametrizing on the -# callable `F` forces specialization so the call and its return value stay unboxed -> a true 0. +# Allocation measured inside a specialized barrier; on 1.10 LTS the result is otherwise boxed. @inline measure_alloc(f::F, args::Vararg{Any, N}) where {F, N} = @allocated f(args...) -# Public-API helpers (vector / scalar dx). `dx` is unused by the kernel but required by the contract. +# Public-API helpers (vector / scalar dx). wf(u, x) = WF(u, P, 0.0, x, diff(x)) wf_scalar(u, x) = WF(u, P, 0.0, x, 0.5) -# Closed-form derivative-ideal weights (mirror of the core; the oracle in the linear-identity set is -# the analytic derivative, so a faulty formula fails against it). +# Closed-form derivative-ideal weights (mirror of the core). function ideal_weights(x) x1, x2, x3, x4, x5 = x d0 = ((x3 - x4) * (x3 - x5)) / ((x1 - x4) * (x1 - x5)) @@ -44,13 +35,13 @@ function ideal_weights(x) return (d0, d1, d2) end -# 3-point sub-stencil first derivative at xt via the shipped Fornberg engine (m = 1 weights). +# 3-point sub-stencil first derivative at xt (Fornberg m = 1 weights). function sub_deriv(α, ua, ub, uc, xt) w1 = MethodOfLines._fornberg3_weights(α, xt)[2] return w1[1] * ua + w1[2] * ub + w1[3] * uc end -# Linear (ideal-weight) reconstruction Σ d_k p_k'(x_i) assembled from the shipped pieces. +# Linear ideal-weight reconstruction Σ d_k p_k'(x_i). function linear_recon(x, u) d0, d1, d2 = ideal_weights(x) r0 = sub_deriv((x[1], x[2], x[3]), u[1], u[2], u[3], x[3]) @@ -62,7 +53,7 @@ end @testset "WENO Non-Uniform Core (Direct Derivative Reconstruction)" begin @testset "Polynomial exactness (degree <= 2)" begin - xs = [0.0, 0.6, 1.4, 2.1, 3.3] # ratios 0.6:0.8:0.7:1.2, center xs[3] = 1.4 + xs = [0.0, 0.6, 1.4, 2.1, 3.3] xc = xs[3] f0(x) = 1.7; df0(x) = 0.0 @@ -73,10 +64,10 @@ end @test wf(f1.(xs), xs) ≈ df1(xc) atol = 1.0e-13 @test wf(f2.(xs), xs) ≈ df2(xc) atol = 1.0e-12 - # Scalar-dx fallback hits the identical core: bit-for-bit equal. + # Scalar-dx fallback hits the identical core. @test wf_scalar(f2.(xs), xs) == wf(f2.(xs), xs) - # The nonlinear scheme is not exact for degree 3 at finite h (ω_k != d_k). + # Not exact for degree 3 at finite h. f3(x) = x^3 @test !isapprox(wf(f3.(xs), xs), 3xc^2; atol = 1.0e-6) end @@ -84,7 +75,7 @@ end @testset "Linear ideal-weight identity (degree <= 4)" begin for xs in ([0.0, 0.6, 1.4, 2.1, 3.3], [-0.3, 0.4, 0.55, 1.9, 2.4]) xc = xs[3] - # Convex partition: Σ d_k = 1 and d_k in [0, 1] (the Shi-Hu-Shu negative branch is inert). + # Convex partition: Σ d_k = 1, d_k in [0, 1]. d0, d1, d2 = ideal_weights(xs) @test d0 + d1 + d2 ≈ 1.0 atol = 1.0e-14 @test all(0 .<= (d0, d1, d2) .<= 1) @@ -103,57 +94,49 @@ end f(x) = sin(1.3x) + 0.5x df(x) = 1.3cos(1.3x) + 0.5 - # Frozen offsets: inner nodes perturbed ~10-20% off uniform, ends fixed. The relative - # geometry is constant under refinement (self-similar), so the asymptotic order is isolated. + # Self-similar refinement: fixed offsets isolate the asymptotic order. o = [-2.0, -1.13, 0.08, 0.91, 2.0] xc = 1.0 hs = [0.2, 0.1, 0.05, 0.025, 0.0125] errs = map(hs) do h xs = xc .+ h .* o - abs(wf(f.(xs), xs) - df(xs[3])) # evaluated at the actual center node xs[3] + abs(wf(f.(xs), xs) - df(xs[3])) end orders = [log2(errs[k] / errs[k + 1]) for k in 1:(length(hs) - 1)] - # Asymptotic order -> 4 (this FD scheme is 4th order on non-uniform grids). @test orders[end] > 3.85 @test orders[end - 1] > 3.7 - # Least-squares slope over the finest region (h <= 0.05) confirms ~4th order. lh = log.(hs[3:end]); le = log.(errs[3:end]); n = length(lh) slope = (n * sum(lh .* le) - sum(lh) * sum(le)) / (n * sum(lh .^ 2) - sum(lh)^2) @test slope > 3.8 - # A 1st/2nd/3rd-order regression would drag every order well below this. @test all(>(3.0), orders) end @testset "Extreme grid stretching (1:1e6)" begin grids = ( - [0.0, 1.0, 1.0 + 1.0e-6, 2.0 + 1.0e-6, 3.0 + 1.0e-6], # 1e-6 center cell - [0.0, 1.0e-6, 2.0e-6, 1.0, 2.0], # extreme clustering near the center + [0.0, 1.0, 1.0 + 1.0e-6, 2.0 + 1.0e-6, 3.0 + 1.0e-6], + [0.0, 1.0e-6, 2.0e-6, 1.0, 2.0], ) for g in grids - # Finiteness under severe ill-conditioning. @test isfinite(wf((x -> x^2).(g), g)) - # A linear field's derivative is exact regardless of the grid (Σ d_k = 1, every r_k - # equals the slope), so the Shi-Hu-Shu recombination must not corrupt it. + # Linear field: exact derivative on any grid. lin(x) = 2.0 + 3.0x @test wf(lin.(g), g) ≈ 3.0 rtol = 1.0e-6 - # Stability (not accuracy) target for a quadratic on a brutally stretched stencil. + # Stability (not accuracy) bound for a quadratic. @test abs(wf((x -> x^2).(g), g) - 2 * g[3]) < 1.0 - # d_k stay a convex partition even at 1:1e6, so d1 < 0 is unreachable. d0, d1, d2 = ideal_weights(g) @test all(isfinite, (d0, d1, d2)) @test d1 >= 0 @test d0 + d1 + d2 ≈ 1.0 atol = 1.0e-10 end - # A batch of deterministic extreme-ratio grids must never produce NaN/Inf. finite_all = true for k in 0:60 - s = 10.0^(k / 10 - 3) # cell-ratio scale spanning 1e-3 .. 1e3 + s = 10.0^(k / 10 - 3) g = cumsum([1.0, s, 1.0, s, 1.0]) .- 1.0 finite_all &= isfinite(wf((x -> sin(x)).(g), g)) end @@ -165,7 +148,7 @@ end u = sin.(xs) dxv = diff(xs) - # Warmup (compile) every form measured below, including the direct invocation. + # Warmup every measured form. bench_weno_f64(u, xs, dxv) bench_weno_f64(u, xs, 0.5) MethodOfLines.weno_f_nonuniform(u, (WENO_EPS_F64,), 0.0, xs, dxv) @@ -188,16 +171,14 @@ end u = f2.(xs) D64 = WF(u, P, 0.0, xs, dxv) - # Static inference: concrete Float64 return. @test (@inferred WF(u, P, 0.0, xs, dxv)) isa Float64 - # Float32 end-to-end. xs32 = Float32.(xs); dxv32 = Float32.(dxv); p32 = Float32[1.0f-6] D32 = WF(Float32.(u), p32, 0.0f0, xs32, dxv32) @test D32 isa Float32 @test D32 ≈ Float32(D64) rtol = 1.0f-4 - # ForwardDiff.Dual: differentiate the output w.r.t. u[3]; compare to a central difference. + # ForwardDiff.Dual w.r.t. u[3] vs central difference. seed = [0.0, 0.0, 1.0, 0.0, 0.0] ud = ForwardDiff.Dual.(u, seed) Dd = WF(ud, P, 0.0, xs, dxv) @@ -209,7 +190,7 @@ end pfd = (WF(up, P, 0.0, xs, dxv) - WF(um, P, 0.0, xs, dxv)) / (2hfd) @test ForwardDiff.partials(Dd)[1] ≈ pfd rtol = 1.0e-5 - # Symbolics.Num: symbolic build then numeric evaluation must equal the direct numeric call. + # Symbolics.Num: symbolic build then numeric evaluation. @variables uu[1:5] usym = collect(uu) Dsym = WF(usym, P, 0.0, xs, dxv) @@ -219,14 +200,13 @@ end end @testset "SubArray view ingestion (production argument types)" begin - # The discretizer passes x as a `@view s.grid[x][itap]` and the non-uniform dx as a view of - # the spacing vector. Verify that compiler views preserve the result, inference, and zero - # allocations, rather than only plain `Vector`s. - global_x = [0.0, 0.3, 0.9, 1.7, 2.2, 3.0, 3.4] # length-7 backing grid - xs_view = @view global_x[2:6] # length-5 SubArray (the interior stencil) - u = sin.(collect(xs_view)) # discvars indexing materializes a Vector + # SubArray views (as passed by the discretizer) must match Vector behavior in result, + # inference, and allocation. + 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] # length-4 SubArray of spacings + dxv_view = @view dx_full[2:5] Dref = WF(u, P, 0.0, collect(xs_view), diff(collect(xs_view))) Dview = WF(u, P, 0.0, xs_view, dxv_view) @@ -234,7 +214,7 @@ end @test Dview isa Float64 @test (@inferred WF(u, P, 0.0, xs_view, dxv_view)) isa Float64 - bench_weno_sub(u, xs_view, dxv_view) # warmup + bench_weno_sub(u, xs_view, dxv_view) measure_alloc(bench_weno_sub, u, xs_view, dxv_view) @test measure_alloc(bench_weno_sub, u, xs_view, dxv_view) == 0 end From 2cef2b3cb8925aaf96256d1468098594d873a2ea Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 03:16:22 +0300 Subject: [PATCH 6/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 e0aa0d364..4b37dfbfb 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -85,6 +85,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 3c745cad43458b110d76d9b1520ea9c01824f826 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 04:40:53 +0300 Subject: [PATCH 7/9] fix: update dispatch test --- src/discretization/schemes/WENO/WENO.jl | 2 -- test/Components/weno_dispatch.jl | 17 +++++++++++------ 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index 4b37dfbfb..e31510da5 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -1,5 +1,3 @@ -include("nonuniform_weno.jl") - """ Implements the WENO scheme of Jiang and Shu. Specified in https://repository.library.brown.edu/studio/item/bdr:297524/PDF/ (Page 8-9) diff --git a/test/Components/weno_dispatch.jl b/test/Components/weno_dispatch.jl index 488cc34a2..836127e8d 100644 --- a/test/Components/weno_dispatch.jl +++ b/test/Components/weno_dispatch.jl @@ -15,16 +15,21 @@ const weno_f_uniform = MethodOfLines.weno_f_uniform @test @allocated(weno_f(u, p, t, discx, dx_scalar)) == 0 end -@testset "WENO dispatch — non-uniform stub" begin +@testset "WENO dispatch — vector dx routes to non-uniform path" begin u = [1.0, 2.0, 3.0, 4.0, 5.0] p = [1.0e-6] t = 0.0 - x_vec = collect(0.0:0.1:0.4) + x_vec = [0.0, 0.1, 0.25, 0.45, 0.7] xv = @view x_vec[1:5] - dx_vec = diff(x_vec[1:5]) - @test_throws ArgumentError("WENO on non-uniform grids is not yet implemented.") weno_f( - u, p, t, xv, dx_vec - ) + dx_vec = diff(x_vec) + + m_vec = which(weno_f, (typeof(u), typeof(p), typeof(t), typeof(xv), typeof(dx_vec))) + m_scalar = which(weno_f, (typeof(u), typeof(p), typeof(t), typeof(xv), typeof(0.1))) + @test m_vec.sig.parameters[end] === AbstractVector + @test m_vec !== m_scalar + + @test weno_f(u, p, t, xv, dx_vec) == MethodOfLines.weno_f_nonuniform(u, p, t, xv, dx_vec) + @test weno_f(u, p, t, xv, dx_vec) isa Real end @testset "WENO dispatch — scalar fallback on vector grid" begin From 41e0af02d1c1ebc599940375751ced3a0341b096 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 05:02:26 +0300 Subject: [PATCH 8/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 2d777fc07f8021841989dd9d5d89afe4066bf086 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 28 Jun 2026 11:00:52 +0300 Subject: [PATCH 9/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