-
Notifications
You must be signed in to change notification settings - Fork 17
Implicit solvers
Some very limited docs of the linear solvers.
Similarly, some docs for AbstractDiffEqOperators.
e.g. as used by Rosenbrock23:
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))
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)
e.g. as used by SBDF (which includes IMEXEuler)
-
Val{true}implies it operates in-place
-
Calls
build_nlsolver, dispatching onAbstractNLSolverAlgorithmtype: these are defined in DiffEqBase. e.g.NLNewton, which would callbuild_nlsolver.a.
nf = nlsolve_f(f, alg), which extracts the implicit part if using aSplitFunctionb.
SciMLBase.islinear(f): this appears to only betrueiff(orf.f1) is anAbstractDiffEqOperatorthatisconstant?-
isconstant(::DiffEqArrayOperator)if it uses the defaultupdate_func -
If so, then it sets
ufandjac_configtonothing, and callslinsolve(Val{:init},nf,u). Otherwise... -
uf = build_uf(alg,nf,t,p,Val(true)): returns aSciMLBase.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 ifislinear(f), 1st if 2nd oralg isa SplitAlgorithmsandislinear(f.f1.f)? -
One of:
-
if
f.jac_prototype isa DiffEqBase.AbstractDiffEqLinearOperator:W = WOperator{IIP}(f, u, dt),J = W.J-
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: unwrapfif split; wrap inDiffEqArrayOperatorto getJ;W = WOperator{IIP}(f.mass_matrix, dt, J, u) -
otherwise
J = ArrayInterface.zeromatrix(u);W = similar(J).
-
-
Returns
J, W
d.
nlcache = NLNewtonConstantCacheobjecte. returns
NLSolverobject. -
-
calc_W!withtransform=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)ifJchanged indo_newJW update_coefficients!(mass_matrix,uprev,p,t)-
jacobian2W!(W, mass_matrix, dtgamma, J, W_transform)ifWchanged indo_newJW set_new_W!(nlsolver, new_W)-
set_W_γdt!(nlsolver, dtgamma)ifWchanged indo_newJW
-
- if
-
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, callsupdate_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)