From 6d6f9c34e05dddbe683d967536efb86056414ee3 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Sun, 14 Jun 2026 12:45:08 +0300 Subject: [PATCH 1/2] feat(bc): implement graceful well-posedness check for missing boundary conditions --- src/MethodOfLines.jl | 1 + src/discretization/boundary_validation.jl | 304 ++++++++++++++++++++ src/system_parsing/interior_map.jl | 1 + test/Components/boundary_validation_test.jl | 259 +++++++++++++++++ test/runtests.jl | 5 +- 5 files changed, 569 insertions(+), 1 deletion(-) create mode 100644 src/discretization/boundary_validation.jl create mode 100644 test/Components/boundary_validation_test.jl diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 0eb00db1f..ce2888e04 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -105,6 +105,7 @@ include("discretization/schemes/WENO/WENO.jl") include("discretization/schemes/integral_expansion/integral_expansion.jl") # System Discretization +include("discretization/boundary_validation.jl") include("discretization/generate_finite_difference_rules.jl") include("discretization/generate_bc_eqs.jl") include("discretization/generate_ic_defaults.jl") diff --git a/src/discretization/boundary_validation.jl b/src/discretization/boundary_validation.jl new file mode 100644 index 000000000..31fc925b6 --- /dev/null +++ b/src/discretization/boundary_validation.jl @@ -0,0 +1,304 @@ +using ModelingToolkit: operation, arguments, iscall, Differential +using Symbolics: unwrap + +""" + validate_system_wellposedness(pdes, bmap, s, _disc) + +Validate the well-posedness of the PDE system before discretization. For every +dependent variable and spatial dimension, the highest spatial derivative order +sets the required number of boundary conditions. + +For order ≥ 2, boundary conditions are counted by effective constraint weight: +each single-point equation contributes 1; a coupled or periodic equation that +evaluates `u` at two endpoints contributes 2. Several independent conditions +may share the same spatial point (e.g. clamped ends `u = 0` and `u_x = 0` at +`x = 0`). + +A first-order spatial derivative whose coefficient is not a compile-time numeric +constant (e.g. Burgers' `-u * u_x`, parameterized `v * u_x`, or spatial `sin(x) * u_x`) +is treated as needing data at *both* ends of the domain, since the upwind direction +cannot be fixed before the solve. +""" +function validate_system_wellposedness(pdes, bmap, s, _disc) + for u_call in s.ū + u_op = operation(unwrap(u_call)) + spatial_ivs = s.x̄ + + for x in spatial_ivs + max_order = 0 + has_dynamic_advection = false + + for eq in pdes + max_order = max(max_order, + get_max_derivative_order(eq.lhs, u_op, x), + get_max_derivative_order(eq.rhs, u_op, x)) + end + + if max_order == 1 + for eq in pdes + if is_dynamic_advection(eq.lhs, u_op, x) || + is_dynamic_advection(eq.rhs, u_op, x) + has_dynamic_advection = true + break + end + end + end + + if max_order > 0 + u_bcs = get(bmap, u_op, nothing) + u_bcs_x = u_bcs !== nothing ? get(u_bcs, x, nothing) : nothing + provided_bcs = u_bcs_x !== nothing ? length(u_bcs_x) : 0 + + if max_order >= 2 + effective_bcs = count_effective_boundary_conditions( + u_bcs_x, u_call, u_op, x) + if effective_bcs < max_order + throw(ArgumentError( + "Ill-posed PDE: $(u_op) in $(x) has order-$(max_order) " * + "spatial derivatives but only $(effective_bcs) effective " * + "boundary constraint(s) provided; need $(max_order)." + )) + end + elseif max_order == 1 && has_dynamic_advection + provided_locations = count_boundary_locations( + u_bcs_x, u_call, u_op, x) + if provided_locations < 2 + throw(ArgumentError( + "Ill-posed PDE: $(u_op) in $(x) has a non-constant advection " * + "coefficient, so the upwind direction is not fixed and boundary " * + "data is required at BOTH ends; only $(provided_locations) provided." + )) + end + elseif max_order == 1 && !has_dynamic_advection && provided_bcs < 1 + throw(ArgumentError( + "Ill-posed PDE: $(u_op) in $(x) has a 1st-order spatial " * + "derivative but no boundary condition was provided." + )) + end + end + end + end +end + +""" + _get_equation_locations(u_bcs_x, u_call, u_op, x) + +Return a vector of pairs mapping each boundary condition equation to a `Set` of +its unique spatial locations for the dimension `x`. Internal helper function. +""" +function _get_equation_locations(u_bcs_x, u_call, u_op, x) + u_bcs_x === nothing && return () + x_index = findfirst(isequal(x), arguments(unwrap(u_call))) + x_index === nothing && return () + + eq_locs = Pair{Any, Set{Any}}[] + for bc_obj in u_bcs_x + bc_eq = boundary_condition_equation(bc_obj) + unique_bounds = Set() + if hasproperty(bc_eq, :lhs) && hasproperty(bc_eq, :rhs) + extract_boundary_locations!(unique_bounds, bc_eq.lhs, u_op, x_index, x) + extract_boundary_locations!(unique_bounds, bc_eq.rhs, u_op, x_index, x) + end + push!(eq_locs, bc_eq => unique_bounds) + end + return eq_locs +end + +""" + count_boundary_locations(u_bcs_x, u_call, u_op, x) + +Return the number of unique spatial locations at which `u_op` appears in the +boundary conditions for dimension `x`. Numeric locations are normalised to +`Float64` so that e.g. `0` and `0.0` hash to the same key. +""" +function count_boundary_locations(u_bcs_x, u_call, u_op, x) + eq_locs = _get_equation_locations(u_bcs_x, u_call, u_op, x) + eq_locs === () && return 0 + + unique_bounds = Set() + for (_, locs) in eq_locs + union!(unique_bounds, locs) + end + return length(unique_bounds) +end + +""" + boundary_condition_equation(bc_obj) + +Return the underlying equation for a boundary condition object. +""" +boundary_condition_equation(bc_obj) = + hasproperty(bc_obj, :eq) ? bc_obj.eq : bc_obj + +""" + count_effective_boundary_conditions(u_bcs_x, u_call, u_op, x) + +Return the effective number of boundary constraints for dimension `x`. Each +single-point equation contributes 1; a coupled or periodic equation that +evaluates `u_op` at two or more unique spatial locations contributes 2. + +Periodic conditions appear twice in the parsed boundary map (once per endpoint); +they are deduplicated by structural equality of the underlying equation before +counting. Every counted equation contributes at least 1 constraint. +""" +function count_effective_boundary_conditions(u_bcs_x, u_call, u_op, x) + eq_locs = _get_equation_locations(u_bcs_x, u_call, u_op, x) + eq_locs === () && return 0 + + function structurally_equal(a, b) + if hasproperty(a, :lhs) && hasproperty(a, :rhs) && + hasproperty(b, :lhs) && hasproperty(b, :rhs) + return isequal(a.lhs, b.lhs) && isequal(a.rhs, b.rhs) + end + return isequal(a, b) + end + + seen_eqs = Any[] + total = 0 + for (bc_eq, unique_bounds) in eq_locs + if any(seen -> structurally_equal(seen, bc_eq), seen_eqs) + continue + end + push!(seen_eqs, bc_eq) + total += length(unique_bounds) >= 2 ? 2 : 1 + end + return total +end + +""" + extract_boundary_locations!(unique_bounds, expr, u_op, x_index, x) + +Walk the AST of a boundary condition and push the spatial location of every +`u_op(..., x_loc, ...)` occurrence into `unique_bounds`. +""" +function extract_boundary_locations!(unique_bounds::Set, expr, u_op, x_index, x) + expr = unwrap(expr) + if iscall(expr) + op = operation(expr) + + if isequal(op, u_op) + args = arguments(expr) + if x_index <= length(args) + raw_val = unwrap(args[x_index]) + normalized_val = raw_val isa Number ? Float64(raw_val) : raw_val + push!(unique_bounds, normalized_val) + else + throw(ArgumentError( + "Malformed boundary condition: $(u_op) needs at least $(x_index) " * + "arguments for dimension $(x), got $(length(args))." + )) + end + return + end + + for arg in arguments(expr) + extract_boundary_locations!(unique_bounds, arg, u_op, x_index, x) + end + end + return +end + +""" + get_max_derivative_order(expr, u_op, x) + +Return the highest order of a spatial derivative of `u_op` with respect to `x` +appearing anywhere in `expr`, or `-1` if `u_op` does not appear at all. +""" +function get_max_derivative_order(expr, u_op, x) + expr = unwrap(expr) + if iscall(expr) + op = operation(expr) + + if isequal(op, u_op) + return 0 + end + + inner_max = -1 + for arg in arguments(expr) + m = get_max_derivative_order(arg, u_op, x) + if m > inner_max + inner_max = m + end + end + + if op isa Differential && isequal(op.x, x) && inner_max >= 0 + return op.order + inner_max + end + + return inner_max + end + return isequal(expr, u_op) ? 0 : -1 +end + +""" + is_static_numeric_constant(expr) + +Return `true` when `expr` is a numeric literal, or a product/quotient of numeric +literals only. +""" +function is_static_numeric_constant(expr) + expr = unwrap(expr) + if expr isa Number + return true + end + val = Symbolics.value(expr) + if val isa Number + return true + end + if iscall(expr) + op = operation(expr) + if isequal(op, *) || isequal(op, /) + return all(is_static_numeric_constant, arguments(expr)) + end + end + return false +end + +""" + is_dynamic_advection(expr, u_op, x) + +Return `true` when `expr` contains a 1st-order spatial derivative of `u_op` +multiplied (or divided) by a coefficient that is not a compile-time numeric +constant. Parameters, states, and independent variables all make the upwind +direction unknown at discretization time. +""" +function is_dynamic_advection(expr, u_op, x) + expr = unwrap(expr) + if iscall(expr) + op = operation(expr) + + if isequal(op, *) || isequal(op, /) + args = arguments(expr) + has_deriv = false + for arg in args + arg_unwrapped = unwrap(arg) + arg_op = iscall(arg_unwrapped) ? operation(arg_unwrapped) : nothing + if arg_op isa Differential && isequal(arg_op.x, x) + if get_max_derivative_order(arg_unwrapped, u_op, arg_op.x) == 1 + has_deriv = true + break + end + end + end + + if has_deriv + for arg in args + arg_unwrapped = unwrap(arg) + arg_op = iscall(arg_unwrapped) ? operation(arg_unwrapped) : nothing + is_deriv = arg_op isa Differential && isequal(arg_op.x, x) && + get_max_derivative_order(arg_unwrapped, u_op, arg_op.x) == 1 + if !is_deriv && !is_static_numeric_constant(arg_unwrapped) + return true + end + end + end + end + + for arg in arguments(expr) + if is_dynamic_advection(arg, u_op, x) + return true + end + end + end + return false +end diff --git a/src/system_parsing/interior_map.jl b/src/system_parsing/interior_map.jl index a9d1cb00e..92e54494e 100644 --- a/src/system_parsing/interior_map.jl +++ b/src/system_parsing/interior_map.jl @@ -31,6 +31,7 @@ function PDEBase.construct_var_equation_mapping( pdes::Vector{Equation}, boundarymap, s::DiscreteSpace{N, M}, discretization::MOLFiniteDifference ) where {N, M} + validate_system_wellposedness(pdes, boundarymap, s, discretization) @assert length(pdes) == M "There must be the same number of equations and unknowns, got $(length(pdes)) equations and $(M) unknowns" m = buildmatrix(pdes, s) varmap = Dict(build_variable_mapping(m, s.ū, pdes)) diff --git a/test/Components/boundary_validation_test.jl b/test/Components/boundary_validation_test.jl new file mode 100644 index 000000000..057dab23c --- /dev/null +++ b/test/Components/boundary_validation_test.jl @@ -0,0 +1,259 @@ +using MethodOfLines, ModelingToolkit, DomainSets, Test, PDEBase +using ModelingToolkit: get_bcs, get_eqs +using PDEBase: VariableMap, cardinalize_eqs!, parse_bcs, construct_discrete_space, + d_orders, all_ivs + +@parameters t x y x1 x2 vel +@variables u(..) v(..) + +function run_boundary_validation(pdesys, disc) + vmap = VariableMap(pdesys, disc) + cardinalize_eqs!(pdesys) + bcorders = Dict(iv => d_orders(iv, get_bcs(pdesys)) for iv in all_ivs(vmap)) + boundarymap = parse_bcs(get_bcs(pdesys), vmap, bcorders) + s = construct_discrete_space(vmap, disc) + MethodOfLines.validate_system_wellposedness(get_eqs(pdesys), boundarymap, s, disc) + return nothing +end + +function standard_1d_disc(; + xspan = (0.0, 1.0), + tspan = (0.0, 1.0), + dx = 0.1, + advection = nothing, +) + t0, t1 = tspan + x0, x1 = xspan + disc_kwargs = advection === nothing ? () : (; advection_scheme = advection) + domains = [t ∈ Interval(t0, t1), x ∈ Interval(x0, x1)] + disc = MOLFiniteDifference([x => dx], t; disc_kwargs...) + return domains, disc +end + +function build_1d_pdesys(eq, bcs; xspan = (0.0, 1.0), tspan = (0.0, 1.0)) + domains, _ = standard_1d_disc(; xspan, tspan) + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + return pdesys +end + +function build_1d_pdesys_multi(eqs, bcs, depvars; xspan = (0.0, 1.0), tspan = (0.0, 1.0)) + domains, _ = standard_1d_disc(; xspan, tspan) + @named pdesys = PDESystem(eqs, bcs, domains, [t, x], depvars) + return pdesys +end + +function standard_2d_disc(; + xspan = (0.0, 1.0), + yspan = (0.0, 1.0), + tspan = (0.0, 1.0), + dx = 0.1, + dy = 0.1, +) + t0, t1 = tspan + x0, x1b = xspan + y0, y1b = yspan + domains = [ + t ∈ Interval(t0, t1), + x ∈ Interval(x0, x1b), + y ∈ Interval(y0, y1b), + ] + disc = MOLFiniteDifference([x => dx, y => dy], t) + return domains, disc +end + +function build_2d_pdesys(eq, bcs; xspan = (0.0, 1.0), yspan = (0.0, 1.0), tspan = (0.0, 1.0)) + domains, _ = standard_2d_disc(; xspan, yspan, tspan) + @named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)]) + return pdesys +end + +@testset "Boundary condition well-posedness" begin + + @testset "1st-order static advection" begin + Dt, Dx = Differential(t), Differential(x) + eq = Dt(u(t, x)) ~ -Dx(u(t, x)) + _, disc = standard_1d_disc(; advection = UpwindScheme()) + + @test run_boundary_validation( + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) === nothing + + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(eq, [u(0, x) ~ 0.0]), disc) + end + + @testset "1st-order dynamic advection" begin + Dt, Dx = Differential(t), Differential(x) + _, disc = standard_1d_disc(; advection = UpwindScheme()) + + eq_burgers = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x)) + @test run_boundary_validation( + build_1d_pdesys(eq_burgers, [u(0, x) ~ 1.0, u(t, 0) ~ 1.0, u(t, 1) ~ 1.0]), + disc, + ) === nothing + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(eq_burgers, [u(0, x) ~ 1.0, u(t, 0) ~ 1.0]), disc) + + eq_param = Dt(u(t, x)) ~ -vel * Dx(u(t, x)) + @test run_boundary_validation( + build_1d_pdesys(eq_param, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) === nothing + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(eq_param, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) + end + + @testset "2nd-order spatial derivatives" begin + Dt, Dxx = Differential(t), Differential(x)^2 + eq = Dt(u(t, x)) ~ Dxx(u(t, x)) + _, disc = standard_1d_disc() + + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) + @test run_boundary_validation( + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) === nothing + end + + @testset "3rd-order spatial derivatives" begin + Dt, Dx, Dx3 = Differential(t), Differential(x), Differential(x)^3 + eq = Dt(u(t, x)) ~ -Dx3(u(t, x)) + _, disc = standard_1d_disc() + + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) + @test run_boundary_validation( + build_1d_pdesys(eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ 0.0, + u(t, 1) ~ 0.0, + Dx(u(t, 0)) ~ 0.0, + ]), disc) === nothing + end + + @testset "4th-order spatial derivatives" begin + Dt, Dx, Dxxxx = Differential(t), Differential(x), Differential(x)^4 + eq = Dt(u(t, x)) ~ -Dxxxx(u(t, x)) + _, disc = standard_1d_disc() + + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ 0.0, + u(t, 1) ~ 0.0, + Dx(u(t, 0)) ~ 0.0, + ]), disc) + @test run_boundary_validation( + build_1d_pdesys(eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ 0.0, + Dx(u(t, 0)) ~ 0.0, + u(t, 1) ~ 0.0, + Dx(u(t, 1)) ~ 0.0, + ]), disc) === nothing + end + + @testset "Periodic boundaries" begin + Dt, Dx, Dxx, Dxxxx = Differential(t), Differential(x), Differential(x)^2, + Differential(x)^4 + _, disc = standard_1d_disc() + + periodic_bcs = [u(0, x) ~ 0.0, u(t, 0) ~ u(t, 1)] + @test run_boundary_validation( + build_1d_pdesys(Dt(u(t, x)) ~ Dxx(u(t, x)), periodic_bcs), disc) === nothing + + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(Dt(u(t, x)) ~ -Dxxxx(u(t, x)), periodic_bcs), disc) + + @test run_boundary_validation( + build_1d_pdesys(Dt(u(t, x)) ~ -Dxxxx(u(t, x)), [ + u(0, x) ~ 0.0, + u(t, 0) ~ u(t, 1), + Dx(u(t, 0)) ~ Dx(u(t, 1)), + ]), disc) === nothing + end + + @testset "Coupled systems" begin + _, disc = standard_1d_disc() + + Dt_u, Dxx = Differential(t), Differential(x)^2 + eqs = [Dt_u(u(t, x)) ~ Dxx(u(t, x)), Dt_u(v(t, x)) ~ Dxx(v(t, x))] + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys_multi(eqs, [u(0, x) ~ 0.0, v(0, x) ~ 0.0, u(t, 0) ~ v(t, 1)], [u(t, x), v(t, x)]), + disc, + ) + + Dx1, Dxx1 = Differential(x1), Differential(x1)^2 + Dx2, Dxx2 = Differential(x2), Differential(x2)^2 + interface_eqs = [ + Dt_u(u(t, x1)) ~ Dxx1(u(t, x1)), + Dt_u(v(t, x2)) ~ Dxx2(v(t, x2)), + ] + interface_bcs = [ + u(0, x1) ~ 1.0, + v(0, x2) ~ 0.0, + u(t, 0.0) ~ 1.0, + v(t, 2.0) ~ 0.0, + u(t, 1.0) ~ v(t, 1.0), + Dx1(u(t, 1.0)) ~ Dx2(v(t, 1.0)), + ] + interface_domains = [ + t ∈ Interval(0.0, 1.0), + x1 ∈ Interval(0.0, 1.0), + x2 ∈ Interval(1.0, 2.0), + ] + @named interface_pdesys = PDESystem( + interface_eqs, interface_bcs, interface_domains, [t, x1, x2], + [u(t, x1), v(t, x2)], + ) + interface_disc = MOLFiniteDifference([x1 => 0.1, x2 => 0.1], t) + @test run_boundary_validation(interface_pdesys, interface_disc) === nothing + end + + @testset "Multidimensional and mixed derivatives" begin + Dt = Differential(t) + Dxx, Dyy = Differential(x)^2, Differential(y)^2 + Dxxy = Differential(x)^2 * Differential(y) + _, disc = standard_2d_disc() + + eq_heat = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) + Dyy(u(t, x, y)) + ic = u(0, x, y) ~ 0.0 + bx0 = u(t, 0, y) ~ 0.0 + bx1 = u(t, 1, y) ~ 0.0 + by0 = u(t, x, 0) ~ 0.0 + by1 = u(t, x, 1) ~ 0.0 + + @test_throws ArgumentError run_boundary_validation( + build_2d_pdesys(eq_heat, [ic, bx0, by0, by1]), disc) + @test_throws ArgumentError run_boundary_validation( + build_2d_pdesys(eq_heat, [ic, bx0, bx1, by0]), disc) + @test run_boundary_validation( + build_2d_pdesys(eq_heat, [ic, bx0, bx1, by0, by1]), disc) === nothing + + eq_mixed = Dt(u(t, x, y)) ~ Dxxy(u(t, x, y)) + @test_throws ArgumentError run_boundary_validation( + build_2d_pdesys(eq_mixed, [ic, bx0, by0]), disc) + @test run_boundary_validation( + build_2d_pdesys(eq_mixed, [ic, bx0, bx1, by0]), disc) === nothing + end + + @testset "Time-dependent boundaries" begin + Dt, Dxx = Differential(t), Differential(x)^2 + eq = Dt(u(t, x)) ~ Dxx(u(t, x)) + _, disc = standard_1d_disc() + + @test run_boundary_validation( + build_1d_pdesys(eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ sin(t) * exp(-t), + u(t, 1) ~ 0.0, + ]), disc) === nothing + end + + @testset "Nested nonlinear advection" begin + Dt, Dx = Differential(t), Differential(x) + eq = Dt(u(t, x)) ~ -sin(u(t, x)^2 + 1) * Dx(u(t, x)) + _, disc = standard_1d_disc(; advection = UpwindScheme()) + + @test_throws ArgumentError run_boundary_validation( + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) + @test run_boundary_validation( + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) === nothing + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 82c4486ff..9b3369d77 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,6 +48,9 @@ run_tests(; @safetestset "MOLFiniteDifference Interface: Staggered constructors" begin include(joinpath(@__DIR__, "Components", "staggered_constructors.jl")) end + @safetestset "Boundary validation" begin + include(joinpath(@__DIR__, "Components", "boundary_validation_test.jl")) + end return @safetestset "Discrete Callbacks" begin include(joinpath(@__DIR__, "Components", "callbacks.jl")) end @@ -76,4 +79,4 @@ run_tests(; ), qa = (; env = joinpath(@__DIR__, "qa"), body = joinpath(@__DIR__, "qa", "qa.jl")), all = FUNCTIONAL_GROUPS, -) +) \ No newline at end of file From ffcb63efa85fccb05a1668f79b9fd17e7478c203 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Utku=20Y=C4=B1lmaz?= Date: Fri, 19 Jun 2026 12:49:20 +0300 Subject: [PATCH 2/2] Fix missing newline in runtests.jl Add missing newline at end of file. --- src/discretization/boundary_validation.jl | 63 +++++---- test/Components/boundary_validation_test.jl | 146 ++++++++++++-------- test/runtests.jl | 2 +- 3 files changed, 128 insertions(+), 83 deletions(-) diff --git a/src/discretization/boundary_validation.jl b/src/discretization/boundary_validation.jl index 31fc925b6..1699d82ec 100644 --- a/src/discretization/boundary_validation.jl +++ b/src/discretization/boundary_validation.jl @@ -29,15 +29,17 @@ function validate_system_wellposedness(pdes, bmap, s, _disc) has_dynamic_advection = false for eq in pdes - max_order = max(max_order, + max_order = max( + max_order, get_max_derivative_order(eq.lhs, u_op, x), - get_max_derivative_order(eq.rhs, u_op, x)) + get_max_derivative_order(eq.rhs, u_op, x) + ) end if max_order == 1 for eq in pdes if is_dynamic_advection(eq.lhs, u_op, x) || - is_dynamic_advection(eq.rhs, u_op, x) + is_dynamic_advection(eq.rhs, u_op, x) has_dynamic_advection = true break end @@ -51,33 +53,42 @@ function validate_system_wellposedness(pdes, bmap, s, _disc) if max_order >= 2 effective_bcs = count_effective_boundary_conditions( - u_bcs_x, u_call, u_op, x) + u_bcs_x, u_call, u_op, x + ) if effective_bcs < max_order - throw(ArgumentError( - "Ill-posed PDE: $(u_op) in $(x) has order-$(max_order) " * - "spatial derivatives but only $(effective_bcs) effective " * - "boundary constraint(s) provided; need $(max_order)." - )) + throw( + ArgumentError( + "Ill-posed PDE: $(u_op) in $(x) has order-$(max_order) " * + "spatial derivatives but only $(effective_bcs) effective " * + "boundary constraint(s) provided; need $(max_order)." + ) + ) end elseif max_order == 1 && has_dynamic_advection provided_locations = count_boundary_locations( - u_bcs_x, u_call, u_op, x) + u_bcs_x, u_call, u_op, x + ) if provided_locations < 2 - throw(ArgumentError( - "Ill-posed PDE: $(u_op) in $(x) has a non-constant advection " * - "coefficient, so the upwind direction is not fixed and boundary " * - "data is required at BOTH ends; only $(provided_locations) provided." - )) + throw( + ArgumentError( + "Ill-posed PDE: $(u_op) in $(x) has a non-constant advection " * + "coefficient, so the upwind direction is not fixed and boundary " * + "data is required at BOTH ends; only $(provided_locations) provided." + ) + ) end elseif max_order == 1 && !has_dynamic_advection && provided_bcs < 1 - throw(ArgumentError( - "Ill-posed PDE: $(u_op) in $(x) has a 1st-order spatial " * - "derivative but no boundary condition was provided." - )) + throw( + ArgumentError( + "Ill-posed PDE: $(u_op) in $(x) has a 1st-order spatial " * + "derivative but no boundary condition was provided." + ) + ) end end end end + return end """ @@ -147,7 +158,7 @@ function count_effective_boundary_conditions(u_bcs_x, u_call, u_op, x) function structurally_equal(a, b) if hasproperty(a, :lhs) && hasproperty(a, :rhs) && - hasproperty(b, :lhs) && hasproperty(b, :rhs) + hasproperty(b, :lhs) && hasproperty(b, :rhs) return isequal(a.lhs, b.lhs) && isequal(a.rhs, b.rhs) end return isequal(a, b) @@ -183,10 +194,12 @@ function extract_boundary_locations!(unique_bounds::Set, expr, u_op, x_index, x) normalized_val = raw_val isa Number ? Float64(raw_val) : raw_val push!(unique_bounds, normalized_val) else - throw(ArgumentError( - "Malformed boundary condition: $(u_op) needs at least $(x_index) " * - "arguments for dimension $(x), got $(length(args))." - )) + throw( + ArgumentError( + "Malformed boundary condition: $(u_op) needs at least $(x_index) " * + "arguments for dimension $(x), got $(length(args))." + ) + ) end return end @@ -286,7 +299,7 @@ function is_dynamic_advection(expr, u_op, x) arg_unwrapped = unwrap(arg) arg_op = iscall(arg_unwrapped) ? operation(arg_unwrapped) : nothing is_deriv = arg_op isa Differential && isequal(arg_op.x, x) && - get_max_derivative_order(arg_unwrapped, u_op, arg_op.x) == 1 + get_max_derivative_order(arg_unwrapped, u_op, arg_op.x) == 1 if !is_deriv && !is_static_numeric_constant(arg_unwrapped) return true end diff --git a/test/Components/boundary_validation_test.jl b/test/Components/boundary_validation_test.jl index 057dab23c..9534ad65b 100644 --- a/test/Components/boundary_validation_test.jl +++ b/test/Components/boundary_validation_test.jl @@ -17,11 +17,11 @@ function run_boundary_validation(pdesys, disc) end function standard_1d_disc(; - xspan = (0.0, 1.0), - tspan = (0.0, 1.0), - dx = 0.1, - advection = nothing, -) + xspan = (0.0, 1.0), + tspan = (0.0, 1.0), + dx = 0.1, + advection = nothing, + ) t0, t1 = tspan x0, x1 = xspan disc_kwargs = advection === nothing ? () : (; advection_scheme = advection) @@ -43,12 +43,12 @@ function build_1d_pdesys_multi(eqs, bcs, depvars; xspan = (0.0, 1.0), tspan = (0 end function standard_2d_disc(; - xspan = (0.0, 1.0), - yspan = (0.0, 1.0), - tspan = (0.0, 1.0), - dx = 0.1, - dy = 0.1, -) + xspan = (0.0, 1.0), + yspan = (0.0, 1.0), + tspan = (0.0, 1.0), + dx = 0.1, + dy = 0.1, + ) t0, t1 = tspan x0, x1b = xspan y0, y1b = yspan @@ -75,10 +75,12 @@ end _, disc = standard_1d_disc(; advection = UpwindScheme()) @test run_boundary_validation( - build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) === nothing + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc + ) === nothing @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(eq, [u(0, x) ~ 0.0]), disc) + build_1d_pdesys(eq, [u(0, x) ~ 0.0]), disc + ) end @testset "1st-order dynamic advection" begin @@ -91,13 +93,16 @@ end disc, ) === nothing @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(eq_burgers, [u(0, x) ~ 1.0, u(t, 0) ~ 1.0]), disc) + build_1d_pdesys(eq_burgers, [u(0, x) ~ 1.0, u(t, 0) ~ 1.0]), disc + ) eq_param = Dt(u(t, x)) ~ -vel * Dx(u(t, x)) @test run_boundary_validation( - build_1d_pdesys(eq_param, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) === nothing + build_1d_pdesys(eq_param, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc + ) === nothing @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(eq_param, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) + build_1d_pdesys(eq_param, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc + ) end @testset "2nd-order spatial derivatives" begin @@ -106,9 +111,11 @@ end _, disc = standard_1d_disc() @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc + ) @test run_boundary_validation( - build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) === nothing + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc + ) === nothing end @testset "3rd-order spatial derivatives" begin @@ -117,14 +124,18 @@ end _, disc = standard_1d_disc() @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc + ) @test run_boundary_validation( - build_1d_pdesys(eq, [ - u(0, x) ~ 0.0, - u(t, 0) ~ 0.0, - u(t, 1) ~ 0.0, - Dx(u(t, 0)) ~ 0.0, - ]), disc) === nothing + build_1d_pdesys( + eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ 0.0, + u(t, 1) ~ 0.0, + Dx(u(t, 0)) ~ 0.0, + ] + ), disc + ) === nothing end @testset "4th-order spatial derivatives" begin @@ -133,20 +144,26 @@ end _, disc = standard_1d_disc() @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(eq, [ - u(0, x) ~ 0.0, - u(t, 0) ~ 0.0, - u(t, 1) ~ 0.0, - Dx(u(t, 0)) ~ 0.0, - ]), disc) + build_1d_pdesys( + eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ 0.0, + u(t, 1) ~ 0.0, + Dx(u(t, 0)) ~ 0.0, + ] + ), disc + ) @test run_boundary_validation( - build_1d_pdesys(eq, [ - u(0, x) ~ 0.0, - u(t, 0) ~ 0.0, - Dx(u(t, 0)) ~ 0.0, - u(t, 1) ~ 0.0, - Dx(u(t, 1)) ~ 0.0, - ]), disc) === nothing + build_1d_pdesys( + eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ 0.0, + Dx(u(t, 0)) ~ 0.0, + u(t, 1) ~ 0.0, + Dx(u(t, 1)) ~ 0.0, + ] + ), disc + ) === nothing end @testset "Periodic boundaries" begin @@ -156,17 +173,22 @@ end periodic_bcs = [u(0, x) ~ 0.0, u(t, 0) ~ u(t, 1)] @test run_boundary_validation( - build_1d_pdesys(Dt(u(t, x)) ~ Dxx(u(t, x)), periodic_bcs), disc) === nothing + build_1d_pdesys(Dt(u(t, x)) ~ Dxx(u(t, x)), periodic_bcs), disc + ) === nothing @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(Dt(u(t, x)) ~ -Dxxxx(u(t, x)), periodic_bcs), disc) + build_1d_pdesys(Dt(u(t, x)) ~ -Dxxxx(u(t, x)), periodic_bcs), disc + ) @test run_boundary_validation( - build_1d_pdesys(Dt(u(t, x)) ~ -Dxxxx(u(t, x)), [ - u(0, x) ~ 0.0, - u(t, 0) ~ u(t, 1), - Dx(u(t, 0)) ~ Dx(u(t, 1)), - ]), disc) === nothing + build_1d_pdesys( + Dt(u(t, x)) ~ -Dxxxx(u(t, x)), [ + u(0, x) ~ 0.0, + u(t, 0) ~ u(t, 1), + Dx(u(t, 0)) ~ Dx(u(t, 1)), + ] + ), disc + ) === nothing end @testset "Coupled systems" begin @@ -220,17 +242,22 @@ end by1 = u(t, x, 1) ~ 0.0 @test_throws ArgumentError run_boundary_validation( - build_2d_pdesys(eq_heat, [ic, bx0, by0, by1]), disc) + build_2d_pdesys(eq_heat, [ic, bx0, by0, by1]), disc + ) @test_throws ArgumentError run_boundary_validation( - build_2d_pdesys(eq_heat, [ic, bx0, bx1, by0]), disc) + build_2d_pdesys(eq_heat, [ic, bx0, bx1, by0]), disc + ) @test run_boundary_validation( - build_2d_pdesys(eq_heat, [ic, bx0, bx1, by0, by1]), disc) === nothing + build_2d_pdesys(eq_heat, [ic, bx0, bx1, by0, by1]), disc + ) === nothing eq_mixed = Dt(u(t, x, y)) ~ Dxxy(u(t, x, y)) @test_throws ArgumentError run_boundary_validation( - build_2d_pdesys(eq_mixed, [ic, bx0, by0]), disc) + build_2d_pdesys(eq_mixed, [ic, bx0, by0]), disc + ) @test run_boundary_validation( - build_2d_pdesys(eq_mixed, [ic, bx0, bx1, by0]), disc) === nothing + build_2d_pdesys(eq_mixed, [ic, bx0, bx1, by0]), disc + ) === nothing end @testset "Time-dependent boundaries" begin @@ -239,11 +266,14 @@ end _, disc = standard_1d_disc() @test run_boundary_validation( - build_1d_pdesys(eq, [ - u(0, x) ~ 0.0, - u(t, 0) ~ sin(t) * exp(-t), - u(t, 1) ~ 0.0, - ]), disc) === nothing + build_1d_pdesys( + eq, [ + u(0, x) ~ 0.0, + u(t, 0) ~ sin(t) * exp(-t), + u(t, 1) ~ 0.0, + ] + ), disc + ) === nothing end @testset "Nested nonlinear advection" begin @@ -252,8 +282,10 @@ end _, disc = standard_1d_disc(; advection = UpwindScheme()) @test_throws ArgumentError run_boundary_validation( - build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc) + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0]), disc + ) @test run_boundary_validation( - build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc) === nothing + build_1d_pdesys(eq, [u(0, x) ~ 0.0, u(t, 0) ~ 0.0, u(t, 1) ~ 0.0]), disc + ) === nothing end end diff --git a/test/runtests.jl b/test/runtests.jl index 9b3369d77..83ef9f0c6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -79,4 +79,4 @@ run_tests(; ), qa = (; env = joinpath(@__DIR__, "qa"), body = joinpath(@__DIR__, "qa", "qa.jl")), all = FUNCTIONAL_GROUPS, -) \ No newline at end of file +)