Skip to content

Implicit solvers

Simon Byrne edited this page Jul 2, 2021 · 15 revisions

Some very limited docs of the linear solvers.

Similarly, some docs for AbstractDiffEqOperators.

Linear

e.g. as used by Rosenbrock23:

Initialization

alg_cache

  • J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  • tf = TimeGradientWrapper(f,uprev,p)
  • uf = UJacobianWrapper(f,t,p)
  • alg.linsolve(Val{:init},uf,u)
  • build_grad_config(alg,f,tf,du1,t)
  • build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2,Val(false))

Running

perform_step!

  • new_W = calc_rosenbrock_differentiation!(integrator, cache, γ, γ, repeat_step, false)
  • calculate_residuals!
  • cache.linsolve(vec(k₁), W, vec(linsolve_tmp), new_W, Pl=DiffEqBase.ScaleVector(weight, true), Pr=DiffEqBase.ScaleVector(weight, false), eltol=opts.reltol)

Running

Nonlinear

e.g. as used by SBDF (which includes IMEXEuler)

Initialization

  1. OrdinaryDiffEq.alg_cache
  • Val{true} implies it operates in-place
  1. 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?

    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

        • WOperator:

          A linear operator that represents the W matrix of an ODEProblem, defined as

          $$W = MM - \\gamma J$$

          or, if transform=true:

          $$W = \\frac{1}{\\gamma}MM - 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)

      • docstring references research note
      • evaluates f
      • if W isa DiffEqBase.AbstractDiffEqLinearOperator, calls update_coefficients!(W, ustep, p, tstep)
      • linsolve(vec(dz), W, b, iter == 1 && new_W; Pl=DiffEqBase.ScaleVector(weight, true), Pr=DiffEqBase.ScaleVector(weight, false), reltol=reltol)
      • checks residuals
    • check divergence

    • apply_step!(nlsolver, integrator)

    • check convergence

  • postamble!(nlsolver, integrator)

Clone this wiki locally