-
Notifications
You must be signed in to change notification settings - Fork 9
Implicit solvers
Some very limited docs of the linear solvers.
Similarly, some docs for AbstractDiffEqOperator
s.
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 onAbstractNLSolverAlgorithm
type: 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 aSplitFunction
b.
SciMLBase.islinear(f)
: this appears to only betrue
iff
(orf.f1
) is anAbstractDiffEqOperator
thatisconstant
?-
isconstant(::DiffEqArrayOperator)
if it uses the defaultupdate_func
-
If so, then it sets
uf
andjac_config
tonothing
, 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 SplitAlgorithms
andislinear(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
: unwrapf
if split; wrap inDiffEqArrayOperator
to getJ
;W = WOperator{IIP}(f.mass_matrix, dt, J, u)
-
otherwise
J = ArrayInterface.zeromatrix(u)
;W = similar(J)
.
-
-
Returns
J, W
d.
nlcache = NLNewtonConstantCache
objecte. returns
NLSolver
object. -
-
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)
ifJ
changed indo_newJW
update_coefficients!(mass_matrix,uprev,p,t)
-
jacobian2W!(W, mass_matrix, dtgamma, J, W_transform)
ifW
changed indo_newJW
set_new_W!(nlsolver, new_W)
-
set_W_γdt!(nlsolver, dtgamma)
ifW
changed 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)