Skip to content

Commit bbda8b1

Browse files
AayushSabharwalAayushSabharwal
committed
test: test SemilinearODEProblem, SemilinearODEFunction 
1 parent d3a5b70 commit bbda8b1

File tree

2 files changed

+212
-0
lines changed

2 files changed

+212
-0
lines changed

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -100,6 +100,7 @@ end
100100
@safetestset "Subsystem replacement" include("substitute_component.jl")
101101
@safetestset "Linearization Tests" include("linearize.jl")
102102
@safetestset "LinearProblem Tests" include("linearproblem.jl")
103+
@safetestset "SemilinearODEProblem tests" include("semilinearodeproblem.jl")
103104
end
104105
end
105106

‎test/semilinearodeproblem.jl

Lines changed: 211 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,211 @@
1+
using ModelingToolkit
2+
using OrdinaryDiffEq
3+
using LinearAlgebra
4+
using Test
5+
using ModelingToolkit: t_nounits as t, D_nounits as D
6+
7+
# from https://docs.sciml.ai/SciMLBenchmarksOutput/dev/AstroChem/nelson/
8+
@testset "Astrochem model" begin
9+
function Nelson!(du, u, p, t)
10+
T, Av, Go, n_H, shield = p
11+
12+
# 1: H2
13+
du[1] = -1.2e-17 * u[1] +
14+
n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
15+
n_H * 4e-16 * u[1] * u[12] -
16+
n_H * 7e-15 * u[1] * u[5] +
17+
n_H * 1.7e-9 * u[10] * u[2] +
18+
n_H * 2e-9 * u[2] * u[6] +
19+
n_H * 2e-9 * u[2] * u[14] +
20+
n_H * 8e-10 * u[2] * u[8] + sin(u[1]) / 1e20 # artificial nonlinear term for testing
21+
22+
# 2: H3+
23+
du[2] = 1.2e-17 * u[1] +
24+
n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
25+
n_H * 1.7e-9 * u[10] * u[2] -
26+
n_H * 2e-9 * u[2] * u[6] -
27+
n_H * 2e-9 * u[2] * u[14] -
28+
n_H * 8e-10 * u[2] * u[8]
29+
30+
# 3: e
31+
du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
32+
n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
33+
n_H * (3.3e-5 * u[11] * u[3]) / T +
34+
1.2e-17 * u[1] -
35+
n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
36+
6.8e-18 * u[4] -
37+
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
38+
3e-10 * Go * exp(-3 * Av) * u[6] +
39+
n_H * 2e-9 * u[2] * u[13]
40+
+2.0e-10 * Go * exp(-1.9 * Av) * u[14]
41+
42+
# 4: He
43+
du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
44+
6.8e-18 * u[4] +
45+
n_H * 7e-15 * u[1] * u[5] +
46+
n_H * 1.6e-9 * u[10] * u[5]
47+
48+
# 5: He+
49+
du[5] = 6.8e-18 * u[4] -
50+
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
51+
n_H * 7e-15 * u[1] * u[5] -
52+
n_H * 1.6e-9 * u[10] * u[5]
53+
54+
# 6: C
55+
du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
56+
n_H * 2e-9 * u[2] * u[6] -
57+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
58+
1e-9 * Go * exp(-1.5 * Av) * u[7] -
59+
3e-10 * Go * exp(-3 * Av) * u[6] +
60+
1e-10 * Go * exp(-3 * Av) * u[10] * shield
61+
62+
# 7: CHx
63+
du[7] = n_H * (-2e-10) * u[7] * u[8] +
64+
n_H * 4e-16 * u[1] * u[12] +
65+
n_H * 2e-9 * u[2] * u[6] -
66+
1e-9 * Go * u[7] * exp(-1.5 * Av)
67+
68+
# 8: O
69+
du[8] = n_H * (-2e-10) * u[7] * u[8] +
70+
n_H * 1.6e-9 * u[10] * u[5] -
71+
n_H * 8e-10 * u[2] * u[8] +
72+
5e-10 * Go * exp(-1.7 * Av) * u[9] +
73+
1e-10 * Go * exp(-3 * Av) * u[10] * shield
74+
75+
# 9: OHx
76+
du[9] = n_H * (-1e-9) * u[9] * u[12] +
77+
n_H * 8e-10 * u[2] * u[8] -
78+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
79+
5e-10 * Go * exp(-1.7 * Av) * u[9]
80+
81+
# 10: CO
82+
du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
83+
n_H * 2e-10 * u[7] * u[8] -
84+
n_H * 1.7e-9 * u[10] * u[2] -
85+
n_H * 1.6e-9 * u[10] * u[5] +
86+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
87+
1e-10 * Go * exp(-3 * Av) * u[10] +
88+
1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
89+
90+
# 11: HCO+
91+
du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
92+
n_H * 1e-9 * u[9] * u[12] +
93+
n_H * 1.7e-9 * u[10] * u[2] -
94+
1.5e-10 * Go * exp(-2.5 * Av) * u[11]
95+
96+
# 12: C+
97+
du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
98+
n_H * 4e-16 * u[1] * u[12] -
99+
n_H * 1e-9 * u[9] * u[12] +
100+
n_H * 1.6e-9 * u[10] * u[5] +
101+
3e-10 * Go * exp(-3 * Av) * u[6]
102+
103+
# 13: M+
104+
du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
105+
n_H * 2e-9 * u[2] * u[14] +
106+
2.0e-10 * Go * exp(-1.9 * Av) * u[14]
107+
108+
# 14: M
109+
du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
110+
n_H * 2e-9 * u[2] * u[14] -
111+
2.0e-10 * Go * exp(-1.9 * Av) * u[14]
112+
end
113+
114+
# Set the Timespan, Parameters, and Initial Conditions
115+
seconds_per_year = 3600 * 24 * 365
116+
tspan = (0.0, 30000 * seconds_per_year) # ~30 thousand yrs
117+
118+
params = (10, # T
119+
2, # Av
120+
1.7, # Go
121+
611, # n_H
122+
1) # shield
123+
124+
u0 = [0.5, # 1: H2
125+
9.059e-9, # 2: H3+
126+
2.0e-4, # 3: e
127+
0.1, # 4: He
128+
7.866e-7, # 5: He+
129+
0.0, # 6: C
130+
0.0, # 7: CHx
131+
0.0004, # 8: O
132+
0.0, # 9: OHx
133+
0.0, # 10: CO
134+
0.0, # 11: HCO+
135+
0.0002, # 12: C+
136+
2.0e-7, # 13: M+
137+
2.0e-7] # 14: M
138+
139+
prob = ODEProblem(Nelson!, u0, tspan, params)
140+
sys = mtkcompile(modelingtoolkitize(prob))
141+
A, B, C = ModelingToolkit.calculate_semiquadratic_form(sys)
142+
@test A !== nothing
143+
@test B !== nothing
144+
@test C !== nothing
145+
x = unknowns(sys)
146+
linear_expr = A * x
147+
linear_fun, = generate_custom_function(sys, linear_expr; expression = Val{false})
148+
quadratic_expr = reduce(vcat, [x' * _B * x for _B in B if _B !== nothing])
149+
quadratic_fun, = generate_custom_function(sys, quadratic_expr; expression = Val{false})
150+
nonlinear_expr = C
151+
nonlinear_fun, = generate_custom_function(sys, nonlinear_expr; expression = Val{false})
152+
prob = ODEProblem(sys, nothing, tspan)
153+
linear_val = linear_fun(prob.u0, prob.p, 0.0)
154+
quadratic_val = quadratic_fun(prob.u0, prob.p, 0.0)
155+
nonlinear_val = nonlinear_fun(prob.u0, prob.p, 0.0)
156+
refsol = solve(prob, Vern9(); abstol = 1e-14, reltol = 1e-14)
157+
158+
@testset "stiff_A: $stiff_A, stiff_B: $stiff_B, stiff_C: $stiff_C" for (stiff_A, stiff_B, stiff_C) in Iterators.product(
159+
[false, true], [false, true], [false, true])
160+
kwargs = (; stiff_A, stiff_B, stiff_C)
161+
if stiff_A == stiff_B == stiff_C
162+
if stiff_A
163+
@test_throws ["All of", "cannot be stiff"] SemilinearODEProblem(
164+
sys, nothing, tspan; kwargs...)
165+
@test_throws ["All of", "cannot be stiff"] SemilinearODEFunction(
166+
sys; kwargs...)
167+
else
168+
@test_throws ["All of", "cannot be non-stiff"] SemilinearODEProblem(
169+
sys, nothing, tspan; kwargs...)
170+
@test_throws ["All of", "cannot be non-stiff"] SemilinearODEFunction(
171+
sys; kwargs...)
172+
end
173+
continue
174+
end
175+
176+
reference_f1 = zeros(length(u0))
177+
reference_f2 = zeros(length(u0))
178+
mul!(stiff_A ? reference_f1 : reference_f2, I, linear_val, true, true)
179+
mul!(stiff_B ? reference_f1 : reference_f2, I, quadratic_val, true, true)
180+
mul!(stiff_C ? reference_f1 : reference_f2, I, nonlinear_val, true, true)
181+
182+
@testset "Standard" begin
183+
prob = SemilinearODEProblem(sys, nothing, tspan; kwargs...)
184+
@test prob.f.f1(prob.u0, prob.p, 0.0)≈reference_f1 atol=1e-10
185+
@test prob.f.f2(prob.u0, prob.p, 0.0)≈reference_f2 atol=1e-10
186+
sol = solve(prob, KenCarp47())
187+
@test SciMLBase.successful_retcode(sol)
188+
@test refsol(sol.t).u≈sol.u atol=1e-8 rtol=1e-8
189+
end
190+
191+
@testset "Symbolic jacobian" begin
192+
prob = SemilinearODEProblem(sys, nothing, tspan; jac = true, kwargs...)
193+
@test prob.f.f1.jac !== nothing
194+
sol = solve(prob, KenCarp47())
195+
@test SciMLBase.successful_retcode(sol)
196+
@test refsol(sol.t).u≈sol.u atol=1e-8 rtol=1e-8
197+
end
198+
199+
@testset "Sparse" begin
200+
prob = SemilinearODEProblem(sys, nothing, tspan; sparse = true, kwargs...)
201+
sol = solve(prob, KenCarp47())
202+
@test SciMLBase.successful_retcode(sol)
203+
@test refsol(sol.t).u≈sol.u atol=1e-8 rtol=1e-8
204+
end
205+
206+
@testset "Sparsejac" begin
207+
@test_throws ["not implemented"] SemilinearODEProblem(
208+
sys, nothing, tspan; jac = true, sparse = true, kwargs...)
209+
end
210+
end
211+
end

0 commit comments

Comments
 (0)