Skip to content

Commit 86d4144

Browse files
vyuduvyudu
committed
test more solvers:
1 parent e74e047 commit 86d4144

File tree

2 files changed

+172
-163
lines changed

2 files changed

+172
-163
lines changed

‎src/systems/diffeqs/abstractodesystem.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -969,7 +969,7 @@ function validate_constraint_syms(eq, constraintsts, constraintps, sts, ps, iv)
969969
for var in constraintsts
970970
if length(arguments(var)) > 1
971971
error("Too many arguments for variable $var.")
972-
elseif arguments(var) == iv
972+
elseif isequal(arguments(var)[1], iv)
973973
var ∈ sts || error("Constraint equation $eq contains a variable $var that is not a variable of the ODESystem.")
974974
error("Constraint equation $eq contains a variable $var that does not have a specified argument. Such equations should be specified as algebraic equations to the ODESystem rather than a boundary constraints.")
975975
else

‎test/bvproblem.jl‎

Lines changed: 171 additions & 162 deletions
Original file line numberDiff line numberDiff line change
@@ -1,183 +1,186 @@
11
### TODO: update when BoundaryValueDiffEqAscher is updated to use the normal boundary condition conventions
22

33
using BoundaryValueDiffEq, OrdinaryDiffEq, BoundaryValueDiffEqAscher
4+
using BenchmarkTools
45
using ModelingToolkit
56
using SciMLBase
67
using ModelingToolkit: t_nounits as t, D_nounits as D
78
import ModelingToolkit: process_constraints
89

910
### Test Collocation solvers on simple problems
10-
solvers = [MIRK4, RadauIIa5]
11+
solvers = [MIRK4, RadauIIa5, LobattoIIIa3]
1112
daesolvers = [Ascher2, Ascher4, Ascher6]
1213

13-
let
14-
@parameters α=7.5 β=4.0 γ=8.0 δ=5.0
15-
@variables x(t)=1.0 y(t)=2.0
16-
17-
eqs = [D(x) ~ α * x - β * x * y,
18-
D(y) ~ -γ * y + δ * x * y]
19-
20-
u0map = [x => 1.0, y => 2.0]
21-
parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0]
22-
tspan = (0.0, 10.0)
23-
24-
@mtkbuild lotkavolterra = ODESystem(eqs, t)
25-
op = ODEProblem(lotkavolterra, u0map, tspan, parammap)
26-
osol = solve(op, Vern9())
27-
28-
bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(
29-
lotkavolterra, u0map, tspan, parammap; eval_expression = true)
30-
31-
for solver in solvers
32-
sol = solve(bvp, solver(), dt = 0.01)
33-
@test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
34-
@test sol.u[1] == [1.0, 2.0]
35-
end
36-
37-
# Test out of place
38-
bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(
39-
lotkavolterra, u0map, tspan, parammap; eval_expression = true)
40-
41-
for solver in solvers
42-
sol = solve(bvp2, solver(), dt = 0.01)
43-
@test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
44-
@test sol.u[1] == [1.0, 2.0]
45-
end
46-
end
47-
48-
### Testing on pendulum
49-
let
50-
@parameters g=9.81 L=1.0
51-
@variables θ(t) = π / 2 θ_t(t)
52-
53-
eqs = [D(θ) ~ θ_t
54-
D(θ_t) ~ -(g / L) * sin(θ)]
55-
56-
@mtkbuild pend = ODESystem(eqs, t)
57-
58-
u0map = [θ => π / 2, θ_t => π / 2]
59-
parammap = [:L => 1.0, :g => 9.81]
60-
tspan = (0.0, 6.0)
61-
62-
op = ODEProblem(pend, u0map, tspan, parammap)
63-
osol = solve(op, Vern9())
64-
65-
bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap)
66-
for solver in solvers
67-
sol = solve(bvp, solver(), dt = 0.01)
68-
@test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
69-
@test sol.u[1] == [π / 2, π / 2]
70-
end
71-
72-
# Test out-of-place
73-
bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap)
74-
75-
for solver in solvers
76-
sol = solve(bvp2, solver(), dt = 0.01)
77-
@test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
78-
@test sol.u[1] == [π / 2, π / 2]
79-
end
80-
end
81-
82-
##################################################################
83-
### ODESystem with constraint equations, DAEs with constraints ###
84-
##################################################################
85-
86-
# Test generation of boundary condition function using `process_constraints`. Compare solutions to manually written boundary conditions
87-
let
88-
@parameters α=1.5 β=1.0 γ=3.0 δ=1.0
89-
@variables x(..) y(..)
90-
eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t),
91-
D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)]
92-
93-
tspan = (0., 1.)
94-
@mtkbuild lksys = ODESystem(eqs, t)
95-
96-
function lotkavolterra!(du, u, p, t)
97-
du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
98-
du[2] = -p[4]*u[2] + p[3]*u[1]*u[2]
99-
end
100-
101-
function lotkavolterra(u, p, t)
102-
[p[1]*u[1] - p[2]*u[1]*u[2], -p[4]*u[2] + p[3]*u[1]*u[2]]
103-
end
104-
# Compare the built bc function to the actual constructed one.
105-
function bc!(resid, u, p, t)
106-
resid[1] = u[1][1] - 1.
107-
resid[2] = u[1][2] - 2.
108-
nothing
109-
end
110-
function bc(u, p, t)
111-
[u[1][1] - 1., u[1][2] - 2.]
112-
end
113-
114-
u0 = [1., 2.]; p = [1.5, 1., 1., 3.]
115-
genbc_iip = ModelingToolkit.process_constraints(lksys, nothing, u0, [1, 2], tspan, true)
116-
genbc_oop = ModelingToolkit.process_constraints(lksys, nothing, u0, [1, 2], tspan, false)
117-
118-
bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1.,2.], tspan, p)
119-
bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1.,2.], tspan, p)
120-
121-
sol1 = solve(bvpi1, MIRK4(), dt = 0.05)
122-
sol2 = solve(bvpi2, MIRK4(), dt = 0.05)
123-
@test sol1 ≈ sol2
124-
125-
bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p)
126-
bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p)
127-
128-
sol1 = solve(bvpo1, MIRK4(), dt = 0.05)
129-
sol2 = solve(bvpo2, MIRK4(), dt = 0.05)
130-
@test sol1 ≈ sol2
131-
132-
# Test with a constraint.
133-
constraints = [y(0.5) ~ 2.]
134-
135-
function bc!(resid, u, p, t)
136-
resid[1] = u(0.0)[1] - 1.
137-
resid[2] = u(0.5)[2] - 2.
138-
end
139-
function bc(u, p, t)
140-
[u(0.0)[1] - 1., u(0.5)[2] - 2.]
141-
end
142-
143-
u0 = [1, 1.]
144-
genbc_iip = ModelingToolkit.process_constraints(lksys, constraints, u0, [1], tspan, true)
145-
genbc_oop = ModelingToolkit.process_constraints(lksys, constraints, u0, [1], tspan, false)
146-
147-
bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, u0, tspan, p)
148-
bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, u0, tspan, p)
149-
bvpi3 = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
150-
bvpi4 = SciMLBase.BVProblem{true, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
151-
152-
sol1 = @btime solve($bvpi1, MIRK4(), dt = 0.01)
153-
sol2 = @btime solve($bvpi2, MIRK4(), dt = 0.01)
154-
sol3 = @btime solve($bvpi3, MIRK4(), dt = 0.01)
155-
sol4 = @btime solve($bvpi4, MIRK4(), dt = 0.01)
156-
@test sol1 ≈ sol2 ≈ sol3 ≈ sol4 # don't get true equality here, not sure why
157-
158-
bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p)
159-
bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p)
160-
bvpo3 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
161-
162-
sol1 = @btime solve($bvpo1, MIRK4(), dt = 0.05)
163-
sol2 = @btime solve($bvpo2, MIRK4(), dt = 0.05)
164-
sol3 = @btime solve($bvpo3, MIRK4(), dt = 0.05)
165-
@test sol1 ≈ sol2 ≈ sol3
166-
end
14+
# let
15+
# @parameters α=7.5 β=4.0 γ=8.0 δ=5.0
16+
# @variables x(t)=1.0 y(t)=2.0
17+
#
18+
# eqs = [D(x) ~ α * x - β * x * y,
19+
# D(y) ~ -γ * y + δ * x * y]
20+
#
21+
# u0map = [x => 1.0, y => 2.0]
22+
# parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0]
23+
# tspan = (0.0, 10.0)
24+
#
25+
# @mtkbuild lotkavolterra = ODESystem(eqs, t)
26+
# op = ODEProblem(lotkavolterra, u0map, tspan, parammap)
27+
# osol = solve(op, Vern9())
28+
#
29+
# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(
30+
# lotkavolterra, u0map, tspan, parammap; eval_expression = true)
31+
#
32+
# for solver in solvers
33+
# sol = solve(bvp, solver(), dt = 0.01)
34+
# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
35+
# @test sol.u[1] == [1.0, 2.0]
36+
# end
37+
#
38+
# # Test out of place
39+
# bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(
40+
# lotkavolterra, u0map, tspan, parammap; eval_expression = true)
41+
#
42+
# for solver in solvers
43+
# sol = solve(bvp2, solver(), dt = 0.01)
44+
# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
45+
# @test sol.u[1] == [1.0, 2.0]
46+
# end
47+
# end
48+
#
49+
# ### Testing on pendulum
50+
# let
51+
# @parameters g=9.81 L=1.0
52+
# @variables θ(t) = π / 2 θ_t(t)
53+
#
54+
# eqs = [D(θ) ~ θ_t
55+
# D(θ_t) ~ -(g / L) * sin(θ)]
56+
#
57+
# @mtkbuild pend = ODESystem(eqs, t)
58+
#
59+
# u0map = [θ => π / 2, θ_t => π / 2]
60+
# parammap = [:L => 1.0, :g => 9.81]
61+
# tspan = (0.0, 6.0)
62+
#
63+
# op = ODEProblem(pend, u0map, tspan, parammap)
64+
# osol = solve(op, Vern9())
65+
#
66+
# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap)
67+
# for solver in solvers
68+
# sol = solve(bvp, solver(), dt = 0.01)
69+
# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
70+
# @test sol.u[1] == [π / 2, π / 2]
71+
# end
72+
#
73+
# # Test out-of-place
74+
# bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap)
75+
#
76+
# for solver in solvers
77+
# sol = solve(bvp2, solver(), dt = 0.01)
78+
# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
79+
# @test sol.u[1] == [π / 2, π / 2]
80+
# end
81+
# end
82+
#
83+
# ##################################################################
84+
# ### ODESystem with constraint equations, DAEs with constraints ###
85+
# ##################################################################
86+
#
87+
# # Test generation of boundary condition function using `process_constraints`. Compare solutions to manually written boundary conditions
88+
# let
89+
# @parameters α=1.5 β=1.0 γ=3.0 δ=1.0
90+
# @variables x(..) y(..)
91+
# eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t),
92+
# D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)]
93+
#
94+
# tspan = (0., 1.)
95+
# @mtkbuild lksys = ODESystem(eqs, t)
96+
#
97+
# function lotkavolterra!(du, u, p, t)
98+
# du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
99+
# du[2] = -p[4]*u[2] + p[3]*u[1]*u[2]
100+
# end
101+
#
102+
# function lotkavolterra(u, p, t)
103+
# [p[1]*u[1] - p[2]*u[1]*u[2], -p[4]*u[2] + p[3]*u[1]*u[2]]
104+
# end
105+
# # Compare the built bc function to the actual constructed one.
106+
# function bc!(resid, u, p, t)
107+
# resid[1] = u[1][1] - 1.
108+
# resid[2] = u[1][2] - 2.
109+
# nothing
110+
# end
111+
# function bc(u, p, t)
112+
# [u[1][1] - 1., u[1][2] - 2.]
113+
# end
114+
#
115+
# u0 = [1., 2.]; p = [1.5, 1., 1., 3.]
116+
# genbc_iip = ModelingToolkit.process_constraints(lksys, nothing, u0, [1, 2], tspan, true)
117+
# genbc_oop = ModelingToolkit.process_constraints(lksys, nothing, u0, [1, 2], tspan, false)
118+
#
119+
# bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1.,2.], tspan, p)
120+
# bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1.,2.], tspan, p)
121+
#
122+
# sol1 = solve(bvpi1, MIRK4(), dt = 0.05)
123+
# sol2 = solve(bvpi2, MIRK4(), dt = 0.05)
124+
# @test sol1 ≈ sol2
125+
#
126+
# bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p)
127+
# bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p)
128+
#
129+
# sol1 = solve(bvpo1, MIRK4(), dt = 0.05)
130+
# sol2 = solve(bvpo2, MIRK4(), dt = 0.05)
131+
# @test sol1 ≈ sol2
132+
#
133+
# # Test with a constraint.
134+
# constraints = [y(0.5) ~ 2.]
135+
#
136+
# function bc!(resid, u, p, t)
137+
# resid[1] = u(0.0)[1] - 1.
138+
# resid[2] = u(0.5)[2] - 2.
139+
# end
140+
# function bc(u, p, t)
141+
# [u(0.0)[1] - 1., u(0.5)[2] - 2.]
142+
# end
143+
#
144+
# u0 = [1, 1.]
145+
# genbc_iip = ModelingToolkit.process_constraints(lksys, constraints, u0, [1], tspan, true)
146+
# genbc_oop = ModelingToolkit.process_constraints(lksys, constraints, u0, [1], tspan, false)
147+
#
148+
# bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, u0, tspan, p)
149+
# bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, u0, tspan, p)
150+
# bvpi3 = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
151+
# bvpi4 = SciMLBase.BVProblem{true, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
152+
#
153+
# sol1 = @btime solve($bvpi1, MIRK4(), dt = 0.01)
154+
# sol2 = @btime solve($bvpi2, MIRK4(), dt = 0.01)
155+
# sol3 = @btime solve($bvpi3, MIRK4(), dt = 0.01)
156+
# sol4 = @btime solve($bvpi4, MIRK4(), dt = 0.01)
157+
# @test sol1 ≈ sol2 ≈ sol3 ≈ sol4 # don't get true equality here, not sure why
158+
#
159+
# bvpo1 = BVProblem(lotkavolterra, bc, u0, tspan, p)
160+
# bvpo2 = BVProblem(lotkavolterra, genbc_oop, u0, tspan, p)
161+
# bvpo3 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
162+
#
163+
# sol1 = @btime solve($bvpo1, MIRK4(), dt = 0.05)
164+
# sol2 = @btime solve($bvpo2, MIRK4(), dt = 0.05)
165+
# sol3 = @btime solve($bvpo3, MIRK4(), dt = 0.05)
166+
# @test sol1 ≈ sol2 ≈ sol3
167+
# end
167168

168-
function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.05, atol = 1e-4)
169+
function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.05, atol = 1e-3)
169170
for solver in solvers
170171
println("Solver: $solver")
171172
sol = @btime solve($prob, $solver(), dt = $dt, abstol = $atol)
172173
@test SciMLBase.successful_retcode(sol.retcode)
173174
p = prob.p; t = sol.t; bc = prob.f.bc
174175
ns = length(prob.u0)
175-
if isinplace(bvp.f)
176+
if isinplace(prob.f)
176177
resid = zeros(ns)
177178
bc(resid, sol, p, t)
178179
@test isapprox(zeros(ns), resid; atol)
180+
@show resid
179181
else
180182
@test isapprox(zeros(ns), bc(sol, p, t); atol)
183+
@show bc(sol, p, t)
181184
end
182185

183186
for (k, v) in u0map
@@ -194,6 +197,12 @@ function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.
194197
end
195198
end
196199

200+
solvers = [RadauIIa3, RadauIIa5, RadauIIa7,
201+
LobattoIIIa2, LobattoIIIa4, LobattoIIIa5,
202+
LobattoIIIb2, LobattoIIIb3, LobattoIIIb4, LobattoIIIb5,
203+
LobattoIIIc2, LobattoIIIc3, LobattoIIIc4, LobattoIIIc5]
204+
weird = [MIRK2, MIRK5, RadauIIa2]
205+
daesolvers = []
197206
# Simple ODESystem with BVP constraints.
198207
let
199208
@parameters α=1.5 β=1.0 γ=3.0 δ=1.0
@@ -213,8 +222,8 @@ let
213222
test_solvers(solvers, bvp, u0map, constraints; dt = 0.05)
214223

215224
# Testing that more complicated constraints give correct solutions.
216-
constraints = [y(.2) + x(.8) ~ 3., y(.3) + x(.5) ~ 5.]
217-
bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, u0map, tspan; guesses, constraints, jac = true)
225+
constraints = [y(.2) + x(.8) ~ 3., y(.3) ~ 2.]
226+
bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, u0map, tspan; guesses, constraints)
218227
test_solvers(solvers, bvp, u0map, constraints; dt = 0.05)
219228

220229
constraints = [α * β - x(.6) ~ 0.0, y(.2) ~ 3.]
@@ -224,14 +233,14 @@ let
224233
# Testing that errors are properly thrown when malformed constraints are given.
225234
@variables bad(..)
226235
constraints = [x(1.) + bad(3.) ~ 10]
227-
@test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses, constraints)
236+
@test_throws ErrorException bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses, constraints)
228237

229238
constraints = [x(t) + y(t) ~ 3]
230-
@test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses, constraints)
239+
@test_throws ErrorException bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses, constraints)
231240

232241
@parameters bad2
233242
constraints = [bad2 + x(0.) ~ 3]
234-
@test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses, constraints)
243+
@test_throws ErrorException bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses, constraints)
235244
end
236245

237246
# Cartesian pendulum from the docs.

0 commit comments

Comments
 (0)