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..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/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index 461d4dcd2..5cd291c42 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -1,6 +1,143 @@ -@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 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; returns the m = 0, 1, 2 operators. +@inline function _fornberg3_weights(α::NTuple{3, T}, xt::T) where {T} + α0, α1, α2 = α + + n1m0 = one(T); n1m1 = zero(T); n1m2 = zero(T) + c1 = one(T) + + 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 + + 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) + m1 = SVector{3, T}(b1m1, b2m1, b3m1) + m2 = SVector{3, T}(b1m2, b2m2, b3m2) + return m0, m1, m2 end -Base.@propagate_inbounds @inline weno_f_nonuniform(u, p, t, x, dx::Real) = - weno_f_uniform(u, p, t, x, dx) +@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), Simpson quadrature over cell 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) + 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) + I2 = Δx * pp^2 + val = Δx * I1 + Δx^3 * I2 + β = _beta_nonneg(val) + return β, r +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) + Tx = eltype(x) + θ = 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 + + xi = x3 + xph = (x3 + x4) / 2 + xmh = (x2 + x3) / 2 + Δx = xph - xmh + 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; 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 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 + σ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 + +""" + weno_f_nonuniform(u, p, t, x, dx) + +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). +""" +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 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_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 diff --git a/test/Components/weno_nonuniform_core.jl b/test/Components/weno_nonuniform_core.jl new file mode 100644 index 000000000..accaa6bb1 --- /dev/null +++ b/test/Components/weno_nonuniform_core.jl @@ -0,0 +1,221 @@ +# Verification of the node-centered WENO-5 core: exact through degree 2 (nonlinear) and degree 4 +# (linear ideal-weight decomposition). + +using Test +using MethodOfLines +using Symbolics +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. +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) + +# 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). +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). +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 (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). +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] + 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 = 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. + @test wf_scalar(f2.(xs), xs) == wf(f2.(xs), xs) + + # 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 + + @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, 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) + + 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 = 1.0e-11 atol = 1.0e-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 + + # 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])) + end + orders = [log2(errs[k] / errs[k + 1]) for k in 1:(length(hs) - 1)] + + @test orders[end] > 3.85 + @test orders[end - 1] > 3.7 + 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 + @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], + [0.0, 1.0e-6, 2.0e-6, 1.0, 2.0], + ) + for g in grids + @test isfinite(wf((x -> x^2).(g), g)) + + # 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) bound for a quadratic. + @test abs(wf((x -> x^2).(g), g) - 2 * g[3]) < 1.0 + + 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 + + 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 + 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) + + # 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) + measure_alloc(bench_weno_f64, u, xs, dxv) + + @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) + measure_alloc(bench_weno_f32, u32, xs32, dxv32) + @test measure_alloc(bench_weno_f32, 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) + + @test (@inferred WF(u, P, 0.0, xs, dxv)) isa Float64 + + 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 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) + @test Dd isa ForwardDiff.Dual + @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 = 1.0e-5 + + # Symbolics.Num: symbolic build then numeric evaluation. + @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 = 1.0e-12 + end + + @testset "SubArray view ingestion (production argument types)" begin + # 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] + + 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 + + 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 +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