Skip to content

Commit 5bb0876

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #2458 from SciML/daeinitialize_error
Better warning on lack of initialization algorithm and slim Rosenbrock
2 parents 7a98e43 + e8e8540 commit 5bb0876

File tree

8 files changed

+77
-57
lines changed

8 files changed

+77
-57
lines changed

‎docs/src/massmatrixdae/Rosenbrock.md

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,13 @@ At medium tolerances `>1e-8` it is recommended you use `Rodas5P` or `Rodas4P`,
1010
the former is more efficient, but the latter is more reliable.
1111
For larger systems look at multistep methods.
1212

13+
!!! warn
14+
15+
In order to use OrdinaryDiffEqRosenbrock with DAEs that require a non-trivial
16+
consistent initialization, a nonlinear solver is required and thus
17+
`using OrdinaryDiffEqNonlinearSolve` is required or you must pass an `initializealg`
18+
with a valid `nlsolve` choice.
19+
1320
## Example usage
1421

1522
```julia

‎lib/OrdinaryDiffEqCore/Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "OrdinaryDiffEqCore"
22
uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
33
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
4-
version = "1.4.1"
4+
version = "1.5.0"
55

66
[deps]
77
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

‎lib/OrdinaryDiffEqCore/src/initialize_dae.jl

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -104,8 +104,65 @@ function _initialize_dae!(integrator, prob::DAEProblem,
104104
end
105105
end
106106

107+
## Nonlinear Solver Defaulting
108+
109+
## If an alg is given use it
110+
default_nlsolve(alg, isinplace, u, initprob, autodiff = false) = alg
111+
112+
## If the initialization is trivial just use nothing alg
113+
function default_nlsolve(
114+
::Nothing, isinplace, u::Nothing, ::NonlinearProblem, autodiff = false)
115+
nothing
116+
end
117+
118+
function default_nlsolve(
119+
::Nothing, isinplace, u::Nothing, ::NonlinearLeastSquaresProblem, autodiff = false)
120+
nothing
121+
end
122+
123+
function OrdinaryDiffEqCore.default_nlsolve(::Nothing, isinplace::Val{true}, u, ::NonlinearProblem, autodiff = false)
124+
error("This ODE requires a DAE initialization and thus a nonlinear solve but no nonlinear solve has been loaded. To solve this problem, do `using OrdinaryDiffEqNonlinearSolve` or pass a custom `nlsolve` choice into the `initializealg`.")
125+
end
126+
127+
function OrdinaryDiffEqCore.default_nlsolve(::Nothing, isinplace::Val{true}, u, ::NonlinearLeastSquaresProblem, autodiff = false)
128+
error("This ODE requires a DAE initialization and thus a nonlinear solve but no nonlinear solve has been loaded. To solve this problem, do `using OrdinaryDiffEqNonlinearSolve` or pass a custom `nlsolve` choice into the `initializealg`.")
129+
end
130+
107131
## NoInit
108132

109133
function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
110134
alg::NoInit, x::Union{Val{true}, Val{false}})
111135
end
136+
137+
## OverrideInit
138+
139+
function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
140+
alg::OverrideInit, isinplace::Union{Val{true}, Val{false}})
141+
initializeprob = prob.f.initializeprob
142+
143+
# If it doesn't have autodiff, assume it comes from symbolic system like ModelingToolkit
144+
# Since then it's the case of not a DAE but has initializeprob
145+
# In which case, it should be differentiable
146+
isAD = if initializeprob.u0 === nothing
147+
AutoForwardDiff
148+
elseif has_autodiff(integrator.alg)
149+
alg_autodiff(integrator.alg) isa AutoForwardDiff
150+
else
151+
true
152+
end
153+
154+
alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD)
155+
nlsol = solve(initializeprob, alg)
156+
if isinplace === Val{true}()
157+
integrator.u .= prob.f.initializeprobmap(nlsol)
158+
elseif isinplace === Val{false}()
159+
integrator.u = prob.f.initializeprobmap(nlsol)
160+
else
161+
error("Unreachable reached. Report this error.")
162+
end
163+
164+
if nlsol.retcode != ReturnCode.Success
165+
integrator.sol = SciMLBase.solution_new_retcode(integrator.sol,
166+
ReturnCode.InitialFailure)
167+
end
168+
end

‎lib/OrdinaryDiffEqNonlinearSolve/Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "OrdinaryDiffEqNonlinearSolve"
22
uuid = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
33
authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <mayingbo5@gmail.com>"]
4-
version = "1.1.0"
4+
version = "1.2.0"
55

66
[deps]
77
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

‎lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,10 @@ import OrdinaryDiffEqCore
2727
import SciMLOperators: islinear
2828
import OrdinaryDiffEqCore: nlsolve_f, set_new_W!, set_W_γdt!
2929

30+
@static if isdefined(OrdinaryDiffEqCore, :default_nlsolve)
31+
import OrdinaryDiffEqCore: default_nlsolve
32+
end
33+
3034
using OrdinaryDiffEqCore: resize_nlsolver!, _initialize_dae!,
3135
AbstractNLSolverAlgorithm, AbstractNLSolverCache,
3236
AbstractNLSolver, NewtonAlgorithm, @unpack,

‎lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl

Lines changed: 6 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -1,64 +1,20 @@
1-
default_nlsolve(alg, isinplace, u, initprob, autodiff = false) = alg
2-
3-
function default_nlsolve(
4-
::Nothing, isinplace, u::Nothing, ::NonlinearProblem, autodiff = false)
5-
nothing
6-
end
7-
function default_nlsolve(::Nothing, isinplace, u, ::NonlinearProblem, autodiff = false)
1+
function default_nlsolve(::Nothing, isinplace::Val{true}, u, ::NonlinearProblem, autodiff = false)
82
FastShortcutNonlinearPolyalg(;
93
autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
104
end
11-
function default_nlsolve(::Nothing, isinplace::Val{false}, u::StaticArray,
12-
::NonlinearProblem, autodiff = false)
13-
SimpleTrustRegion(autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
14-
end
15-
165
function default_nlsolve(
17-
::Nothing, isinplace, u::Nothing, ::NonlinearLeastSquaresProblem, autodiff = false)
18-
nothing
19-
end
20-
function default_nlsolve(
21-
::Nothing, isinplace, u, ::NonlinearLeastSquaresProblem, autodiff = false)
6+
::Nothing, isinplace::Val{true}, u, ::NonlinearLeastSquaresProblem, autodiff = false)
227
FastShortcutNLLSPolyalg(; autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
238
end
9+
function default_nlsolve(::Nothing, isinplace::Val{false}, u::StaticArray,
10+
::NonlinearProblem, autodiff = false)
11+
SimpleTrustRegion(autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
12+
end
2413
function default_nlsolve(::Nothing, isinplace::Val{false}, u::StaticArray,
2514
::NonlinearLeastSquaresProblem, autodiff = false)
2615
SimpleGaussNewton(autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
2716
end
2817

29-
## OverrideInit
30-
31-
function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
32-
alg::OverrideInit, isinplace::Union{Val{true}, Val{false}})
33-
initializeprob = prob.f.initializeprob
34-
35-
# If it doesn't have autodiff, assume it comes from symbolic system like ModelingToolkit
36-
# Since then it's the case of not a DAE but has initializeprob
37-
# In which case, it should be differentiable
38-
isAD = if initializeprob.u0 === nothing
39-
AutoForwardDiff
40-
elseif has_autodiff(integrator.alg)
41-
alg_autodiff(integrator.alg) isa AutoForwardDiff
42-
else
43-
true
44-
end
45-
46-
alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD)
47-
nlsol = solve(initializeprob, alg)
48-
if isinplace === Val{true}()
49-
integrator.u .= prob.f.initializeprobmap(nlsol)
50-
elseif isinplace === Val{false}()
51-
integrator.u = prob.f.initializeprobmap(nlsol)
52-
else
53-
error("Unreachable reached. Report this error.")
54-
end
55-
56-
if nlsol.retcode != ReturnCode.Success
57-
integrator.sol = SciMLBase.solution_new_retcode(integrator.sol,
58-
ReturnCode.InitialFailure)
59-
end
60-
end
61-
6218
## ShampineCollocationInit
6319

6420
#=

‎lib/OrdinaryDiffEqRosenbrock/Project.toml

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "OrdinaryDiffEqRosenbrock"
22
uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
33
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
4-
version = "1.1.1"
4+
version = "1.2.0"
55

66
[deps]
77
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -15,7 +15,6 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
1515
MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
1616
OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
1717
OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b"
18-
OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
1918
Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
2019
PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
2120
Preferences = "21216c6a-2e73-6563-6e65-726566657250"
@@ -37,7 +36,6 @@ MuladdMacro = "0.2.4"
3736
ODEProblemLibrary = "0.1.8"
3837
OrdinaryDiffEqCore = "1.1"
3938
OrdinaryDiffEqDifferentiation = "<0.0.1, 1"
40-
OrdinaryDiffEqNonlinearSolve = "<0.0.1, 1"
4139
Polyester = "0.7.16"
4240
PrecompileTools = "1.2.1"
4341
Preferences = "1.4.3"

‎lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -33,8 +33,6 @@ using OrdinaryDiffEqDifferentiation: TimeDerivativeWrapper, TimeGradientWrapper,
3333
calc_W, calc_rosenbrock_differentiation!, build_J_W,
3434
UJacobianWrapper, dolinsolve
3535

36-
import OrdinaryDiffEqNonlinearSolve # Required for DAE initialization
37-
3836
using Reexport
3937
@reexport using DiffEqBase
4038

0 commit comments

Comments
 (0)