fix(tests): resolve O(dx^2) initial condition residual in DAE tests

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.
 

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

src/discretization/dae_detection.jl

@@ -0,0 +1,26 @@
+"""
+    is_implicit_dae(prob::ODEProblem) -> Bool
+
+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)
+    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
+        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
+end

src/discretization/dae_init_fallback.jl

@@ -0,0 +1,35 @@
+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
+    if is_implicit_dae(prob)
+        return SciMLBase.remake(prob; initializealg = BrownFullBasicInit())
+    end
+    return prob
+end
+
+apply_dae_initialization_fallback(prob, ::MOLFiniteDifference; kwargs...) = prob

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

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)

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

test/Components/dae_init_fallback.jl

@@ -0,0 +1,89 @@
+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
+    @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

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