Implicit solvers

Some docs here: https://diffeq.sciml.ai/dev/features/linear_nonlinear/

nonlinear

initialization

1. OrdinaryDiffEq.alg_cache

• Val{true} implies it operates in-place

2. Calls build_nlsolver, dispatching on AbstractNLSolverAlgorithm type: these are defined in DiffEqBase. e.g. NLNewton, which would call build_nlsolver.

   a. nf = nlsolve_f(f, alg), which extracts the implicit part if using a SplitFunction

   b. SciMLBase.islinear(f): this appears to only be true if f (or f.f1) is an AbstractDiffEqOperator that isconstant?

   • If so, then it sets uf and jac_config to nothing, and calls linsolve(Val{:init},nf,u). Otherwise...

   • uf = build_uf(alg,nf,t,p,Val(true)): returns a SciMLBase.UJacobianWrapper, which captures the function, parameter and time, and can be called with just (dest, src) args (to make it easier to use the AD tools I assume?).

   • jac_config = build_jac_config(alg,nf,uf,du1,uprev,u,ztmp,dz): I think it allocates the cache for forming the Jacobian matrix?

   c. build_J_W:

   • islinearfunction(f, alg): "return the tuple (is_linear_wrt_odealg, islinearodefunction).": 2nd is true if islinear(f), 1st if 2nd or alg isa SplitAlgorithms and islinear(f.f1.f)?

   • One of:

     • if f.jac_prototype isa DiffEqBase.AbstractDiffEqLinearOperator: W = WOperator{IIP}(f, u, dt), J = W.J

     • if f.jac_prototype !== nothing: J = similar(f.jac_prototype); W = similar(J)

     • if islin: unwrap f if split; wrap in DiffEqArrayOperator to get J; W = WOperator{IIP}(f.mass_matrix, dt, J, u)

     • otherwise J = ArrayInterface.zeromatrix(u); W = similar(J).

   • Returns J, W

   d. nlcache = NLNewtonConstantCache object

   e. returns NLSolver object.

running

1. perform_step!:

2. nlsolve!

a. update_W!

- [`calc_W!`](https://github.com/SciML/OrdinaryDiffEq.jl/blob/v5.59.1/src/derivative_utils.jl#L422) with `transform=true`

  - if `DiffEqBase.has_Wfact_t(f)`

    - `set_W_γdt!(nlsolver, dtgamma)`

  - otherwise:

    - `do_newJW`

    - if `W isa WOperator`:

      - `DiffEqBase.update_coefficients!(W,uprev,p,t)`

    - otherwise

      - `calc_J!(J, integrator, lcache)` if `J` changed in `do_newJW`

      - `update_coefficients!(mass_matrix,uprev,p,t)`

      - `jacobian2W!(W, mass_matrix, dtgamma, J, W_transform)` if `W` changed in `do_newJW`

      - `set_new_W!(nlsolver, new_W)`

      - `set_W_γdt!(nlsolver, dtgamma)` if `W` changed in `do_newJW`