Extensions

HDF5 files of Runes2P: the output of a TimeSolution, the checkpoints of a run and their readers. Loaded when HDF5 is, so that Runes2P does not depend on HDF5.

Interpolators of Interpolations.jl as radial profiles: evaluated on the points of the radial grid, and recorded by their points, the values there, the method and the rule beyond the last point, from which the run can be rebuilt.

Intel MKL PARDISO as the linear solver MKLPardiso of the Newton iterations, through Pardiso.jl.

MAX_RELATIVE_RESIDUAL = 1e-6

The largest relative residual $\|b - A x\|/\|b\|$ of a solve of MKLPardiso before the Woodbury update: a solve above it returns :singular_matrix, as SparseLU does on a zero pivot, so that the Newton iteration ends and the stepper retries the step shorter. PARDISO perturbs small pivots instead of failing, so a singular or nearly singular matrix gives a finite solution with a large residual. 1e-6 is a choice of this package (see the docs page Choices of this package): three orders of magnitude above the residuals of PARDISO on the disruption case, $10^{-13}$ to $10^{-9}$, and four or more below those of singular matrices, $10^{-2}$ to $10$.

SciML form of an EquationSystem: loaded when SciMLBase is, so that OrdinaryDiffEq and NonlinearSolve can drive Runes2P systems without making SciML a dependency of Runes2P.

ode_function(system::EquationSystem)

SciMLBase.ODEFunction of M du/dt = rhs(u, t): in-place right-hand side f(du, u, p, t) with the time t [s], the mass matrix, and the Jacobian-vector product jvp(Jv, v, u, p, t) (stored sparse Jacobian plus matrix-free terms). The sparse Jacobian itself (jac, jac_prototype) is included only when it is the whole Jacobian, i.e. when the system has no matrix-free and no explicit terms, and no term stores an approximate block (a frozen limiter, or the Jacobian of another scheme). In the ODE form explicit terms are part of f but not of the Jacobian, which then only approximates ∂f/∂u (as DREAM's explicit knock-on term does). A system with a term that lags part of its residual to the start of the step (lags_in_ode_form) gets neither jac nor jvp: in the ODE form the start of the step is the state itself, and the stored Jacobian leaves out the derivative of the lagged part. The solver then takes its own derivatives. The parameter p is ignored.

Cost: jvp assembles the whole stored Jacobian at u on every call (one AD pass per color, about 23 operator applications for the Fokker-Planck stencil) before the product. Fine for occasional products; a Krylov solver driven through SciML should reuse an assembled Jacobian instead.