From 54708cac419e38a824fed2f9862ffdbafe42a5a8 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Wed, 8 Jul 2026 02:48:49 +0300 Subject: [PATCH 1/3] fix(tests): resolve O(dx^2) initial condition residual in DAE tests --- test/DAE/MOL_1D_PDAE.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/DAE/MOL_1D_PDAE.jl b/test/DAE/MOL_1D_PDAE.jl index a46a33d76..b65bd1986 100644 --- a/test/DAE/MOL_1D_PDAE.jl +++ b/test/DAE/MOL_1D_PDAE.jl @@ -2,6 +2,7 @@ # Packages and inclusions using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets +using DiffEqBase: BrownFullBasicInit using OrdinaryDiffEqRosenbrock: Rodas4 using ModelingToolkit: Differential @@ -51,7 +52,9 @@ using ModelingToolkit: Differential # Convert the PDE problem into an ODE problem prob = discretize(pdesys, discretization) # Solve ODE problem - sol = solve(prob, Rodas4(), saveat = 0.1) + # Semi-discrete algebraic IC v(0,x)=sin(x) satisfies Dxx(v)+f only to O(dx^2); + # BrownFullBasicInit projects algebraic variables onto the consistent manifold. + sol = solve(prob, Rodas4(), saveat = 0.1, initializealg = BrownFullBasicInit()) x_sol = sol[x] t_sol = sol[t] From 71fffbf4c022dff8d3697afa4f06ce40af0bae69 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Wed, 8 Jul 2026 20:03:12 +0300 Subject: [PATCH 2/3] feat(discretization): implement automatic BrownFullBasicInit fallback for implicit DAEs --- docs/src/MOLFiniteDifference.md | 17 +++++ src/MethodOfLines.jl | 3 + src/discretization/dae_detection.jl | 23 ++++++ src/discretization/dae_init_fallback.jl | 35 +++++++++ src/discretization/discretize.jl | 83 ++++++++++++++++++++++ src/discretization/staggered_discretize.jl | 3 +- test/Complex/schroedinger.jl | 5 +- test/Components/dae_init_fallback.jl | 76 ++++++++++++++++++++ test/DAE/MOL_1D_PDAE.jl | 5 +- test/runtests.jl | 3 + 10 files changed, 245 insertions(+), 8 deletions(-) create mode 100644 src/discretization/dae_detection.jl create mode 100644 src/discretization/dae_init_fallback.jl create mode 100644 src/discretization/discretize.jl create mode 100644 test/Components/dae_init_fallback.jl diff --git a/docs/src/MOLFiniteDifference.md b/docs/src/MOLFiniteDifference.md index ad418fc4a..ba57c7d66 100644 --- a/docs/src/MOLFiniteDifference.md +++ b/docs/src/MOLFiniteDifference.md @@ -47,5 +47,22 @@ 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` produces an implicit DAE (`0 ~ ...` constraints, over-determined boundary +conditions, etc.), MethodOfLines automatically attaches `initializealg = BrownFullBasicInit()` +to the returned `ODEProblem.kwargs`. This makes `solve(prob, alg)` work under OrdinaryDiffEq +v7's default `CheckInit` path without passing `initializealg` explicitly. + +Override at any level (highest priority wins at `solve` time): + +- `solve(prob, alg; initializealg = MyInit())` +- `discretize(pdesys, disc; initializealg = MyInit())` +- `MOLFiniteDifference(dxs, t; initializealg = MyInit())` + +Pure ODE discretizations are unchanged. The fallback is not applied when building an +`ODEFunction` directly via `ODEFunction(pdesys, discretization)`; construct an +`ODEProblem` from that function 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..0f344aaab 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -108,6 +108,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..aad2b0ca3 --- /dev/null +++ b/src/discretization/dae_detection.jl @@ -0,0 +1,23 @@ +""" + is_implicit_dae(prob::ODEProblem) -> Bool + +Return `true` when `prob` represents an implicit ODE/DAE system that requires +consistent initialization (singular mass matrix or MTK initialization data). +Matches the DAE detection used in OrdinaryDiffEq when `initialize_dae!` runs. +""" +function is_implicit_dae(prob::SciMLBase.ODEProblem) + SciMLBase.has_initializeprob(prob.f) && return true + return is_singular_mass_matrix(prob.f.mass_matrix) +end + +function is_singular_mass_matrix(M) + M === I && return false + M isa UniformScaling && return false + M isa Tuple && return false + if M isa Diagonal + return any(iszero, M.diag) + elseif M isa AbstractMatrix + return any(iszero, diag(M)) + end + return false +end diff --git a/src/discretization/dae_init_fallback.jl b/src/discretization/dae_init_fallback.jl new file mode 100644 index 000000000..004d21687 --- /dev/null +++ b/src/discretization/dae_init_fallback.jl @@ -0,0 +1,35 @@ +using DiffEqBase: BrownFullBasicInit + +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...) + +If `prob` is an implicit DAE and the user did not supply `initializealg`, +attach `BrownFullBasicInit()` to `prob.kwargs` so `solve(prob, alg)` succeeds +under OrdinaryDiffEq v7's default `CheckInit` path. +""" +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 + is_implicit_dae(prob) || return prob + return SciMLBase.remake(prob; initializealg = BrownFullBasicInit()) +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..f2897eb48 --- /dev/null +++ b/src/discretization/discretize.jl @@ -0,0 +1,83 @@ +using ModelingToolkit: get_ps, mtkcompile, ProblemTypeCtx, getmetadata + +function _safe_unwrap(x) + return x isa Num ? unwrap(x) : x +end + +function _build_ode_problem( + simpsys, mol_metadata, tspan, discretization::MOLFiniteDifference, kwargs... + ) + u0 = hasproperty(mol_metadata, :u0) ? mol_metadata.u0 : [] + pdesys_ic = mol_metadata.pdesys.initial_conditions + ps_raw = get_ps(mol_metadata.pdesys) + prob_kwargs = (; + build_initializeprob = false, + discretization.kwargs..., + kwargs..., + ) + if ps_raw !== nothing && ps_raw !== SciMLBase.NullParameters() && !isempty(ps_raw) + param_vals = Dict{Any, Any}() + if first(ps_raw) isa Pair + for p in ps_raw + param_vals[first(p)] = last(p) + end + else + ps_unwrapped = [_safe_unwrap(p) for p in ps_raw] + for (k, v) in pairs(pdesys_ic) + k_unwrapped = _safe_unwrap(k) + if any(p -> isequal(k_unwrapped, _safe_unwrap(p)), ps_unwrapped) + v_numeric = try + Symbolics.value(v) + catch + _safe_unwrap(v) + end + param_vals[k] = v_numeric + end + end + end + if !isempty(param_vals) + op = merge(Dict(u0), param_vals) + return ODEProblem(simpsys, op, tspan; prob_kwargs...) + end + end + return ODEProblem(simpsys, u0, tspan; prob_kwargs...) +end + +function SciMLBase.discretize( + pdesys::PDESystem, + discretization::MOLFiniteDifference{G, D}; + analytic = nothing, checks = true, kwargs... + ) where {G, D <: ScalarizedDiscretization} + sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization; checks = checks) + return try + simpsys = mtkcompile(sys) + if tspan === nothing + add_metadata!(getmetadata(sys, ProblemTypeCtx, nothing), sys) + unknowns_list = ModelingToolkit.unknowns(simpsys) + u0_guess = Dict(u => 1.0 for u in unknowns_list) + return NonlinearProblem( + simpsys, u0_guess; + discretization.kwargs..., kwargs... + ) + else + mol_metadata = getmetadata(simpsys, ProblemTypeCtx, nothing) + add_metadata!(mol_metadata, sys) + prob = _build_ode_problem(simpsys, mol_metadata, tspan, discretization, kwargs...) + if analytic === nothing + return apply_dae_initialization_fallback(prob, discretization; kwargs...) + else + f = SciMLBase.ODEFunction( + pdesys, discretization, analytic = analytic, + discretization.kwargs..., kwargs... + ) + prob = ODEProblem( + f, prob.u0, prob.tspan, prob.p; + discretization.kwargs..., kwargs... + ) + return apply_dae_initialization_fallback(prob, discretization; kwargs...) + end + end + catch e + error_analysis(sys, e) + end +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..0c77fe76c --- /dev/null +++ b/test/Components/dae_init_fallback.jl @@ -0,0 +1,76 @@ +using MethodOfLines, ModelingToolkit, DomainSets, Test, OrdinaryDiffEq, SciMLBase +using DiffEqBase: BrownFullBasicInit, CheckInit +using ModelingToolkit: Differential + +@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/DAE/MOL_1D_PDAE.jl b/test/DAE/MOL_1D_PDAE.jl index b65bd1986..a46a33d76 100644 --- a/test/DAE/MOL_1D_PDAE.jl +++ b/test/DAE/MOL_1D_PDAE.jl @@ -2,7 +2,6 @@ # Packages and inclusions using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets -using DiffEqBase: BrownFullBasicInit using OrdinaryDiffEqRosenbrock: Rodas4 using ModelingToolkit: Differential @@ -52,9 +51,7 @@ using ModelingToolkit: Differential # Convert the PDE problem into an ODE problem prob = discretize(pdesys, discretization) # Solve ODE problem - # Semi-discrete algebraic IC v(0,x)=sin(x) satisfies Dxx(v)+f only to O(dx^2); - # BrownFullBasicInit projects algebraic variables onto the consistent manifold. - sol = solve(prob, Rodas4(), saveat = 0.1, initializealg = BrownFullBasicInit()) + sol = solve(prob, Rodas4(), saveat = 0.1) x_sol = sol[x] t_sol = sol[t] 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 From 6c9b740b76dfdc005674577071e061613df5dd98 Mon Sep 17 00:00:00 2001 From: utkuyilmaz1903 Date: Thu, 9 Jul 2026 11:05:56 +0300 Subject: [PATCH 3/3] fix(discretize): resolve ambiguity via invoke wrapper and use numeric-only DAE detection --- src/MethodOfLines.jl | 1 + src/discretization/dae_detection.jl | 11 ++-- src/discretization/dae_init_fallback.jl | 8 +-- src/discretization/discretize.jl | 86 +++---------------------- test/Components/dae_init_fallback.jl | 15 ++++- 5 files changed, 34 insertions(+), 87 deletions(-) diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 0f344aaab..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, diff --git a/src/discretization/dae_detection.jl b/src/discretization/dae_detection.jl index aad2b0ca3..3b5b840c5 100644 --- a/src/discretization/dae_detection.jl +++ b/src/discretization/dae_detection.jl @@ -1,12 +1,11 @@ """ is_implicit_dae(prob::ODEProblem) -> Bool -Return `true` when `prob` represents an implicit ODE/DAE system that requires -consistent initialization (singular mass matrix or MTK initialization data). -Matches the DAE detection used in OrdinaryDiffEq when `initialize_dae!` runs. +Return `true` when `prob` represents an implicit DAE with a singular numeric mass +matrix. Intended for MTK-compiled numeric `ODEProblem`s from `discretize`; symbolic +or non-numeric mass matrices conservatively return `false`. """ function is_implicit_dae(prob::SciMLBase.ODEProblem) - SciMLBase.has_initializeprob(prob.f) && return true return is_singular_mass_matrix(prob.f.mass_matrix) end @@ -15,8 +14,12 @@ function is_singular_mass_matrix(M) M isa UniformScaling && return false M isa Tuple && return false if M isa Diagonal + any(d -> d isa Num, M.diag) && return false + eltype(M.diag) <: Number || return false return any(iszero, M.diag) elseif M isa AbstractMatrix + any(d -> d isa Num, M) && return false + eltype(M) <: Number || return false return any(iszero, diag(M)) end return false diff --git a/src/discretization/dae_init_fallback.jl b/src/discretization/dae_init_fallback.jl index 004d21687..1944d7570 100644 --- a/src/discretization/dae_init_fallback.jl +++ b/src/discretization/dae_init_fallback.jl @@ -1,5 +1,3 @@ -using DiffEqBase: BrownFullBasicInit - function _disc_kwargs_nt(disc::MOLFiniteDifference) kw = disc.kwargs return kw isa NamedTuple ? kw : NamedTuple() @@ -28,8 +26,10 @@ function apply_dae_initialization_fallback( if _user_provided_initializealg(discretization, discretize_kwargs, prob) return prob end - is_implicit_dae(prob) || return prob - return SciMLBase.remake(prob; initializealg = BrownFullBasicInit()) + 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 index f2897eb48..53db1f8ef 100644 --- a/src/discretization/discretize.jl +++ b/src/discretization/discretize.jl @@ -1,83 +1,13 @@ -using ModelingToolkit: get_ps, mtkcompile, ProblemTypeCtx, getmetadata - -function _safe_unwrap(x) - return x isa Num ? unwrap(x) : x -end - -function _build_ode_problem( - simpsys, mol_metadata, tspan, discretization::MOLFiniteDifference, kwargs... - ) - u0 = hasproperty(mol_metadata, :u0) ? mol_metadata.u0 : [] - pdesys_ic = mol_metadata.pdesys.initial_conditions - ps_raw = get_ps(mol_metadata.pdesys) - prob_kwargs = (; - build_initializeprob = false, - discretization.kwargs..., - kwargs..., - ) - if ps_raw !== nothing && ps_raw !== SciMLBase.NullParameters() && !isempty(ps_raw) - param_vals = Dict{Any, Any}() - if first(ps_raw) isa Pair - for p in ps_raw - param_vals[first(p)] = last(p) - end - else - ps_unwrapped = [_safe_unwrap(p) for p in ps_raw] - for (k, v) in pairs(pdesys_ic) - k_unwrapped = _safe_unwrap(k) - if any(p -> isequal(k_unwrapped, _safe_unwrap(p)), ps_unwrapped) - v_numeric = try - Symbolics.value(v) - catch - _safe_unwrap(v) - end - param_vals[k] = v_numeric - end - end - end - if !isempty(param_vals) - op = merge(Dict(u0), param_vals) - return ODEProblem(simpsys, op, tspan; prob_kwargs...) - end - end - return ODEProblem(simpsys, u0, tspan; prob_kwargs...) -end - function SciMLBase.discretize( pdesys::PDESystem, discretization::MOLFiniteDifference{G, D}; analytic = nothing, checks = true, kwargs... - ) where {G, D <: ScalarizedDiscretization} - sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization; checks = checks) - return try - simpsys = mtkcompile(sys) - if tspan === nothing - add_metadata!(getmetadata(sys, ProblemTypeCtx, nothing), sys) - unknowns_list = ModelingToolkit.unknowns(simpsys) - u0_guess = Dict(u => 1.0 for u in unknowns_list) - return NonlinearProblem( - simpsys, u0_guess; - discretization.kwargs..., kwargs... - ) - else - mol_metadata = getmetadata(simpsys, ProblemTypeCtx, nothing) - add_metadata!(mol_metadata, sys) - prob = _build_ode_problem(simpsys, mol_metadata, tspan, discretization, kwargs...) - if analytic === nothing - return apply_dae_initialization_fallback(prob, discretization; kwargs...) - else - f = SciMLBase.ODEFunction( - pdesys, discretization, analytic = analytic, - discretization.kwargs..., kwargs... - ) - prob = ODEProblem( - f, prob.u0, prob.tspan, prob.p; - discretization.kwargs..., kwargs... - ) - return apply_dae_initialization_fallback(prob, discretization; kwargs...) - end - end - catch e - error_analysis(sys, e) - end + ) 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/test/Components/dae_init_fallback.jl b/test/Components/dae_init_fallback.jl index 0c77fe76c..5749242bc 100644 --- a/test/Components/dae_init_fallback.jl +++ b/test/Components/dae_init_fallback.jl @@ -1,6 +1,19 @@ -using MethodOfLines, ModelingToolkit, DomainSets, Test, OrdinaryDiffEq, SciMLBase +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) +end @testset "is_implicit_dae" begin @parameters t x