Skip to content

Commit b8c0e2c

Browse files
ChrisRackauckasChrisRackauckas
committed
Add stiff initial step size algorithm for implicit solvers
Implement a component-wise iterative initial step size algorithm (StiffInitDt) based on the CVHin algorithm from SUNDIALS CVODE (Hindmarsh et al., 2005). This replaces the Hairer-Wanner algorithm for all implicit/stiff solvers. The Hairer algorithm uses RMS norms where problematic components get diluted by 1/sqrt(N), producing catastrophically small initial dt when species start at zero with tiny abstol (e.g., dt ~ 5e-25 for the reactive-transport problem in issue #1496). The new algorithm uses component-wise max-norm bounds and iterative refinement, producing ~10 orders of magnitude larger initial dt. Key changes: - `InitDtAlg` type hierarchy with `DefaultInitDt` and `StiffInitDt` - `initdt_alg(alg)` trait: all implicit solvers default to `StiffInitDt`, explicit solvers use `DefaultInitDt` - Both in-place and out-of-place implementations - Handles mass matrices via linsolve with try/catch fallback for singular matrices (DAEs) Fixes #1496 Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
1 parent 6f57229 commit b8c0e2c

6 files changed

Lines changed: 625 additions & 1 deletion

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("sundials_initdt_tests.jl")
1617
end
1718

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

‎lib/OrdinaryDiffEqBDF/test/sundials_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: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -176,6 +176,15 @@ isimplicit(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = true
176176
isimplicit(alg::OrdinaryDiffEqImplicitAlgorithm) = true
177177
isimplicit(alg::CompositeAlgorithm) = any(isimplicit.(alg.algs))
178178

179+
# Initial dt selection algorithm trait.
180+
# DefaultInitDt uses the Hairer-Wanner algorithm (suitable for non-stiff and general problems).
181+
# StiffInitDt uses a component-wise iterative algorithm (more robust for stiff multi-scale problems).
182+
# All implicit solvers default to StiffInitDt; explicit solvers use DefaultInitDt.
183+
initdt_alg(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = DefaultInitDt()
184+
initdt_alg(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = StiffInitDt()
185+
initdt_alg(alg::OrdinaryDiffEqImplicitAlgorithm) = StiffInitDt()
186+
initdt_alg(alg::CompositeAlgorithm) = initdt_alg(alg.algs[1])
187+
179188
isdtchangeable(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = true
180189
isdtchangeable(alg::CompositeAlgorithm) = all(isdtchangeable.(alg.algs))
181190

‎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)