Skip to content

Commit ec4aed2

Browse files
ChrisRackauckasChrisRackauckas
committed
Add stiff initial step size algorithm for implicit solvers
Implement CVHin algorithm from SUNDIALS CVODE for initial step size selection in stiff solvers. This component-wise iterative algorithm is more robust than Hairer-Wanner for multi-scale problems where some state variables start at zero with tiny absolute tolerances (fixes #1496). Key changes: - Add StiffInitDt algorithm type with initdt_alg trait dispatch - All implicit solvers (OrdinaryDiffEqAdaptiveImplicitAlgorithm, OrdinaryDiffEqImplicitAlgorithm, DAEAlgorithm) use StiffInitDt - Explicit solvers continue using DefaultInitDt (Hairer-Wanner) - DAE problems (mass-matrix and AbstractDAEProblem) use IDA-style h = 0.001 * tdist for initial step size - Fall back to DefaultInitDt for non-Array types (GPU) and non-AbstractFloat element types (ForwardDiff Duals, Complex) Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
1 parent 77e805d commit ec4aed2

6 files changed

Lines changed: 670 additions & 11 deletions

File tree

‎lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,8 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
2323
get_fsalfirstlast, generic_solver_docstring, _bool_to_ADType,
2424
_process_AD_choice,
2525
_ode_interpolant, _ode_interpolant!, has_stiff_interpolation,
26-
_ode_addsteps!, DerivativeOrderNotPossibleError
26+
_ode_addsteps!, DerivativeOrderNotPossibleError,
27+
initdt_alg
2728
using OrdinaryDiffEqSDIRK: ImplicitEulerConstantCache, ImplicitEulerCache
2829

2930
using TruncatedStacktraces: @truncate_stacktrace

‎lib/OrdinaryDiffEqBDF/test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ if TEST_GROUP != "QA"
1313
@time @safetestset "BDF Inference Tests" include("inference_tests.jl")
1414
@time @safetestset "BDF Convergence Tests" include("bdf_convergence_tests.jl")
1515
@time @safetestset "BDF Regression Tests" include("bdf_regression_tests.jl")
16+
@time @safetestset "StiffInitDt Tests" include("stiff_initdt_tests.jl")
1617
end
1718

1819
# Run QA tests (AllocCheck, JET, Aqua) - skip on pre-release Julia

‎lib/OrdinaryDiffEqBDF/test/stiff_initdt_tests.jl‎

Lines changed: 137 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
using Test
2+
using OrdinaryDiffEqBDF
3+
using OrdinaryDiffEqSDIRK
4+
using OrdinaryDiffEqCore: initdt_alg, DefaultInitDt, StiffInitDt
5+
6+
@testset "StiffInitDt Algorithm" begin
7+
8+
@testset "Trait dispatch" begin
9+
# All implicit methods should use StiffInitDt
10+
@test initdt_alg(FBDF()) isa StiffInitDt
11+
@test initdt_alg(QNDF()) isa StiffInitDt
12+
@test initdt_alg(QNDF1()) isa StiffInitDt
13+
@test initdt_alg(QNDF2()) isa StiffInitDt
14+
@test initdt_alg(ABDF2()) isa StiffInitDt
15+
@test initdt_alg(ImplicitEuler()) isa StiffInitDt
16+
end
17+
18+
@testset "Simple exponential decay (in-place)" begin
19+
function f_exp!(du, u, p, t)
20+
du[1] = -u[1]
21+
end
22+
prob = ODEProblem(f_exp!, [1.0], (0.0, 10.0))
23+
sol = solve(prob, FBDF())
24+
@test sol.retcode == ReturnCode.Success
25+
@test isapprox(sol.u[end][1], exp(-10.0), rtol = 1.0e-2)
26+
end
27+
28+
@testset "Simple exponential decay (out-of-place)" begin
29+
f_exp(u, p, t) = [-u[1]]
30+
prob = ODEProblem(f_exp, [1.0], (0.0, 10.0))
31+
sol = solve(prob, FBDF())
32+
@test sol.retcode == ReturnCode.Success
33+
@test isapprox(sol.u[end][1], exp(-10.0), rtol = 1.0e-2)
34+
end
35+
36+
@testset "Multi-scale with zero states and tiny abstol (in-place)" begin
37+
# This is the pathological case from issue #1496:
38+
# Some species start at 0 with tiny abstol, causing the Hairer algorithm
39+
# to produce catastrophically small initial dt.
40+
function f_ms!(du, u, p, t)
41+
du[1] = -100.0 * u[1] + 1.0
42+
du[2] = 0.1 * u[1]
43+
du[3] = 1.84 # mimics TH2S
44+
end
45+
u0 = [1.0, 0.0, 0.0]
46+
prob = ODEProblem(f_ms!, u0, (0.0, 100.0))
47+
48+
# With FBDF (StiffInitDt) - should succeed efficiently
49+
sol_fbdf = solve(prob, FBDF(), abstol = 1.0e-15, reltol = 1.0e-8)
50+
@test sol_fbdf.retcode == ReturnCode.Success
51+
@test sol_fbdf.t[end] == 100.0
52+
53+
# The StiffInitDt should produce an initial dt much larger than the
54+
# catastrophically small ~5e-25 that the Hairer algorithm would give.
55+
# With abstol=1e-15 and u0=0, the StiffInitDt computes ~3.5e-15
56+
# (10 orders of magnitude larger).
57+
dt_fbdf = sol_fbdf.t[2] - sol_fbdf.t[1]
58+
@test dt_fbdf > 1.0e-20 # Much larger than the ~5e-25 Hairer would give
59+
end
60+
61+
@testset "Multi-scale with zero states and tiny abstol (out-of-place)" begin
62+
function f_ms(u, p, t)
63+
return [-100.0 * u[1] + 1.0, 0.1 * u[1], 1.84]
64+
end
65+
u0 = [1.0, 0.0, 0.0]
66+
prob = ODEProblem(f_ms, u0, (0.0, 100.0))
67+
68+
sol = solve(prob, FBDF(), abstol = 1.0e-15, reltol = 1.0e-8)
69+
@test sol.retcode == ReturnCode.Success
70+
@test sol.t[end] == 100.0
71+
end
72+
73+
@testset "Stiff Robertson problem" begin
74+
# Classic stiff test problem
75+
function robertson!(du, u, p, t)
76+
du[1] = -0.04 * u[1] + 1.0e4 * u[2] * u[3]
77+
du[2] = 0.04 * u[1] - 1.0e4 * u[2] * u[3] - 3.0e7 * u[2]^2
78+
du[3] = 3.0e7 * u[2]^2
79+
end
80+
u0 = [1.0, 0.0, 0.0]
81+
prob = ODEProblem(robertson!, u0, (0.0, 1.0e5))
82+
83+
sol = solve(prob, FBDF(), abstol = 1.0e-8, reltol = 1.0e-8)
84+
@test sol.retcode == ReturnCode.Success
85+
@test sol.t[end] == 1.0e5
86+
87+
# Conservation law: u1 + u2 + u3 = 1
88+
@test isapprox(sum(sol.u[end]), 1.0, atol = 1.0e-6)
89+
end
90+
91+
@testset "Backward integration" begin
92+
function f_back!(du, u, p, t)
93+
du[1] = -u[1]
94+
end
95+
prob = ODEProblem(f_back!, [exp(-10.0)], (10.0, 0.0))
96+
sol = solve(prob, FBDF())
97+
@test sol.retcode == ReturnCode.Success
98+
@test isapprox(sol.u[end][1], 1.0, rtol = 2.0e-2)
99+
end
100+
101+
@testset "User-specified dt still works" begin
102+
function f_dt!(du, u, p, t)
103+
du[1] = -u[1]
104+
end
105+
prob = ODEProblem(f_dt!, [1.0], (0.0, 1.0))
106+
# When user specifies dt, StiffInitDt should not be called
107+
sol = solve(prob, FBDF(), dt = 0.01)
108+
@test sol.retcode == ReturnCode.Success
109+
end
110+
111+
@testset "QNDF with multi-scale problem" begin
112+
function f_qndf!(du, u, p, t)
113+
du[1] = -100.0 * u[1] + 1.0
114+
du[2] = 0.0 # starts at 0, stays at 0
115+
end
116+
u0 = [1.0, 0.0]
117+
prob = ODEProblem(f_qndf!, u0, (0.0, 1.0))
118+
119+
sol = solve(prob, QNDF(), abstol = 1.0e-12, reltol = 1.0e-8)
120+
@test sol.retcode == ReturnCode.Success
121+
@test sol.t[end] == 1.0
122+
end
123+
124+
@testset "ImplicitEuler with StiffInitDt" begin
125+
# ImplicitEuler now also uses StiffInitDt since all implicit methods do
126+
function f_ie!(du, u, p, t)
127+
du[1] = -100.0 * u[1] + 1.0
128+
du[2] = 0.0
129+
end
130+
u0 = [1.0, 0.0]
131+
prob = ODEProblem(f_ie!, u0, (0.0, 1.0))
132+
133+
sol = solve(prob, ImplicitEuler(), abstol = 1.0e-12, reltol = 1.0e-8)
134+
@test sol.retcode == ReturnCode.Success
135+
@test sol.t[end] == 1.0
136+
end
137+
end

‎lib/OrdinaryDiffEqCore/src/alg_utils.jl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -186,6 +186,16 @@ isimplicit(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = true
186186
isimplicit(alg::OrdinaryDiffEqImplicitAlgorithm) = true
187187
isimplicit(alg::CompositeAlgorithm) = any(isimplicit.(alg.algs))
188188

189+
# Initial dt selection algorithm trait.
190+
# DefaultInitDt uses the Hairer-Wanner algorithm (suitable for non-stiff and general problems).
191+
# StiffInitDt uses a component-wise iterative algorithm (more robust for stiff multi-scale problems).
192+
# All implicit and DAE solvers default to StiffInitDt; explicit solvers use DefaultInitDt.
193+
initdt_alg(alg::OrdinaryDiffEqAlgorithm) = DefaultInitDt()
194+
initdt_alg(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = StiffInitDt()
195+
initdt_alg(alg::OrdinaryDiffEqImplicitAlgorithm) = StiffInitDt()
196+
initdt_alg(alg::DAEAlgorithm) = StiffInitDt()
197+
initdt_alg(alg::CompositeAlgorithm) = initdt_alg(alg.algs[1])
198+
189199
isdtchangeable(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = true
190200
isdtchangeable(alg::CompositeAlgorithm) = all(isdtchangeable.(alg.algs))
191201

‎lib/OrdinaryDiffEqCore/src/algorithms.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,8 @@
1+
# Initial dt selection algorithms
2+
abstract type InitDtAlg end
3+
struct DefaultInitDt <: InitDtAlg end
4+
struct StiffInitDt <: InitDtAlg end
5+
16
abstract type OrdinaryDiffEqAlgorithm <: SciMLBase.AbstractODEAlgorithm end
27
abstract type OrdinaryDiffEqAdaptiveAlgorithm <: OrdinaryDiffEqAlgorithm end
38
abstract type OrdinaryDiffEqCompositeAlgorithm <: OrdinaryDiffEqAlgorithm end

0 commit comments

Comments
 (0)