Adaptive IDSolve #2881
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -4,38 +4,40 @@ authors = ["vyudu <[email protected]>"] | |
| version = "1.2.0" | ||
|
|
||
| [deps] | ||
| ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" | ||
| DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" | ||
| NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0" | ||
| NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d" | ||
| OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" | ||
| Reexport = "189a3867-3050-52da-a836-e630ba90ab69" | ||
| SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" | ||
| SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" | ||
|
|
||
| [sources] | ||
| OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"} | ||
|
|
||
| [compat] | ||
| AllocCheck = "0.2" | ||
| Aqua = "0.8.11" | ||
| ConcreteStructs = "0.2.3" | ||
| DiffEqBase = "6.176" | ||
| JET = "0.9.18, 0.10.4" | ||
| NonlinearSolveBase = "2.0.0" | ||
| NonlinearSolveFirstOrder = "1.9.0" | ||
| OrdinaryDiffEqCore = "1.29.0" | ||
| OrdinaryDiffEqSDIRK = "1.6.0" | ||
| Reexport = "1.2" | ||
| SciMLBase = "2.99" | ||
| SymbolicIndexingInterface = "0.3.38" | ||
| Test = "1.10.0" | ||
| julia = "1.10" | ||
|
|
||
| [extras] | ||
| AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" | ||
| Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" | ||
| JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" | ||
| OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" | ||
| Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" | ||
|
|
||
| [targets] | ||
| test = ["OrdinaryDiffEqSDIRK", "Test", "JET", "Aqua", "AllocCheck"] | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,17 @@ | ||
| """ | ||
| IDSolve() | ||
|
|
||
| First order solver for `ImplicitDiscreteSystems`. | ||
| """ | ||
| struct IDSolve{NLS} <: | ||
| OrdinaryDiffEqAlgorithm | ||
| nlsolve::NLS | ||
| extrapolant::Symbol | ||
| end | ||
|
|
||
| function IDSolve(; | ||
| nlsolve = NewtonRaphson(), | ||
| extrapolant = :constant, | ||
| ) | ||
| IDSolve{typeof(nlsolve)}(nlsolve, extrapolant) | ||
| end |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,36 +1,66 @@ | ||
| struct ImplicitDiscreteState{uType, pType, tType} | ||
| u::uType | ||
| p::pType | ||
| t::tType | ||
| end | ||
|
|
||
| mutable struct IDSolveCache{uType, cType, thetaType} <: OrdinaryDiffEqMutableCache | ||
| u::uType | ||
| uprev::uType | ||
| z::uType | ||
| nlcache::cType | ||
| Θks::thetaType | ||
| end | ||
|
|
||
| function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, | ||
| ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, | ||
| dt, reltol, p, calck, | ||
| ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} | ||
| state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t) | ||
| f_nl = (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t) | ||
|
There was a problem hiding this comment. What is the reasoning here to include the current There was a problem hiding this comment.
There was a problem hiding this comment. I guess my question is simply, why is |
||
|
|
||
| u_len = isnothing(u) ? 0 : length(u) | ||
| nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len) | ||
| unl = isnothing(u) ? Float64[] : u # FIXME nonlinear solve cannot handle nothing for u | ||
| prob = if nlls | ||
| NonlinearLeastSquaresProblem{isinplace(f)}( | ||
| NonlinearFunction(f_nl; resid_prototype = f.resid_prototype), | ||
| unl, state) | ||
| else | ||
| NonlinearProblem{isinplace(f)}(f_nl, unl, state) | ||
| end | ||
|
|
||
| nlcache = init(prob, alg.nlsolve) | ||
|
|
||
| IDSolveCache(u, uprev, state.u, nlcache, uBottomEltypeNoUnits[]) | ||
| end | ||
|
|
||
| isdiscretecache(cache::IDSolveCache) = true | ||
|
|
||
| function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, | ||
| ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, | ||
| dt, reltol, p, calck, | ||
| ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} | ||
| @assert !isnothing(u) "Empty u not supported with out of place functions yet." | ||
|
|
||
| state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t) | ||
| f_nl = (u_next, p) -> f(u_next, p.u, p.p, p.t) | ||
|
|
||
| u_len = isnothing(u) ? 0 : length(u) | ||
| nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len) | ||
| unl = isnothing(u) ? Float64[] : u # FIXME nonlinear solve cannot handle nothing for u | ||
| prob = if nlls | ||
| NonlinearLeastSquaresProblem{isinplace(f)}( | ||
| NonlinearFunction(f_nl; resid_prototype = f.resid_prototype), | ||
| unl, state) | ||
| else | ||
| NonlinearProblem{isinplace(f)}(f_nl, unl, state) | ||
| end | ||
|
|
||
| nlcache = init(prob, alg.nlsolve) | ||
|
|
||
| # FIXME Use IDSolveConstantCache? | ||
| IDSolveCache(u, uprev, state.u, nlcache, uBottomEltypeNoUnits[]) | ||
| end | ||
|
|
||
| get_fsalfirstlast(cache::IDSolveCache, rate_prototype) = (nothing, nothing) | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,89 @@ | ||
| """ | ||
| KantorovichTypeController | ||
|
|
||
| Default controller for implicit discrete solvers. Assuming a Newton method is used | ||
| to solve the nonlinear problem, this controller uses convergence rate estimates to | ||
| adapt the step size based on an posteriori estimate of the change in the Newton | ||
| convergence radius between steps as described below. | ||
|
|
||
| Given the convergence rate estimate Θ₀ for the first iteration, the step size controller | ||
| adapts the time step as `dtₙ₊₁ = γ * (g(Θbar) / (g(Θ₀)))^(1 / p) dtₙ` | ||
| with `g(x) = √(1 + 4x) - 1`. `p` denotes the order of the solver -- i.e. the order of | ||
| the extrapolation algorithm to compute the initial guess for the solve at `tₙ` given a | ||
| solution at `tₙ₋₁` -- and `Θbar` denotes the target convergence rate. `γ` is a safety factor. | ||
|
|
||
| A factor `Θreject` controls the rejection of a time step if any `Θₖ > Θreject`. | ||
| In this case the first Θₖ violating this criterion is taken and the time step is adapted | ||
| such that `dtₙ₊₁ = γ * (g(Θbar) / (g(Θk)))^(1 / p) dtₙ`. This behavior can be changed | ||
| by setting `strict = false`. In this case the step is accepted whenever the Newton | ||
| solver converges. | ||
|
|
||
| The controller furthermore limits the growth and shrinkage of the time step by a factor | ||
| between `qmin` and `qmax`. | ||
|
|
||
| The baseline algorithm has been derived in Peter Deuflhard's book "Newton Methods for | ||
| Nonlinear Problems" in Section 5.1.3 (Adaptive pathfollowing algorithms). Please note | ||
| that some implementation details deviate from the original algorithm. | ||
| """ | ||
| Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractController | ||
|
There was a problem hiding this comment. TODO reference and documentation There was a problem hiding this comment. Yeah I'm not sure what this is. There was a problem hiding this comment. Oh, sorry. I will write the docs, no worries. I just left it here so I do not forget before merging. This is a controller derived from a posteriori estimates on how much the convergence radius in the Newton-Kantorovich theorem changes for some increment There was a problem hiding this comment. I have added a docstring with a brief explanation and reference now. |
||
| Θmin::Float64 | ||
| p::Int64 | ||
| Θreject::Float64 = 0.95 | ||
| Θbar::Float64 = 0.5 | ||
| γ::Float64 = 0.95 | ||
| qmin::Float64 = 1 / 5 | ||
| qmax::Float64 = 5.0 | ||
| strict::Bool = true | ||
| end | ||
|
|
||
| function OrdinaryDiffEqCore.default_controller( | ||
| alg::IDSolve, cache::IDSolveCache, _1, _2, _3) | ||
| return KantorovichTypeController(; Θmin = 1 // 8, p = 1) | ||
| end | ||
|
|
||
| function OrdinaryDiffEqCore.stepsize_controller!( | ||
| integrator, controller::KantorovichTypeController, alg::IDSolve | ||
| ) | ||
| @inline g(x) = √(1 + 4x) - 1 | ||
|
|
||
| # Adapt dt with a priori estimate (Eq. 5.24) | ||
| (; Θks) = integrator.cache | ||
| (; Θbar, γ, Θmin, qmin, qmax, p) = controller | ||
|
|
||
| Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin | ||
| q = clamp(γ * (g(Θbar) / (g(Θ₀)))^(1 / p), qmin, qmax) | ||
|
|
||
| return q | ||
termi-official marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| end | ||
|
|
||
| function OrdinaryDiffEqCore.step_accept_controller!( | ||
| integrator, controller::KantorovichTypeController, alg::IDSolve, q | ||
| ) | ||
| return q * integrator.dt | ||
| end | ||
|
|
||
| function OrdinaryDiffEqCore.step_reject_controller!( | ||
| integrator, controller::KantorovichTypeController, alg::IDSolve | ||
| ) | ||
| @inline g(x) = √(1 + 4x) - 1 | ||
|
|
||
| # Shorten dt according to (Eq. 5.24) | ||
| (; Θks) = integrator.cache | ||
| (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller | ||
| for Θk in Θks | ||
| if Θk > Θreject | ||
| q = clamp(γ * (g(Θbar) / g(Θk))^(1 / p), qmin, qmax) | ||
| integrator.dt = q * integrator.dt | ||
| return | ||
| end | ||
| end | ||
| end | ||
|
|
||
| function OrdinaryDiffEqCore.accept_step_controller(integrator, controller::KantorovichTypeController) | ||
| (; Θks) = integrator.cache | ||
| if controller.strict | ||
| return all(controller.Θreject .< Θks) | ||
| else | ||
| return true | ||
| end | ||
| end | ||
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
why 1 here?
Uh oh!
There was an error while loading. Please reload this page.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This comes from the analysis shown in Deuflhard, Newton Methods for Nonlinear Problems (Section 5.1.1, see Equation 5.6 and the surrounding definition).
The algorithms here have an associated order in the sense that for a given$dt_n = t_{n+1} - t_n$ we have for some solution $\hat{u}(t_n+1)$ derived from an initial guess given at $u({t_n})$ . Now we can define an associated ODE (Davidenko differential equation*) for each nonlinear problem with a "time parameter" by taking the time derivative of the time parameter, which has an analytical solution $\bar{u}(t_{n+1})$ , given the same initial guess ($u({t_n})$). The solver now has order $p$ if we have $$||\hat{u}(t_{n+1}) - \bar{u}(t_{n+1})|| \leq C dt^p_n . $$
Does that explain it?
*The Davidenko differential equation for$F(u,t)$ is simply $du/dt = - dF/dx (u,t)^{-1} * dF/dt (u,t)$ .
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess I'm confused. It's a discrete time problem, it's exact?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay, let me try to explain it differently then. We have the parametric function$F(u,t)$ and want to find the solution of $F(u_2,t_2) = 0$ given a $u_1$ such that $F(u_1, t_1) = 0$ . With with $t_1 < t_2$ we want to find some initial guess for $u^0_2$ given $u_1$ , such that the initial guess $u^0_2$ is inside the convergence radius of a Newton method to solve $F(u_2,t_2) = 0$ . The obvious choice is that we can simply say our initial guess is simply $u_1$ , but this is typically not really great, as $t_2 - t_1$ is often quite small. However, we can use additional information contained in $F(u,t) = 0$ . To be specific the derivative with respect to the parameter $t$ contains some extra information which we can use to improve the initial guess. Here we can observe that analytically solving the associated Davidenko differential equation with initial condition $u_1$ on $[t_1, t_2]$ is equivalent to solving $F(u_2,t_2) = 0$ . Furthermore, we can exploit this information to inform how large $t_2$ can be chosen, such that the Newton method is guaranteed converge for a given initial guess. Now, the order of this extrapolation polynomial is directly related to the order of the implicit discrete solver. I hope that helps.