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 6237896690e55a12171fa8e0e143d8f64035bf00 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 +++ src/discretization/schemes/WENO/nonuniform_weno.jl | 8 ++++++++ 2 files changed, 11 insertions(+) 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 fbdc74576..b718554e4 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -196,3 +196,11 @@ 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 + +struct WENONonUniformBoundary{T} <: Function end + +Base.@propagate_inbounds @inline function (::WENONonUniformBoundary{T})( + u, p, t, x, dx::AbstractVector + ) where {T} + return weno_f_nonuniform(u, p, t, x, dx, Val(T)) +end From cea34861c686ecaed869172898543086ed38e462 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 9ce24aa332c79d897f96633e1621d1fb6f8310b5 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 a542d439f..692e256bf 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -100,7 +100,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 6ca86d68af2b99dc9410604393edfaf6e6e90ae2 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 692e256bf..a542d439f 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -100,6 +100,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 af29ac3e2e5efc74463629cc1beca2111694876c Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Tue, 30 Jun 2026 15:06:27 +0300 Subject: [PATCH 6/9] feat(WENO): integrate non-uniform grids --- src/MethodOfLines.jl | 2 +- src/discretization/schemes/WENO/WENO.jl | 12 +++- src/interface/scheme_types.jl | 2 + src/system_parsing/interior_map.jl | 4 +- test/Components/weno_boundary_integration.jl | 64 ++++++++++++++++++++ test/runtests.jl | 3 + 6 files changed, 81 insertions(+), 6 deletions(-) create mode 100644 test/Components/weno_boundary_integration.jl diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index a542d439f..db3b30695 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -99,8 +99,8 @@ include("discretization/schemes/upwind_difference/upwind_difference.jl") 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/WENO/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 e31510da5..b92fed721 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -78,8 +78,14 @@ end - `epsilon`: A quantity used to prevent vanishing denominators in the scheme, defaults to `1e-6`. More sensitive problems will benefit from a smaller value. It is defined as a functional scheme. """ function WENOScheme(; epsilon = 1.0e-6) - boundary_f = [nothing, nothing] - return FunctionalScheme{5, 0}( - weno_f, boundary_f, boundary_f, true, [epsilon], name = "WENO" + lower = WENONonUniformBoundary[WENONonUniformBoundary{1}(), WENONonUniformBoundary{2}()] + upper = WENONonUniformBoundary[WENONonUniformBoundary{5}(), WENONonUniformBoundary{4}()] + return FunctionalScheme{5, 5}( + weno_f, lower, upper, true, [epsilon], name = "WENO" ) end + +# extent: uniform dx -> legacy extrapolation buffer (2); vector dx -> Val boundaries (0). +const _WENOBoundaryVec = AbstractVector{<:WENONonUniformBoundary} +extent(::FunctionalScheme{<:Any, <:_WENOBoundaryVec}, dorder, dx::Number) = 2 +extent(::FunctionalScheme{<:Any, <:_WENOBoundaryVec}, dorder, dx::AbstractVector) = 0 diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 20615efd3..7d65f45d5 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -9,6 +9,8 @@ end extent(scheme::UpwindScheme, dorder) = 0 # dorder + scheme.order - 1 +extent(scheme::AbstractScheme, dorder, dx) = extent(scheme, dorder) + # Functional Schemes """ diff --git a/src/system_parsing/interior_map.jl b/src/system_parsing/interior_map.jl index a9d1cb00e..db6aecda3 100644 --- a/src/system_parsing/interior_map.jl +++ b/src/system_parsing/interior_map.jl @@ -118,10 +118,10 @@ function calculate_stencil_extents(s, u, discretization, orders, bcmap) for dorder in filter(isodd, orders[x]) ascheme = dorder == 1 ? advection_scheme : UpwindScheme() if !haslower - lowerextents[j] = max(lowerextents[j], extent(ascheme, dorder)) + lowerextents[j] = max(lowerextents[j], extent(ascheme, dorder, s.dxs[x])) end if !hasupper - upperextents[j] = max(upperextents[j], extent(ascheme, dorder)) + upperextents[j] = max(upperextents[j], extent(ascheme, dorder, s.dxs[x])) end end end diff --git a/test/Components/weno_boundary_integration.jl b/test/Components/weno_boundary_integration.jl new file mode 100644 index 000000000..4c0211013 --- /dev/null +++ b/test/Components/weno_boundary_integration.jl @@ -0,0 +1,64 @@ +# WENO discretizer integration: extent dispatch and interior map. + +using Test +using ModelingToolkit, MethodOfLines, DomainSets + +const M = MethodOfLines + +@parameters t x +@variables u(..) +Dt = Differential(t) +Dx = Differential(x) + +@testset "Grid-dispatched extent" begin + s = WENOScheme() + @test M.extent(s, 1, 1.0) == 2 # uniform dx + @test M.extent(s, 1, [0.1, 0.2, 0.3]) == 0 # vector dx + + # WENO specialization; not AbstractScheme fallback. + m_nonuniform = which(M.extent, (typeof(s), Int, Vector{Float64})) + m_fallback = which(M.extent, (M.AbstractScheme, Int, Any)) + @test m_nonuniform !== m_fallback + + # Grid-agnostic schemes ignore dx. + @test M.extent(UpwindScheme(), 1, [0.1, 0.2]) == M.extent(UpwindScheme(), 1) +end + +@testset "Interior map routes boundary nodes to Val WENO (non-uniform only)" begin + eq = Dt(u(t, x)) ~ -Dx(u(t, x)) + bcs = [u(0, x) ~ sinpi(x), u(t, 0.0) ~ 0.0, u(t, 1.0) ~ 0.0] + domains = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0)] + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + + function build_interiormap(disc) + v = M.VariableMap(pdesys, disc) + bcorders = Dict( + map(xx -> xx => M.d_orders(xx, M.get_bcs(pdesys)), M.PDEBase.all_ivs(v)) + ) + bmap = M.PDEBase.parse_bcs(M.get_bcs(pdesys), v, bcorders) + s = M.construct_discrete_space(v, disc) + eqs = M.get_eqs(pdesys) + eqs = eqs isa AbstractVector ? Vector{Equation}(eqs) : Equation[eqs] + return M.PDEBase.construct_var_equation_mapping(eqs, bmap, s, disc) + end + + nonuniform_grid = collect(range(0.0, 1.0, length = 21)) + nonuniform_grid[2:(end - 1)] .+= 0.004 .* sin.(1:19) + + im_nu = build_interiormap(MOLFiniteDifference([x => nonuniform_grid], t; advection_scheme = WENOScheme())) + im_u = build_interiormap(MOLFiniteDifference([x => 1 / 20], t; advection_scheme = WENOScheme())) + + pde_nu = first(keys(im_nu.I)) + pde_u = first(keys(im_u.I)) + + @test im_nu.stencil_extents[pde_nu] == ([0], [0]) + @test im_u.stencil_extents[pde_u] == ([2], [2]) + + interior_nu = vec(collect(im_nu.I[pde_nu])) + @test CartesianIndex(2) in interior_nu # Val(2) + @test CartesianIndex(20) in interior_nu # Val(4) + + interior_u = vec(collect(im_u.I[pde_u])) + @test !(CartesianIndex(2) in interior_u) + @test !(CartesianIndex(20) in interior_u) +end diff --git a/test/runtests.jl b/test/runtests.jl index 40f711166..033ea1a58 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -51,6 +51,9 @@ run_tests(; @safetestset "WENO Non-Uniform Boundary" begin include(joinpath(@__DIR__, "Components", "weno_nonuniform_boundary.jl")) end + @safetestset "WENO Boundary Integration" begin + include(joinpath(@__DIR__, "Components", "weno_boundary_integration.jl")) + end @safetestset "ODEFunction" begin include(joinpath(@__DIR__, "Components", "ODEFunction_test.jl")) end From 9e931d6dd545a4debdd6af48c0ed49d55d59eeff Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 5 Jul 2026 22:06:17 +0300 Subject: [PATCH 7/9] fix: address architectural edge cases and missing runtime test --- src/discretization/schemes/WENO/WENO.jl | 3 +- .../schemes/WENO/nonuniform_weno.jl | 6 ++ test/Components/weno_boundary_integration.jl | 59 ++++++++++++++++++- 3 files changed, 64 insertions(+), 4 deletions(-) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index b92fed721..c8a6a31d6 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -85,7 +85,8 @@ function WENOScheme(; epsilon = 1.0e-6) ) end -# extent: uniform dx -> legacy extrapolation buffer (2); vector dx -> Val boundaries (0). +# extent dispatch keys on typeof(lower): must stay Vector{<:WENONonUniformBoundary}. +# Replacing lower/upper with nothing restores 2-arg extent=0 and breaks uniform routing. const _WENOBoundaryVec = AbstractVector{<:WENONonUniformBoundary} extent(::FunctionalScheme{<:Any, <:_WENOBoundaryVec}, dorder, dx::Number) = 2 extent(::FunctionalScheme{<:Any, <:_WENOBoundaryVec}, dorder, dx::AbstractVector) = 0 diff --git a/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index b718554e4..e16ffdc5a 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -204,3 +204,9 @@ Base.@propagate_inbounds @inline function (::WENONonUniformBoundary{T})( ) where {T} return weno_f_nonuniform(u, p, t, x, dx, Val(T)) end + +Base.@propagate_inbounds @inline function (::WENONonUniformBoundary{T})( + u, p, t, x, dx::Number + ) where {T} + return weno_f_nonuniform(u, p, t, x, dx, Val(T)) +end diff --git a/test/Components/weno_boundary_integration.jl b/test/Components/weno_boundary_integration.jl index 4c0211013..f9a1a09fa 100644 --- a/test/Components/weno_boundary_integration.jl +++ b/test/Components/weno_boundary_integration.jl @@ -1,7 +1,7 @@ # WENO discretizer integration: extent dispatch and interior map. using Test -using ModelingToolkit, MethodOfLines, DomainSets +using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, SciMLBase const M = MethodOfLines @@ -10,6 +10,14 @@ const M = MethodOfLines Dt = Differential(t) Dx = Differential(x) +@testset "WENONonUniformBoundary scalar dx callable" begin + xs = collect(range(0.0, 1.0, length = 5)) + u = sin.(2π .* xs) + result = @inferred M.WENONonUniformBoundary{1}()(u, (1.0e-6,), 0.0, xs, 0.1) + @test result isa Float64 + @test isfinite(result) +end + @testset "Grid-dispatched extent" begin s = WENOScheme() @test M.extent(s, 1, 1.0) == 2 # uniform dx @@ -55,10 +63,55 @@ end @test im_u.stencil_extents[pde_u] == ([2], [2]) interior_nu = vec(collect(im_nu.I[pde_nu])) - @test CartesianIndex(2) in interior_nu # Val(2) - @test CartesianIndex(20) in interior_nu # Val(4) + @test !(CartesianIndex(1) in interior_nu) # wall: Dirichlet BC + @test CartesianIndex(2) in interior_nu # penultimate lower + @test CartesianIndex(20) in interior_nu # penultimate upper + @test !(CartesianIndex(21) in interior_nu) # wall: Dirichlet BC interior_u = vec(collect(im_u.I[pde_u])) + @test !(CartesianIndex(1) in interior_u) @test !(CartesianIndex(2) in interior_u) @test !(CartesianIndex(20) in interior_u) + @test !(CartesianIndex(21) in interior_u) +end + +@testset "Non-uniform WENO discretize + solve smoke (MMS)" begin + T_END = 0.2 + u_exact(x, t) = sin(2π * (x - t)) + + nonuniform_grid = collect(range(0.0, 1.0, length = 21)) + nonuniform_grid[2:(end - 1)] .+= 0.004 .* sin.(1:19) + x0, xL = nonuniform_grid[1], nonuniform_grid[end] + + eq = Dt(u(t, x)) ~ -Dx(u(t, x)) + bcs = [ + u(0.0, x) ~ sin(2π * x), + u(t, x0) ~ u_exact(x0, t), + u(t, xL) ~ u_exact(xL, t), + ] + domains = [t ∈ Interval(0.0, T_END), x ∈ Interval(x0, xL)] + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + + disc = MOLFiniteDifference([x => nonuniform_grid], t; advection_scheme = WENOScheme()) + prob = discretize(pdesys, disc) + sol = solve(prob, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat = [T_END]) + + @test SciMLBase.successful_retcode(sol) + + xs = sol[x] + u_num = sol[u(t, x)][end, :] + @test all(isfinite, u_num) + + u_ref = u_exact.(xs, T_END) + w = similar(xs, Float64) + w[1] = (xs[2] - xs[1]) / 2 + w[end] = (xs[end] - xs[end - 1]) / 2 + for i in 2:(length(xs) - 1) + w[i] = (xs[i + 1] - xs[i - 1]) / 2 + end + err = sqrt(sum(w .* abs2.(u_num .- u_ref))) + ref = max(sqrt(sum(w .* abs2.(u_ref))), eps()) + rel_l2 = err / ref + # measured rel_l2 ≈ 3.2e-3 at N=21, t=0.2 + @test rel_l2 < 0.01 end From ed0f18780f63a76fc6175d7c2cba9f80a8de13be Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Mon, 6 Jul 2026 00:06:59 +0300 Subject: [PATCH 8/9] test: implement genuine discretizer execution coverage and finalize PR polish --- src/discretization/schemes/WENO/WENO.jl | 5 +- .../schemes/WENO/nonuniform_weno.jl | 4 +- test/Components/weno_boundary_integration.jl | 119 ++++++++++++++---- 3 files changed, 104 insertions(+), 24 deletions(-) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index c8a6a31d6..9448d50e3 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -73,7 +73,10 @@ end end """ -`WENOScheme` of Jiang and Shu + WENOScheme(; epsilon = 1.0e-6) + +Jiang-Shu WENO-5 advection scheme for uniform and non-uniform grids. + ## Keyword Arguments - `epsilon`: A quantity used to prevent vanishing denominators in the scheme, defaults to `1e-6`. More sensitive problems will benefit from a smaller value. It is defined as a functional scheme. """ diff --git a/src/discretization/schemes/WENO/nonuniform_weno.jl b/src/discretization/schemes/WENO/nonuniform_weno.jl index e16ffdc5a..5d0af0d28 100644 --- a/src/discretization/schemes/WENO/nonuniform_weno.jl +++ b/src/discretization/schemes/WENO/nonuniform_weno.jl @@ -183,12 +183,12 @@ Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::Abst return _weno_f_nonuniform_core(u, p[1], x, Val(3)) end -# Scalar-dx method required by the FunctionalScheme{5,0} contract. +# Scalar-dx method required by the FunctionalScheme{5,5} contract. Base.@propagate_inbounds @inline function weno_f_nonuniform(u, p, t, x, dx::Number) return _weno_f_nonuniform_core(u, p[1], x, Val(3)) end -# 6-arg: explicit Val{Target}; dx unused (FunctionalScheme{5,0} contract). +# 6-arg: explicit Val{Target}; dx unused (FunctionalScheme{5,5} 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 diff --git a/test/Components/weno_boundary_integration.jl b/test/Components/weno_boundary_integration.jl index f9a1a09fa..73f558900 100644 --- a/test/Components/weno_boundary_integration.jl +++ b/test/Components/weno_boundary_integration.jl @@ -1,4 +1,4 @@ -# WENO discretizer integration: extent dispatch and interior map. +# WENO discretizer integration: extent dispatch, interior metadata, boundary execution. using Test using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, SciMLBase @@ -10,12 +10,34 @@ const M = MethodOfLines Dt = Differential(t) Dx = Differential(x) +function build_discrete_system(pdesys, disc) + v = M.VariableMap(pdesys, disc) + bcorders = Dict( + map(xx -> xx => M.d_orders(xx, M.get_bcs(pdesys)), M.PDEBase.all_ivs(v)) + ) + bmap = M.PDEBase.parse_bcs(M.get_bcs(pdesys), v, bcorders) + s = M.construct_discrete_space(v, disc) + eqs = M.get_eqs(pdesys) + eqs = eqs isa AbstractVector ? Vector{Equation}(eqs) : Equation[eqs] + im = M.PDEBase.construct_var_equation_mapping(eqs, bmap, s, disc) + return im, s, bmap +end + +function perturbed_nu_grid() + g = collect(range(0.0, 1.0, length = 21)) + g[2:(end - 1)] .+= 0.004 .* sin.(1:19) + return g +end + @testset "WENONonUniformBoundary scalar dx callable" begin xs = collect(range(0.0, 1.0, length = 5)) u = sin.(2π .* xs) result = @inferred M.WENONonUniformBoundary{1}()(u, (1.0e-6,), 0.0, xs, 0.1) @test result isa Float64 @test isfinite(result) + val5 = @inferred M.WENONonUniformBoundary{5}()(u, (1.0e-6,), 0.0, xs, 0.1) + @test val5 isa Float64 + @test isfinite(val5) end @testset "Grid-dispatched extent" begin @@ -23,6 +45,10 @@ end @test M.extent(s, 1, 1.0) == 2 # uniform dx @test M.extent(s, 1, [0.1, 0.2, 0.3]) == 0 # vector dx + # 2-arg extent counts nothing placeholders; legacy nothing-vector returned 2, Val callables return 0. + @test M.extent(WENOScheme(), 1) == 0 + @test M.extent(WENOScheme(), 1, 1.0) == 2 # 3-arg uniform path unchanged + # WENO specialization; not AbstractScheme fallback. m_nonuniform = which(M.extent, (typeof(s), Int, Vector{Float64})) m_fallback = which(M.extent, (M.AbstractScheme, Int, Any)) @@ -32,29 +58,20 @@ end @test M.extent(UpwindScheme(), 1, [0.1, 0.2]) == M.extent(UpwindScheme(), 1) end -@testset "Interior map routes boundary nodes to Val WENO (non-uniform only)" begin +@testset "Interior map extent dispatch (non-uniform vs uniform)" begin eq = Dt(u(t, x)) ~ -Dx(u(t, x)) bcs = [u(0, x) ~ sinpi(x), u(t, 0.0) ~ 0.0, u(t, 1.0) ~ 0.0] domains = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0)] @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) - function build_interiormap(disc) - v = M.VariableMap(pdesys, disc) - bcorders = Dict( - map(xx -> xx => M.d_orders(xx, M.get_bcs(pdesys)), M.PDEBase.all_ivs(v)) - ) - bmap = M.PDEBase.parse_bcs(M.get_bcs(pdesys), v, bcorders) - s = M.construct_discrete_space(v, disc) - eqs = M.get_eqs(pdesys) - eqs = eqs isa AbstractVector ? Vector{Equation}(eqs) : Equation[eqs] - return M.PDEBase.construct_var_equation_mapping(eqs, bmap, s, disc) - end - - nonuniform_grid = collect(range(0.0, 1.0, length = 21)) - nonuniform_grid[2:(end - 1)] .+= 0.004 .* sin.(1:19) + nonuniform_grid = perturbed_nu_grid() - im_nu = build_interiormap(MOLFiniteDifference([x => nonuniform_grid], t; advection_scheme = WENOScheme())) - im_u = build_interiormap(MOLFiniteDifference([x => 1 / 20], t; advection_scheme = WENOScheme())) + im_nu, _, _ = build_discrete_system( + pdesys, MOLFiniteDifference([x => nonuniform_grid], t; advection_scheme = WENOScheme()) + ) + im_u, _, _ = build_discrete_system( + pdesys, MOLFiniteDifference([x => 1 / 20], t; advection_scheme = WENOScheme()) + ) pde_nu = first(keys(im_nu.I)) pde_u = first(keys(im_u.I)) @@ -75,12 +92,72 @@ end @test !(CartesianIndex(21) in interior_u) end +@testset "Discretizer executes Val boundary callables at interior boundary nodes" begin + T_END = 0.1 + u_exact(x, t) = sin(2π * (x - t)) + + nonuniform_grid = perturbed_nu_grid() + x0, xL = nonuniform_grid[1], nonuniform_grid[end] + + eq = Dt(u(t, x)) ~ -Dx(u(t, x)) + bcs = [ + u(0.0, x) ~ sin(2π * x), + u(t, x0) ~ u_exact(x0, t), + Dx(u(t, xL)) ~ 0.0, + ] + domains = [t ∈ Interval(0.0, T_END), x ∈ Interval(x0, xL)] + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + + disc = MOLFiniteDifference([x => nonuniform_grid], t; advection_scheme = WENOScheme()) + im, s, bmap = build_discrete_system(pdesys, disc) + + pde = first(keys(im.I)) + interior_indices = vec(collect(im.I[pde])) + @test CartesianIndex(2) in interior_indices + @test CartesianIndex(20) in interior_indices + @test !(CartesianIndex(1) in interior_indices) + @test !(CartesianIndex(21) in interior_indices) + + F = WENOScheme() + u_field = im.var[pde] + j = M.x2i(s, u_field, x) + u_op = only(keys(bmap)) + bs = M.PDEBase.filter_interfaces(bmap[u_op][x]) + + f_lo, Itap_lo = M.get_f_and_taps(F, CartesianIndex(2), s, bs, (j, x), u_field) + f_hi, Itap_hi = M.get_f_and_taps(F, CartesianIndex(20), s, bs, (j, x), u_field) + @test f_lo === M.WENONonUniformBoundary{2}() + @test f_hi === M.WENONonUniformBoundary{4}() + + xs = collect(s.grid[x]) + uvals = sin.(2π .* xs) + dx_vec = collect(s.dxs[x]) + ps = (F.ps[1],) + itap_lo = map(I -> I[j], Itap_lo) + itap_hi = map(I -> I[j], Itap_hi) + discx_lo = @view xs[itap_lo] + discx_hi = @view xs[itap_hi] + dx_lo = @view dx_vec[itap_lo[1:(end - 1)]] + dx_hi = @view dx_vec[itap_hi[1:(end - 1)]] + u_disc_lo = uvals[itap_lo] + u_disc_hi = uvals[itap_hi] + + val_lo = @inferred f_lo(u_disc_lo, ps, 0.0, discx_lo, dx_lo) + val_hi = @inferred f_hi(u_disc_hi, ps, 0.0, discx_hi, dx_hi) + @test val_lo isa Float64 && isfinite(val_lo) + @test val_hi isa Float64 && isfinite(val_hi) + + prob = discretize(pdesys, disc) + sol = solve(prob, Tsit5(); abstol = 1.0e-8, reltol = 1.0e-8, saveat = [T_END]) + @test SciMLBase.successful_retcode(sol) + @test all(isfinite, sol[u(t, x)][end, :]) +end + @testset "Non-uniform WENO discretize + solve smoke (MMS)" begin T_END = 0.2 u_exact(x, t) = sin(2π * (x - t)) - nonuniform_grid = collect(range(0.0, 1.0, length = 21)) - nonuniform_grid[2:(end - 1)] .+= 0.004 .* sin.(1:19) + nonuniform_grid = perturbed_nu_grid() x0, xL = nonuniform_grid[1], nonuniform_grid[end] eq = Dt(u(t, x)) ~ -Dx(u(t, x)) @@ -94,7 +171,7 @@ end disc = MOLFiniteDifference([x => nonuniform_grid], t; advection_scheme = WENOScheme()) prob = discretize(pdesys, disc) - sol = solve(prob, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat = [T_END]) + sol = solve(prob, Tsit5(); abstol = 1.0e-8, reltol = 1.0e-8, saveat = [T_END]) @test SciMLBase.successful_retcode(sol) From 15c05b1dee98a253fedd9b6c4a7d6455c8c75d1b Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Mon, 6 Jul 2026 09:55:42 +0300 Subject: [PATCH 9/9] test: integration tests to extent dispatch and callable wiring --- test/Components/weno_boundary_integration.jl | 91 ++++---------------- 1 file changed, 17 insertions(+), 74 deletions(-) diff --git a/test/Components/weno_boundary_integration.jl b/test/Components/weno_boundary_integration.jl index 73f558900..42b9c359e 100644 --- a/test/Components/weno_boundary_integration.jl +++ b/test/Components/weno_boundary_integration.jl @@ -1,4 +1,4 @@ -# WENO discretizer integration: extent dispatch, interior metadata, boundary execution. +# WENO discretizer integration: extent dispatch, interior map, boundary callable wiring. using Test using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, SciMLBase @@ -29,32 +29,25 @@ function perturbed_nu_grid() return g end -@testset "WENONonUniformBoundary scalar dx callable" begin +@testset "WENONonUniformBoundary scalar dx" begin xs = collect(range(0.0, 1.0, length = 5)) u = sin.(2π .* xs) result = @inferred M.WENONonUniformBoundary{1}()(u, (1.0e-6,), 0.0, xs, 0.1) @test result isa Float64 @test isfinite(result) - val5 = @inferred M.WENONonUniformBoundary{5}()(u, (1.0e-6,), 0.0, xs, 0.1) - @test val5 isa Float64 - @test isfinite(val5) end @testset "Grid-dispatched extent" begin s = WENOScheme() - @test M.extent(s, 1, 1.0) == 2 # uniform dx - @test M.extent(s, 1, [0.1, 0.2, 0.3]) == 0 # vector dx - - # 2-arg extent counts nothing placeholders; legacy nothing-vector returned 2, Val callables return 0. + @test M.extent(s, 1, 1.0) == 2 + @test M.extent(s, 1, [0.1, 0.2, 0.3]) == 0 @test M.extent(WENOScheme(), 1) == 0 - @test M.extent(WENOScheme(), 1, 1.0) == 2 # 3-arg uniform path unchanged + @test M.extent(WENOScheme(), 1, 1.0) == 2 - # WENO specialization; not AbstractScheme fallback. m_nonuniform = which(M.extent, (typeof(s), Int, Vector{Float64})) m_fallback = which(M.extent, (M.AbstractScheme, Int, Any)) @test m_nonuniform !== m_fallback - # Grid-agnostic schemes ignore dx. @test M.extent(UpwindScheme(), 1, [0.1, 0.2]) == M.extent(UpwindScheme(), 1) end @@ -80,10 +73,10 @@ end @test im_u.stencil_extents[pde_u] == ([2], [2]) interior_nu = vec(collect(im_nu.I[pde_nu])) - @test !(CartesianIndex(1) in interior_nu) # wall: Dirichlet BC - @test CartesianIndex(2) in interior_nu # penultimate lower - @test CartesianIndex(20) in interior_nu # penultimate upper - @test !(CartesianIndex(21) in interior_nu) # wall: Dirichlet BC + @test !(CartesianIndex(1) in interior_nu) + @test CartesianIndex(2) in interior_nu + @test CartesianIndex(20) in interior_nu + @test !(CartesianIndex(21) in interior_nu) interior_u = vec(collect(im_u.I[pde_u])) @test !(CartesianIndex(1) in interior_u) @@ -92,7 +85,7 @@ end @test !(CartesianIndex(21) in interior_u) end -@testset "Discretizer executes Val boundary callables at interior boundary nodes" begin +@testset "Discretizer boundary callable routing and execution" begin T_END = 0.1 u_exact(x, t) = sin(2π * (x - t)) @@ -112,17 +105,10 @@ end im, s, bmap = build_discrete_system(pdesys, disc) pde = first(keys(im.I)) - interior_indices = vec(collect(im.I[pde])) - @test CartesianIndex(2) in interior_indices - @test CartesianIndex(20) in interior_indices - @test !(CartesianIndex(1) in interior_indices) - @test !(CartesianIndex(21) in interior_indices) - F = WENOScheme() u_field = im.var[pde] j = M.x2i(s, u_field, x) - u_op = only(keys(bmap)) - bs = M.PDEBase.filter_interfaces(bmap[u_op][x]) + bs = M.PDEBase.filter_interfaces(bmap[only(keys(bmap))][x]) f_lo, Itap_lo = M.get_f_and_taps(F, CartesianIndex(2), s, bs, (j, x), u_field) f_hi, Itap_hi = M.get_f_and_taps(F, CartesianIndex(20), s, bs, (j, x), u_field) @@ -135,15 +121,12 @@ end ps = (F.ps[1],) itap_lo = map(I -> I[j], Itap_lo) itap_hi = map(I -> I[j], Itap_hi) - discx_lo = @view xs[itap_lo] - discx_hi = @view xs[itap_hi] - dx_lo = @view dx_vec[itap_lo[1:(end - 1)]] - dx_hi = @view dx_vec[itap_hi[1:(end - 1)]] - u_disc_lo = uvals[itap_lo] - u_disc_hi = uvals[itap_hi] - - val_lo = @inferred f_lo(u_disc_lo, ps, 0.0, discx_lo, dx_lo) - val_hi = @inferred f_hi(u_disc_hi, ps, 0.0, discx_hi, dx_hi) + val_lo = @inferred f_lo( + uvals[itap_lo], ps, 0.0, @view(xs[itap_lo]), @view(dx_vec[itap_lo[1:(end - 1)]]) + ) + val_hi = @inferred f_hi( + uvals[itap_hi], ps, 0.0, @view(xs[itap_hi]), @view(dx_vec[itap_hi[1:(end - 1)]]) + ) @test val_lo isa Float64 && isfinite(val_lo) @test val_hi isa Float64 && isfinite(val_hi) @@ -152,43 +135,3 @@ end @test SciMLBase.successful_retcode(sol) @test all(isfinite, sol[u(t, x)][end, :]) end - -@testset "Non-uniform WENO discretize + solve smoke (MMS)" begin - T_END = 0.2 - u_exact(x, t) = sin(2π * (x - t)) - - nonuniform_grid = perturbed_nu_grid() - x0, xL = nonuniform_grid[1], nonuniform_grid[end] - - eq = Dt(u(t, x)) ~ -Dx(u(t, x)) - bcs = [ - u(0.0, x) ~ sin(2π * x), - u(t, x0) ~ u_exact(x0, t), - u(t, xL) ~ u_exact(xL, t), - ] - domains = [t ∈ Interval(0.0, T_END), x ∈ Interval(x0, xL)] - @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) - - disc = MOLFiniteDifference([x => nonuniform_grid], t; advection_scheme = WENOScheme()) - prob = discretize(pdesys, disc) - sol = solve(prob, Tsit5(); abstol = 1.0e-8, reltol = 1.0e-8, saveat = [T_END]) - - @test SciMLBase.successful_retcode(sol) - - xs = sol[x] - u_num = sol[u(t, x)][end, :] - @test all(isfinite, u_num) - - u_ref = u_exact.(xs, T_END) - w = similar(xs, Float64) - w[1] = (xs[2] - xs[1]) / 2 - w[end] = (xs[end] - xs[end - 1]) / 2 - for i in 2:(length(xs) - 1) - w[i] = (xs[i + 1] - xs[i - 1]) / 2 - end - err = sqrt(sum(w .* abs2.(u_num .- u_ref))) - ref = max(sqrt(sum(w .* abs2.(u_ref))), eps()) - rel_l2 = err / ref - # measured rel_l2 ≈ 3.2e-3 at N=21, t=0.2 - @test rel_l2 < 0.01 -end