diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index b7f781658..08b06532a 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: package: - - {user: SciML, repo: PDESystemLibrary.jl, group: MOL} + - {user: SciML, repo: PDESystemLibrary.jl, group: Core} uses: "SciML/.github/.github/workflows/downstream.yml@v1" with: owner: "${{ matrix.package.user }}" diff --git a/docs/src/MOLFiniteDifference.md b/docs/src/MOLFiniteDifference.md index ad418fc4a..b29bd02c4 100644 --- a/docs/src/MOLFiniteDifference.md +++ b/docs/src/MOLFiniteDifference.md @@ -47,5 +47,23 @@ Currently supported options are `grid_align`: `center_align` and `edge_align`. E `use_ODAE`: MethodOfLines will automatically make use of `ODAEProblem` where relevant, which improves performance for DAEs (as discretized PDEs are in general), if this is set to true. Defaults to false. +### DAE / PDAE initialization + +When `discretize` returns a numeric `ODEProblem` with a singular mass matrix or MTK +initialization metadata, MethodOfLines attaches `initializealg = BrownFullBasicInit()` +to `ODEProblem.kwargs`. This matches the DAE signals used in OrdinaryDiffEq v7 and allows +`solve(prob, alg)` under the default `CheckInit` path without passing `initializealg` +explicitly. + +Override at any level (highest priority at `solve`): + +- `solve(prob, alg; initializealg = MyInit())` +- `discretize(pdesys, disc; initializealg = MyInit())` +- `MOLFiniteDifference(dxs, t; initializealg = MyInit())` + +Pure ODE discretizations (`mass_matrix = I`, no initialization metadata) are unchanged. +The fallback is not applied for `ODEFunction(pdesys, discretization)`; construct an +`ODEProblem` if you need the same behavior. + Any unrecognized keyword arguments will be passed to the `ODEProblem` constructor, see [its documentation](https://docs.sciml.ai/ModelingToolkit/stable/API/problems/#Dynamical-systems) for available options. diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 3feb8c796..26adfa1d6 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -2,6 +2,7 @@ module MethodOfLines using LinearAlgebra using SciMLBase using DiffEqBase +using DiffEqBase: BrownFullBasicInit using ModelingToolkit using ModelingToolkit: get_unknowns, get_eqs, get_bcs, get_dvs, @@ -108,6 +109,9 @@ include("discretization/schemes/integral_expansion/integral_expansion.jl") include("discretization/generate_finite_difference_rules.jl") include("discretization/generate_bc_eqs.jl") include("discretization/generate_ic_defaults.jl") +include("discretization/dae_detection.jl") +include("discretization/dae_init_fallback.jl") +include("discretization/discretize.jl") include("discretization/staggered_discretize.jl") # Main diff --git a/src/discretization/dae_detection.jl b/src/discretization/dae_detection.jl new file mode 100644 index 000000000..2f6752e73 --- /dev/null +++ b/src/discretization/dae_detection.jl @@ -0,0 +1,84 @@ +function safe_is_numeric_zero(x) + x isa Num && return false + x isa Bool && return !x + x isa Number && return iszero(x) + return false +end + +_is_numeric_mass_entry(x) = x isa Number || x isa Bool + +function _diagonal_mass_is_singular(diag) + @inbounds for x in diag + if safe_is_numeric_zero(x) + return true + end + end + return false +end + +function _state_vector_numeric(u0) + u0 isa Num && return false + u0 isa Number && return true + u0 isa AbstractArray || return false + @inbounds for u in u0 + u isa Num && return false + if !(u isa Number) && !(u isa Bool) + return false + end + end + return true +end + +function _mass_matrix_numeric_or_identity(M) + M === I && return true + M isa UniformScaling && return _is_numeric_mass_entry(M.λ) + M isa Tuple && return false + if M isa Diagonal + @inbounds for x in M.diag + _is_numeric_mass_entry(x) || return false + end + return true + elseif M isa AbstractMatrix + @inbounds for x in M + _is_numeric_mass_entry(x) || return false + end + return true + end + return false +end + +function numeric_contract(prob::SciMLBase.ODEProblem) + return _state_vector_numeric(prob.u0) && + _mass_matrix_numeric_or_identity(prob.f.mass_matrix) +end + +function safe_has_initializeprob(f) + return SciMLBase.has_initializeprob(f) === true +end + +function safe_isdae_mass(M) + M === I && return false + M isa UniformScaling && return false + M isa Tuple && return false + if M isa Diagonal + return _diagonal_mass_is_singular(M.diag) + elseif M isa AbstractMatrix + return _diagonal_mass_is_singular(diag(M)) + end + return false +end + +is_singular_mass_matrix(M) = safe_isdae_mass(M) + +""" + is_implicit_dae(prob::ODEProblem) -> Bool + +`true` for numeric problems with singular mass matrix or MTK initialization metadata. +""" +function is_implicit_dae(prob::SciMLBase.ODEProblem) + numeric_contract(prob) || return false + if safe_has_initializeprob(prob.f) + return true + end + return safe_isdae_mass(prob.f.mass_matrix) +end diff --git a/src/discretization/dae_init_fallback.jl b/src/discretization/dae_init_fallback.jl new file mode 100644 index 000000000..2072c73c2 --- /dev/null +++ b/src/discretization/dae_init_fallback.jl @@ -0,0 +1,33 @@ +function _disc_kwargs_nt(disc::MOLFiniteDifference) + kw = disc.kwargs + return kw isa NamedTuple ? kw : NamedTuple() +end + +function _user_provided_initializealg( + disc::MOLFiniteDifference, discretize_kwargs::NamedTuple, prob + ) + merged = merge(_disc_kwargs_nt(disc), discretize_kwargs) + return haskey(merged, :initializealg) || haskey(prob.kwargs, :initializealg) +end + +""" + apply_dae_initialization_fallback(prob, discretization; kwargs...) + +Remake with `BrownFullBasicInit()` when `is_implicit_dae(prob)` and no user `initializealg`. +""" +function apply_dae_initialization_fallback( + prob::SciMLBase.ODEProblem, + discretization::MOLFiniteDifference; + kwargs... + ) + discretize_kwargs = NamedTuple(kwargs) + if _user_provided_initializealg(discretization, discretize_kwargs, prob) + return prob + end + if is_implicit_dae(prob) + return SciMLBase.remake(prob; initializealg = BrownFullBasicInit()) + end + return prob +end + +apply_dae_initialization_fallback(prob, ::MOLFiniteDifference; kwargs...) = prob diff --git a/src/discretization/discretize.jl b/src/discretization/discretize.jl new file mode 100644 index 000000000..53db1f8ef --- /dev/null +++ b/src/discretization/discretize.jl @@ -0,0 +1,13 @@ +function SciMLBase.discretize( + pdesys::PDESystem, + discretization::MOLFiniteDifference{G, D}; + analytic = nothing, checks = true, kwargs... + ) where {G <: Union{CenterAlignedGrid, EdgeAlignedGrid}, D <: ScalarizedDiscretization} + prob = invoke( + SciMLBase.discretize, + Tuple{PDESystem, AbstractEquationSystemDiscretization}, + pdesys, discretization; + analytic, checks, kwargs... + ) + return apply_dae_initialization_fallback(prob, discretization; kwargs...) +end diff --git a/src/discretization/staggered_discretize.jl b/src/discretization/staggered_discretize.jl index 87331d041..efc8bcd4d 100644 --- a/src/discretization/staggered_discretize.jl +++ b/src/discretization/staggered_discretize.jl @@ -23,7 +23,8 @@ function SciMLBase.discretize( discretization.kwargs..., kwargs... ) - return symbolic_trace(prob, simpsys) + prob = symbolic_trace(prob, simpsys) + return apply_dae_initialization_fallback(prob, discretization; kwargs...) end catch e error_analysis(sys, e) diff --git a/test/Complex/schroedinger.jl b/test/Complex/schroedinger.jl index af84ff6cf..fc9fd6c90 100644 --- a/test/Complex/schroedinger.jl +++ b/test/Complex/schroedinger.jl @@ -1,6 +1,5 @@ using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test using SciMLBase -using DiffEqBase: BrownFullBasicInit @testset "Schroedinger" begin @parameters t, x @@ -32,7 +31,7 @@ using DiffEqBase: BrownFullBasicInit prob = discretize(sys, disc) - sol = solve(prob, FBDF(), saveat = 0.01, initializealg = BrownFullBasicInit()) + sol = solve(prob, FBDF(), saveat = 0.01) discx = sol[x] disct = sol[t] @@ -91,7 +90,7 @@ end prob = discretize(sys, disc) - sol = solve(prob, FBDF(), saveat = 0.01, initializealg = BrownFullBasicInit()) + sol = solve(prob, FBDF(), saveat = 0.01) @test SciMLBase.successful_retcode(sol) end diff --git a/test/Components/dae_init_fallback.jl b/test/Components/dae_init_fallback.jl new file mode 100644 index 000000000..45e86ac54 --- /dev/null +++ b/test/Components/dae_init_fallback.jl @@ -0,0 +1,92 @@ +using MethodOfLines, ModelingToolkit, DomainSets, Test, OrdinaryDiffEq, SciMLBase, LinearAlgebra +using DiffEqBase: BrownFullBasicInit, CheckInit +using ModelingToolkit: Differential +using Symbolics + +@testset "is_singular_mass_matrix numeric and symbolic safety" begin + @test !MethodOfLines.is_singular_mass_matrix(I) + @test !MethodOfLines.is_singular_mass_matrix(UniformScaling(true)) + @test !MethodOfLines.is_singular_mass_matrix(UniformScaling(2.0)) + @test MethodOfLines.is_singular_mass_matrix(Diagonal([1.0, 0.0, 2.0])) + @test !MethodOfLines.is_singular_mass_matrix(Diagonal([1.0, 2.0, 3.0])) + + @variables m1 m2 + M_sym = Diagonal([m1, 0, m2]) + @test !MethodOfLines.is_singular_mass_matrix(M_sym) + @test MethodOfLines.is_singular_mass_matrix(M_sym) isa Bool + @test MethodOfLines.safe_is_numeric_zero(0.0) + @test !MethodOfLines.safe_is_numeric_zero(m1) +end + +@testset "is_implicit_dae" begin + @parameters t x + @variables u(..) v(..) + Dt = Differential(t) + Dxx = Differential(x)^2 + pdae_eqs = [ + Dt(u(t, x)) ~ Dxx(u(t, x)), + 0 ~ Dxx(v(t, x)) + exp(-t) * sin(x), + ] + pdae_bcs = [ + u(0, x) ~ cos(x), + v(0, x) ~ sin(x), + u(t, 0) ~ exp(-t), + Differential(x)(u(t, 1)) ~ -exp(-t) * sin(1), + Differential(x)(v(t, 0)) ~ exp(-t), + v(t, 1) ~ exp(-t) * sin(1), + ] + ode_eq = Dt(u(t, x)) ~ Dxx(u(t, x)) + ode_bcs = [u(0, x) ~ sin(x), u(t, 0) ~ 0.0, u(t, 1) ~ 0.0] + domains = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0)] + + @named pdae_sys = PDESystem(pdae_eqs, pdae_bcs, domains, [t, x], [u(t, x), v(t, x)]) + @named ode_sys = PDESystem(ode_eq, ode_bcs, domains, [t, x], [u(t, x)]) + + dx = 1 / 19 + pdae_prob = discretize(pdae_sys, MOLFiniteDifference([x => dx], t)) + ode_prob = discretize(ode_sys, MOLFiniteDifference([x => dx], t)) + + @test MethodOfLines.is_implicit_dae(pdae_prob) + @test !MethodOfLines.is_implicit_dae(ode_prob) +end + +@testset "apply_dae_initialization_fallback" begin + @parameters t x + @variables u(..) v(..) + Dt = Differential(t) + Dxx = Differential(x)^2 + eqs = [ + Dt(u(t, x)) ~ Dxx(u(t, x)), + 0 ~ Dxx(v(t, x)) + exp(-t) * sin(x), + ] + bcs = [ + u(0, x) ~ cos(x), + v(0, x) ~ sin(x), + u(t, 0) ~ exp(-t), + Differential(x)(u(t, 1)) ~ -exp(-t) * sin(1), + Differential(x)(v(t, 0)) ~ exp(-t), + v(t, 1) ~ exp(-t) * sin(1), + ] + domains = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0)] + @named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), v(t, x)]) + disc = MOLFiniteDifference([x => 1 / 19], t) + + prob = discretize(pdesys, disc) + @test haskey(prob.kwargs, :initializealg) + @test prob.kwargs[:initializealg] isa BrownFullBasicInit + + sol = solve(prob, FBDF(), saveat = 0.1) + @test SciMLBase.successful_retcode(sol) + + @named ode_only = PDESystem( + [Dt(u(t, x)) ~ Dxx(u(t, x))], + [u(0, x) ~ sin(x), u(t, 0) ~ 0.0, u(t, 1) ~ 0.0], + domains, [t, x], [u(t, x)] + ) + prob_ode = discretize(ode_only, disc) + @test !haskey(prob_ode.kwargs, :initializealg) + + prob_check = discretize(pdesys, disc; initializealg = CheckInit()) + @test prob_check.kwargs[:initializealg] isa CheckInit + @test_throws Exception solve(prob_check, FBDF()) +end diff --git a/test/qa/qa.jl b/test/qa/qa.jl index 6a1450820..efe83e84a 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -53,7 +53,7 @@ run_qa( :AbstractCartesianIndex, # non-public in SciMLBase: :AbstractDiscretizationMetadata, :AbstractODESolution, :AbstractPDESolution, - :PDESolution, :observed, + :PDESolution, :observed, :has_initializeprob, # non-public in PDEBase: :EquationState, :cardinalize_eqs!, :check_boundarymap, :construct_differential_discretizer, :construct_disc_state, diff --git a/test/runtests.jl b/test/runtests.jl index 033ea1a58..cd652e7f5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,6 +57,9 @@ run_tests(; @safetestset "ODEFunction" begin include(joinpath(@__DIR__, "Components", "ODEFunction_test.jl")) end + @safetestset "DAE initialization fallback" begin + include(joinpath(@__DIR__, "Components", "dae_init_fallback.jl")) + end @safetestset "MOLFiniteDifference Interface: Staggered constructors" begin include(joinpath(@__DIR__, "Components", "staggered_constructors.jl")) end