Skip to content

Commit 37e9e11

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 ffd62a6 commit 37e9e11

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
@@ -973,25 +973,46 @@ end
973973
end
974974

975975
@inline function __apply_algebraic_constraint!(
976-
residᵢ, ::Nothing, f!, yᵢ₊₁, p, t, tmp)
976+
residᵢ, ::Nothing, f!, yᵢ₊₁, p, t, tmp
977+
)
977978
return nothing
978979
end
979980

980981
@inline function __apply_algebraic_constraint!(
981-
residᵢ, algebraic_indices::Vector{Int}, f!, yᵢ₊₁, p, t, tmp)
982+
residᵢ, algebraic_indices::Vector{Int}, f!, yᵢ₊₁, p, t, tmp
983+
)
982984
f!(tmp, yᵢ₊₁, p, t)
983985
residᵢ[algebraic_indices] .= tmp[algebraic_indices]
984986
return nothing
985987
end
986988

987989
@inline function __apply_algebraic_constraint_oop!(
988-
residᵢ, ::Nothing, f, yᵢ₊₁, p, t)
990+
residᵢ, ::Nothing, f, yᵢ₊₁, p, t
991+
)
989992
return nothing
990993
end
991994

992995
@inline function __apply_algebraic_constraint_oop!(
993-
residᵢ, algebraic_indices::Vector{Int}, f, yᵢ₊₁, p, t)
996+
residᵢ, algebraic_indices::Vector{Int}, f, yᵢ₊₁, p, t
997+
)
994998
tmp = f(yᵢ₊₁, p, t)
995999
residᵢ[algebraic_indices] .= tmp[algebraic_indices]
9961000
return nothing
9971001
end
1002+
1003+
@inline __check_dae_adaptivity(::Nothing, adaptive::Bool) = nothing
1004+
1005+
@inline function __check_dae_adaptivity(::Vector{Int}, adaptive::Bool)
1006+
if adaptive
1007+
throw(
1008+
ArgumentError(
1009+
"Adaptive mesh refinement is not supported for DAE problems (mass " *
1010+
"matrices with zero rows): the collocation interpolant is inaccurate " *
1011+
"for algebraic variables, so the defect estimate cannot converge. " *
1012+
"Pass `adaptive = false`, or use a solver from " *
1013+
"BoundaryValueDiffEqAscher.jl, which supports mesh adaptivity for DAEs."
1014+
)
1015+
)
1016+
end
1017+
return nothing
1018+
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: BoundaryValueDiffEqCore,
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
@@ -44,6 +44,7 @@ run_tests(;
4444
@time @safetestset "FIRK Expanded NLLS Tests" include("expanded/nlls_tests.jl")
4545
@time @safetestset "FIRK Expanded Ensemble Tests" include("expanded/ensemble_tests.jl")
4646
@time @safetestset "FIRK Expanded Singular BVP Tests" include("expanded/singular_bvp_tests.jl")
47+
@time @safetestset "FIRK Expanded DAE Tests" include("expanded/dae_tests.jl")
4748
return @time @safetestset "FIRK Expanded VectorOfVector Initials Tests" include("expanded/vectorofvector_initials_tests.jl")
4849
end,
4950
# The nested formulation is split for the same reason as the expanded one, but its
@@ -85,6 +86,7 @@ run_tests(;
8586
@time @safetestset "FIRK Nested NLLS Tests" include("nested/nlls_tests.jl")
8687
@time @safetestset "FIRK Nested Underconstrained NLLS Tests" include("nested/nlls_underconstrained_tests.jl")
8788
@time @safetestset "FIRK Nested Ensemble Tests" include("nested/ensemble_tests.jl")
89+
@time @safetestset "FIRK Nested DAE Tests" include("nested/dae_tests.jl")
8890
return @time @safetestset "FIRK Nested VectorOfVector Initials Tests" include("nested/vectorofvector_initials_tests.jl")
8991
end,
9092
# AD: the different-AD-backend compatibility tests. Enzyme and Mooncake are

‎lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ using BoundaryValueDiffEqCore: BoundaryValueDiffEqCore,
2929
__default_sparsity_detector, __build_cost, __add_singular_term!,
3030
__apply_mass_matrix!,
3131
__get_algebraic_indices, __apply_algebraic_constraint!,
32-
__apply_algebraic_constraint_oop!
32+
__apply_algebraic_constraint_oop!, __check_dae_adaptivity
3333

3434
using ConcreteStructs: @concrete
3535
using DifferentiationInterface: DifferentiationInterface, Constant

‎lib/BoundaryValueDiffEqMIRK/src/mirk.jl‎

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

260260
algebraic_indices = __get_algebraic_indices(prob.f.mass_matrix)
261+
__check_dae_adaptivity(algebraic_indices, adaptive)
261262

262263
return MIRKCache{iip, T, use_both, typeof(diffcache), tune_parameters}(
263264
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)