Adaptive IDSolve #2881

termi-official · 2025-10-06T17:17:14Z

Checklist

Appropriate tests were added
Any code changes were done in a way that does not break public API
All documentation related to code changes were updated
The new code follows the
contributor guidelines, in particular the SciML Style Guide and
COLPRAC.
Any new documentation only uses public API

Additional context

Switch from hard-coded SimpleNonlinearSolve to any NonlinearSolve algorithm which supports the iterator interface.
Add adaptive homotopy path increments via modified Deuflhard's estimates

…ds generic adaptivity

termi-official · 2025-10-06T17:17:53Z

lib/ImplicitDiscreteSolve/test/runtests.jl

-    function empty(u_next, u, p, t)
-        nothing
-    end
+# @testset "Handle nothing in u0" begin


I have trouble understanding the purpose of this test. What exactly is the practical scenario here?

we use u=nothing to represent length 0 u in a type stable manner (as op0osed to a Vector which is length 0 at runtime)

Sorry, I think I do not understand the point - do you imply here that a zero-length solution vector not type stable and this is why we resort to "nothing"?

When I try to get this test running, I get errors from an alias analysis function not being dispatched for nothing in NonlinearSolve (https://github.com/SciML/NonlinearSolve.jl/blob/ac9344f9359833282e443c4479427ad9ce3311dd/lib/NonlinearSolveFirstOrder/src/solve.jl#L157).

Inplace functions with empty return types are also not correctly handled.

using NonlinearSolveFirstOrder function f(u, p) nothing end prob = NonlinearProblem{false}(f, Float64[], nothing) iter = init(prob, NewtonRaphson())

errors when the termination cache is built, because the increment type cannot be derived.

the point is that zero length u tends to cause problems (e.g. solvers will try to look at the first element of the state to find a type), so this way you can skip the solve process since nothing interesting can happen if you don't have any state. See https://github.com/SciML/DiffEqBase.jl/blob/3667bdbdc85489f7b296316df7f4c440519e82f6/src/solve.jl#L31 for how this gets handled for ODEs/DAEs.

Wouldn't it make more sense from a software engineering stand point to replace such isa statements with dispatchable functions to make the code more extensible?

quite possibly. DiffEqBase is not my favorite code organization.

…timates.

termi-official · 2025-10-06T21:23:49Z

lib/ImplicitDiscreteSolve/src/alg_utils.jl

+# @concrete struct ConvergenceRateTracing
+#     inner_tracing
+# end
+
+# @concrete struct ConvergenceRateTraceTrick
+#     incrementL2norms
+#     residualL2norms
+#     trace_wrapper
+# end
+
+# function NonlinearSolveBase.init_nonlinearsolve_trace(
+#         prob, alg::IDSolve, u, fu, J, δu;
+#         trace_level::ConvergenceRateTracing, kwargs... # This kind of dispatch does not work. Need to figure out a different way.
+# )
+#     inner_trace = NonlinearSolveBase.init_nonlinearsolve_trace(
+#         prob, alg, u, fu, J, δu;
+#         trace_level.inner_tracing, kwargs...
+#     )
+
+#     return ConvergenceRateTraceTrick(eltype(δu)[], eltype(fu)[], inner_trace)
+# end


From my understanding, it should be possible to use NonlinearSolveBase.init_nonlinearsolve_trace to query convergence rate estimates. However, in the current design I cannot add a new dispatch. Should I make a PR to pull the trace level before the kwargs (i.e. NonlinearSolveBase.init_nonlinearsolve_trace(prob, alg, u, fu, J, δu, trace_level; kwargs...) for custom dispatches?

termi-official · 2025-10-06T21:24:55Z

lib/ImplicitDiscreteSolve/src/cache.jl

-    state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t)
-    IDSolveCache(u, uprev, state, nothing)
+    state = ImplicitDiscreteState(zero(u), p, t)
+    f_nl = (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t)


What is the reasoning here to include the current u in the signature of this function, but no information on dt?

dt is just a parameter?

I guess my question is simply, why is $dt$ (or tprev) not a parameter, but $uprev$ is part of the function signature?

lib/ImplicitDiscreteSolve/src/cache.jl

lib/ImplicitDiscreteSolve/src/controller.jl

termi-official · 2025-10-06T21:28:47Z

lib/ImplicitDiscreteSolve/src/controller.jl

@@ -0,0 +1,52 @@
+Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractController


TODO reference and documentation

Yeah I'm not sure what this is.

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 $dt_n$ and a solution given at $t_n$.

termi-official · 2025-10-06T21:29:40Z

lib/ImplicitDiscreteSolve/src/solve.jl

+    else # :constant
+        cache.z .= integrator.u
+    end
+    state = ImplicitDiscreteState(cache.z, p, t+dt)


On master we solve at time $t$. From my understanding, we should solve at $t+dt$ to obtain the next solution. Can someone confirm or reject this?

the somewhat complicated part here is that we want to match a DiscreteSolve. I forget which is the right way here. @jClugstor might remember.

From what I can tell from the docs https://docs.sciml.ai/DiffEqDocs/stable/types/discrete_types/
if you put in u_n, p, t_(n+1) to the function you get u_(n+1) out, so in this case you would get u_(t + dt) out I guess. So if you want u_(t + dt) I think that's correct. Not entirely sure though.

ChrisRackauckas · 2025-10-23T22:57:56Z

lib/ImplicitDiscreteSolve/src/alg_utils.jl


 isfsal(alg::IDSolve) = false
-alg_order(alg::IDSolve) = 0
+alg_order(alg::IDSolve) = 1


why 1 here?

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)$.

I guess I'm confused. It's a discrete time problem, it's exact?

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.

ChrisRackauckas · 2025-10-23T23:02:13Z

lib/ImplicitDiscreteSolve/src/solve.jl

+    resize!(Θks, 0)
+    residualnormprev = zero(eltype(u))
+    while NonlinearSolveBase.not_terminated(nlcache)
+        step!(nlcache)


should it also be trying jacobian reuse in here?

Yes. However, shoudln't the reuse strategy should be part of the NonlinearSolve algorithm instead of some intermediate layer?

if you expand out the iterator then you're taking over the iteration strategy.

Correct. I would like to remove that later. Right now I cannot use the high level API because I cannot access the convergence rates which I need in the controller.

termi-official added 2 commits October 6, 2025 19:10

Prototype for using NonlinearSolve in IDSolve as the first step towar…

d84627c

…ds generic adaptivity

Reenable JET tests

f02492f

termi-official commented Oct 6, 2025

View reviewed changes

termi-official added 4 commits October 6, 2025 21:50

Add poor mans adaptivity to get the ball rolling

c96581c

Add a modified variant of Deuflhard's controller using Kantorovich es…

e54472a

…timates.

Cleanup debris and rejection implementation

cd0fbf6

:)

c25238c

termi-official commented Oct 6, 2025

View reviewed changes

termi-official marked this pull request as ready for review October 6, 2025 21:31

oscardssmith self-requested a review October 7, 2025 05:54

termi-official added 2 commits October 7, 2025 16:29

Add partial handling of empty u

e9a58fa

Remove debug message

bdfad62

termi-official mentioned this pull request Oct 9, 2025

User-defined Traces SciML/NonlinearSolve.jl#715

Open

Format

4b597fb

oscardssmith requested a review from ChrisRackauckas October 20, 2025 13:34

ChrisRackauckas reviewed Oct 23, 2025

View reviewed changes

		@@ -0,0 +1,52 @@
		Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractController

Uh oh!

Uh oh!

Adaptive IDSolve #2881

Are you sure you want to change the base?

Adaptive IDSolve #2881

Conversation

termi-official commented Oct 6, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Checklist

Additional context

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

termi-official Oct 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

termi-official Oct 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

termi-official commented Oct 6, 2025 •

edited

Loading

termi-official Oct 24, 2025 •

edited

Loading

termi-official Oct 24, 2025 •

edited

Loading