Skip to content

Commit bd2fc6a

Browse files
SebastianM-Cclaude
SebastianM-C
and
claude
committed
Error on adaptive solves for DAEs and add index-1 DAE tests
Defect-based mesh adaptivity cannot converge for DAEs since the collocation interpolant is inaccurate for algebraic variables, so requesting it now throws an ArgumentError pointing users to `adaptive = false` or BoundaryValueDiffEqAscher.jl. Also fix the nested FIRK out-of-place collocation kernels calling the in-place algebraic constraint helper on the user's out-of-place f. Tests cover MIRK2-6, RadauIIa and LobattoIIIc (expanded and nested) on a simple index-1 DAE and on example problem 1 from Ascher & Spiteri (1994), plus the adaptivity error. LobattoIIIa/IIIb are excluded: their tableau structure cannot determine the algebraic stage components. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
1 parent 94a7cb6 commit bd2fc6a

11 files changed

Lines changed: 344 additions & 12 deletions

File tree

‎lib/BoundaryValueDiffEqCore/src/utils.jl‎

Lines changed: 25 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -798,25 +798,46 @@ end
798798
end
799799

800800
@inline function __apply_algebraic_constraint!(
801-
residᵢ, ::Nothing, f!, yᵢ₊₁, p, t, tmp)
801+
residᵢ, ::Nothing, f!, yᵢ₊₁, p, t, tmp
802+
)
802803
return nothing
803804
end
804805

805806
@inline function __apply_algebraic_constraint!(
806-
residᵢ, algebraic_indices::Vector{Int}, f!, yᵢ₊₁, p, t, tmp)
807+
residᵢ, algebraic_indices::Vector{Int}, f!, yᵢ₊₁, p, t, tmp
808+
)
807809
f!(tmp, yᵢ₊₁, p, t)
808810
residᵢ[algebraic_indices] .= tmp[algebraic_indices]
809811
return nothing
810812
end
811813

812814
@inline function __apply_algebraic_constraint_oop!(
813-
residᵢ, ::Nothing, f, yᵢ₊₁, p, t)
815+
residᵢ, ::Nothing, f, yᵢ₊₁, p, t
816+
)
814817
return nothing
815818
end
816819

817820
@inline function __apply_algebraic_constraint_oop!(
818-
residᵢ, algebraic_indices::Vector{Int}, f, yᵢ₊₁, p, t)
821+
residᵢ, algebraic_indices::Vector{Int}, f, yᵢ₊₁, p, t
822+
)
819823
tmp = f(yᵢ₊₁, p, t)
820824
residᵢ[algebraic_indices] .= tmp[algebraic_indices]
821825
return nothing
822826
end
827+
828+
@inline __check_dae_adaptivity(::Nothing, adaptive::Bool) = nothing
829+
830+
@inline function __check_dae_adaptivity(::Vector{Int}, adaptive::Bool)
831+
if adaptive
832+
throw(
833+
ArgumentError(
834+
"Adaptive mesh refinement is not supported for DAE problems (mass " *
835+
"matrices with zero rows): the collocation interpolant is inaccurate " *
836+
"for algebraic variables, so the defect estimate cannot converge. " *
837+
"Pass `adaptive = false`, or use a solver from " *
838+
"BoundaryValueDiffEqAscher.jl, which supports mesh adaptivity for DAEs."
839+
)
840+
)
841+
end
842+
return nothing
843+
end

‎lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
2727
__internal_solve, __default_sparsity_detector, __build_cost,
2828
__tunable_part, __add_singular_term!, __apply_mass_matrix!,
2929
__get_algebraic_indices, __subtract_mass_stage!, __apply_algebraic_constraint!,
30-
__apply_algebraic_constraint_oop!
30+
__apply_algebraic_constraint_oop!, __check_dae_adaptivity
3131

3232
using ConcreteStructs: @concrete
3333
using DifferentiationInterface: DifferentiationInterface, Constant

‎lib/BoundaryValueDiffEqFIRK/src/collocation.jl‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -302,7 +302,7 @@ function Φ(cache::FIRKCacheNested, y, u, trait)
302302
cache.fᵢ_cache, cache.k_discrete, cache.f, cache.TU, y, u,
303303
cache.p, cache.mass_matrix, cache.algebraic_indices,
304304
cache.mesh, cache.mesh_dt, cache.stage, cache, trait
305-
)
305+
)
306306
end
307307

308308
@views function Φ(
@@ -398,7 +398,6 @@ end
398398
)
399399
(; b) = TU
400400
(; nest_prob, alg) = cache
401-
tmp1 = similar(get_tmp(y[1], u))
402401

403402
residuals = [safe_similar(yᵢ) for yᵢ in y[1:(end - 1)]]
404403

@@ -422,7 +421,7 @@ end
422421

423422
@. residᵢ = yᵢ₊₁ - yᵢ
424423
__maybe_matmul!(residᵢ, nestsol.u, b, -h, T(1))
425-
__apply_algebraic_constraint!(residᵢ, algebraic_indices, f!, yᵢ₊₁, p, mesh[i + 1], tmp1)
424+
__apply_algebraic_constraint_oop!(residᵢ, algebraic_indices, f!, yᵢ₊₁, p, mesh[i + 1])
426425
end
427426
return residuals
428427
end
@@ -434,7 +433,6 @@ end
434433
)
435434
(; b) = TU
436435
(; nest_prob, alg) = cache
437-
tmp1 = similar(y[1])
438436

439437
residuals = [safe_similar(yᵢ) for yᵢ in y[1:(end - 1)]]
440438

@@ -458,7 +456,7 @@ end
458456

459457
@. residᵢ = yᵢ₊₁ - yᵢ
460458
__maybe_matmul!(residᵢ, nestsol.u, b, -h, T(1))
461-
__apply_algebraic_constraint!(residᵢ, algebraic_indices, f!, yᵢ₊₁, p, mesh[i + 1], tmp1)
459+
__apply_algebraic_constraint_oop!(residᵢ, algebraic_indices, f!, yᵢ₊₁, p, mesh[i + 1])
462460
end
463461
return residuals
464462
end

‎lib/BoundaryValueDiffEqFIRK/src/firk.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -278,6 +278,7 @@ function init_nested(
278278
nestprob_p = zeros(T, M + 2)
279279

280280
algebraic_indices = __get_algebraic_indices(prob.f.mass_matrix)
281+
__check_dae_adaptivity(algebraic_indices, adaptive)
281282

282283
mm = prob.f.mass_matrix
283284

@@ -453,7 +454,7 @@ function init_expanded(
453454
end
454455

455456
algebraic_indices = __get_algebraic_indices(prob.f.mass_matrix)
456-
457+
__check_dae_adaptivity(algebraic_indices, adaptive)
457458

458459
return FIRKCacheExpand{iip, T, typeof(diffcache), tune_parameters}(
459460
alg_order(alg), stage, M, size(u0), f, prob.f.mass_matrix, algebraic_indices, bc, prob_, prob.problem_type, prob.p,

‎lib/BoundaryValueDiffEqFIRK/test/expanded/dae_tests.jl‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
using BoundaryValueDiffEqFIRK
2+
using Test
3+
4+
# Index-1 DAE BVPs solved via unprojected collocation (Ascher & Spiteri 1994).
5+
# The algebraic constraint is enforced exactly at the mesh points, so accuracy for
6+
# algebraic variables is only checked there: the continuous interpolant is not
7+
# accurate for algebraic components, which is also why mesh adaptivity is not
8+
# supported for DAEs (the defect estimate cannot converge).
9+
#
10+
# LobattoIIIa and LobattoIIIb are excluded: their tableau structure leaves the
11+
# algebraic components of the stages underdetermined, so they cannot solve DAEs.
12+
13+
@testset "Simple index-1 DAE" begin
14+
using BoundaryValueDiffEqFIRK, SciMLBase
15+
using LinearAlgebra
16+
17+
nested = false
18+
19+
# u1' = u2, 0 = u2 - cos(t) with u1(0) = 0
20+
# Analytic solution: u1 = sin(t), u2 = cos(t)
21+
function f1!(du, u, p, t)
22+
du[1] = u[2]
23+
du[2] = u[2] - cos(t)
24+
end
25+
f1(u, p, t) = [u[2], u[2] - cos(t)]
26+
function bc1!(res, u, p, t)
27+
res[1] = u(0.0)[1]
28+
res[2] = u(0.0)[2] - 1.0
29+
end
30+
bc1(u, p, t) = [u(0.0)[1], u(0.0)[2] - 1.0]
31+
32+
mass_matrix = [1.0 0.0; 0.0 0.0]
33+
tspan = (0.0, pi / 2)
34+
prob_iip = BVProblem(BVPFunction(f1!, bc1!; mass_matrix), [0.0, 1.0], tspan)
35+
prob_oop = BVProblem(BVPFunction(f1, bc1; mass_matrix), [0.0, 1.0], tspan)
36+
37+
@testset "$(nameof(typeof(alg))), $(SciMLBase.isinplace(prob) ? "iip" : "oop")" for alg in
38+
(
39+
RadauIIa3(; nested_nlsolve = nested), RadauIIa5(; nested_nlsolve = nested),
40+
RadauIIa7(; nested_nlsolve = nested), LobattoIIIc4(; nested_nlsolve = nested),
41+
),
42+
prob in (prob_iip, prob_oop)
43+
44+
sol = solve(prob, alg; dt = 0.01, adaptive = false)
45+
@test SciMLBase.successful_retcode(sol)
46+
@test maximum(abs(sol.u[i][1] - sin(sol.t[i])) for i in eachindex(sol.t)) < 1.0e-10
47+
@test maximum(abs(sol.u[i][2] - cos(sol.t[i])) for i in eachindex(sol.t)) < 1.0e-12
48+
end
49+
end
50+
51+
@testset "Ascher & Spiteri example problem 1" begin
52+
using BoundaryValueDiffEqFIRK, SciMLBase
53+
using LinearAlgebra
54+
55+
nested = false
56+
57+
# Singular index-1 BVDAE from the Ascher & Spiteri paper.
58+
# Analytic solution: [sin(t), sin(t), 1, 0]
59+
function f2!(du, u, p, t)
60+
e = 2.7
61+
du[1] = (1 + u[2] - sin(t)) * u[4] + cos(t)
62+
du[2] = cos(t)
63+
du[3] = u[4]
64+
du[4] = (u[1] - sin(t)) * (u[4] - e^t)
65+
end
66+
function f2(u, p, t)
67+
e = 2.7
68+
return [
69+
(1 + u[2] - sin(t)) * u[4] + cos(t), cos(t),
70+
u[4], (u[1] - sin(t)) * (u[4] - e^t),
71+
]
72+
end
73+
function bc2!(res, u, p, t)
74+
res[1] = u(0.0)[1]
75+
res[2] = u(0.0)[3] - 1.0
76+
res[3] = u(1.0)[2] - sin(1.0)
77+
res[4] = u(0.0)[4]
78+
end
79+
bc2(u, p, t) = [u(0.0)[1], u(0.0)[3] - 1.0, u(1.0)[2] - sin(1.0), u(0.0)[4]]
80+
f2_analytic(t) = [sin(t), sin(t), 1.0, 0.0]
81+
82+
mass_matrix = [
83+
1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0;
84+
0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0
85+
]
86+
tspan = (0.0, 1.0)
87+
prob_iip = BVProblem(BVPFunction(f2!, bc2!; mass_matrix), zeros(4), tspan)
88+
prob_oop = BVProblem(BVPFunction(f2, bc2; mass_matrix), zeros(4), tspan)
89+
90+
@testset "RadauIIa5, $(SciMLBase.isinplace(prob) ? "iip" : "oop")" for prob in
91+
(prob_iip, prob_oop)
92+
93+
sol = solve(prob, RadauIIa5(; nested_nlsolve = nested); dt = 0.01, adaptive = false)
94+
@test SciMLBase.successful_retcode(sol)
95+
err = maximum(
96+
maximum(abs.(sol.u[i] .- f2_analytic(sol.t[i]))) for i in eachindex(sol.t)
97+
)
98+
@test err < 1.0e-8
99+
end
100+
end
101+
102+
@testset "Mesh adaptivity is not supported for DAEs" begin
103+
using BoundaryValueDiffEqFIRK, SciMLBase
104+
using LinearAlgebra
105+
106+
nested = false
107+
108+
function f3!(du, u, p, t)
109+
du[1] = u[2]
110+
du[2] = u[2] - cos(t)
111+
end
112+
function bc3!(res, u, p, t)
113+
res[1] = u(0.0)[1]
114+
res[2] = u(0.0)[2] - 1.0
115+
end
116+
mass_matrix = [1.0 0.0; 0.0 0.0]
117+
prob = BVProblem(BVPFunction(f3!, bc3!; mass_matrix), [0.0, 1.0], (0.0, pi / 2))
118+
119+
@test_throws ArgumentError solve(prob, RadauIIa5(; nested_nlsolve = nested); dt = 0.05)
120+
@test_throws ArgumentError solve(
121+
prob, RadauIIa5(; nested_nlsolve = nested); dt = 0.05, adaptive = true
122+
)
123+
end

‎lib/BoundaryValueDiffEqFIRK/test/nested/dae_tests.jl‎

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
using BoundaryValueDiffEqFIRK
2+
using Test
3+
4+
# Index-1 DAE BVPs solved via unprojected collocation (Ascher & Spiteri 1994).
5+
# The algebraic constraint is enforced exactly at the mesh points, so accuracy for
6+
# algebraic variables is only checked there: the continuous interpolant is not
7+
# accurate for algebraic components, which is also why mesh adaptivity is not
8+
# supported for DAEs (the defect estimate cannot converge).
9+
#
10+
# LobattoIIIa and LobattoIIIb are excluded: their tableau structure leaves the
11+
# algebraic components of the stages underdetermined, so they cannot solve DAEs.
12+
# Problems whose constraint Jacobian is singular on the solution (e.g. the Ascher
13+
# & Spiteri example problem 1) are also excluded: the nested nonlinear solve of
14+
# the stage equations hits the singular constraint directly.
15+
16+
@testset "Simple index-1 DAE" begin
17+
using BoundaryValueDiffEqFIRK, SciMLBase
18+
using LinearAlgebra
19+
20+
nested = true
21+
22+
# u1' = u2, 0 = u2 - cos(t) with u1(0) = 0
23+
# Analytic solution: u1 = sin(t), u2 = cos(t)
24+
function f1!(du, u, p, t)
25+
du[1] = u[2]
26+
du[2] = u[2] - cos(t)
27+
end
28+
f1(u, p, t) = [u[2], u[2] - cos(t)]
29+
function bc1!(res, u, p, t)
30+
res[1] = u(0.0)[1]
31+
res[2] = u(0.0)[2] - 1.0
32+
end
33+
bc1(u, p, t) = [u(0.0)[1], u(0.0)[2] - 1.0]
34+
35+
mass_matrix = [1.0 0.0; 0.0 0.0]
36+
tspan = (0.0, pi / 2)
37+
prob_iip = BVProblem(BVPFunction(f1!, bc1!; mass_matrix), [0.0, 1.0], tspan)
38+
prob_oop = BVProblem(BVPFunction(f1, bc1; mass_matrix), [0.0, 1.0], tspan)
39+
40+
@testset "$(nameof(typeof(alg))), $(SciMLBase.isinplace(prob) ? "iip" : "oop")" for alg in
41+
(
42+
RadauIIa3(; nested_nlsolve = nested), RadauIIa5(; nested_nlsolve = nested),
43+
LobattoIIIc4(; nested_nlsolve = nested),
44+
),
45+
prob in (prob_iip, prob_oop)
46+
47+
sol = solve(prob, alg; dt = 0.01, adaptive = false)
48+
@test SciMLBase.successful_retcode(sol)
49+
@test maximum(abs(sol.u[i][1] - sin(sol.t[i])) for i in eachindex(sol.t)) < 1.0e-10
50+
@test maximum(abs(sol.u[i][2] - cos(sol.t[i])) for i in eachindex(sol.t)) < 1.0e-12
51+
end
52+
end
53+
54+
@testset "Mesh adaptivity is not supported for DAEs" begin
55+
using BoundaryValueDiffEqFIRK, SciMLBase
56+
using LinearAlgebra
57+
58+
nested = true
59+
60+
function f3!(du, u, p, t)
61+
du[1] = u[2]
62+
du[2] = u[2] - cos(t)
63+
end
64+
function bc3!(res, u, p, t)
65+
res[1] = u(0.0)[1]
66+
res[2] = u(0.0)[2] - 1.0
67+
end
68+
mass_matrix = [1.0 0.0; 0.0 0.0]
69+
prob = BVProblem(BVPFunction(f3!, bc3!; mass_matrix), [0.0, 1.0], (0.0, pi / 2))
70+
71+
@test_throws ArgumentError solve(prob, RadauIIa5(; nested_nlsolve = nested); dt = 0.05)
72+
@test_throws ArgumentError solve(
73+
prob, RadauIIa5(; nested_nlsolve = nested); dt = 0.05, adaptive = true
74+
)
75+
end

‎lib/BoundaryValueDiffEqFIRK/test/runtests.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,12 +18,14 @@ run_tests(;
1818
@time @safetestset "FIRK Expanded AD Tests" include("expanded/ad_tests.jl")
1919
@time @safetestset "FIRK Expanded Ensemble Tests" include("expanded/ensemble_tests.jl")
2020
@time @safetestset "FIRK Expanded Singular BVP Tests" include("expanded/singular_bvp_tests.jl")
21+
@time @safetestset "FIRK Expanded DAE Tests" include("expanded/dae_tests.jl")
2122
return @time @safetestset "FIRK Expanded VectorOfVector Initials Tests" include("expanded/vectorofvector_initials_tests.jl")
2223
end,
2324
"NESTED" => function ()
2425
@time @safetestset "FIRK Nested Basic Tests" include("nested/firk_basic_tests.jl")
2526
@time @safetestset "FIRK Nested NLLS Tests" include("nested/nlls_tests.jl")
2627
@time @safetestset "FIRK Nested Ensemble Tests" include("nested/ensemble_tests.jl")
28+
@time @safetestset "FIRK Nested DAE Tests" include("nested/dae_tests.jl")
2729
return @time @safetestset "FIRK Nested VectorOfVector Initials Tests" include("nested/vectorofvector_initials_tests.jl")
2830
end,
2931
),

‎lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
2727
__default_sparsity_detector, __build_cost, __add_singular_term!,
2828
__apply_mass_matrix!,
2929
__get_algebraic_indices, __apply_algebraic_constraint!,
30-
__apply_algebraic_constraint_oop!
30+
__apply_algebraic_constraint_oop!, __check_dae_adaptivity
3131

3232
using ConcreteStructs: @concrete
3333
using DifferentiationInterface: DifferentiationInterface, Constant, prepare_jacobian

‎lib/BoundaryValueDiffEqMIRK/src/mirk.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -249,6 +249,7 @@ function SciMLBase.__init(
249249
end
250250

251251
algebraic_indices = __get_algebraic_indices(prob.f.mass_matrix)
252+
__check_dae_adaptivity(algebraic_indices, adaptive)
252253

253254
return MIRKCache{iip, T, use_both, typeof(diffcache), tune_parameters}(
254255
alg_order(alg), stage, N, size(u0), f, prob.f.mass_matrix, algebraic_indices, bc, prob_, prob.problem_type, prob.p, alg,

0 commit comments

Comments
 (0)