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!

• update_W!

• calc_W! 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

• initialize!(nlsolver, integrator)

• if get_new_W!(nlsolver)

  • initial_η(nlsolver, integrator)

• for iter in 1:maxiters

  • compute_step!(nlsolver, integrator)
  • check divergence
  • apply_step!(nlsolver, integrator)
  • check convergence

• postamble!(nlsolver, integrator)