Linear solvers and Newton

LinearSolver

Supertype of the linear solvers used by NewtonSolver. Each Newton iteration solves A x = b with the scaled Newton matrix A = W (M/dt - J) S (row weights W, unknown scales S, see NewtonSolver), where J is the stored Jacobian of the implicit :sparse terms plus, for solvers that support them, the implicit matrix-free terms. A solver implements

• linear_solver_cache(solver, matrix, system, row_weights, scales): the state reused between solves; matrix is the scaled stored matrix (its pattern is fixed, its values change), row_weights and scales are the vectors W and S, which the Newton iteration updates in place;
• solve!(x, cache, matrix, b, u, t, dt, iteration): solve with the matrix at the state u, time t [s] and step dt [s], in Newton iteration iteration of the step; returns :ok, :singular_matrix or :linear_solve_failed (an iterative solve that did not reach its tolerance);
• requires_stored_jacobian(solver): whether every implicit term must be in the stored matrix, apart from matrix-free terms with a low-rank Jacobian (has_low_rank_jacobian), which the solver adds by a Woodbury update;
• krylov_iterations(cache), preconditioner_updates(cache): Krylov iterations and preconditioner factorizations of the last solve (0 for direct solvers);
• factorization_time(cache): seconds spent factorizing in the last solve (the LU of a direct solver, the preconditioner blocks of a Krylov solver);
• linear_tolerance(cache): bound on the relative residual ‖b - A x‖/‖b‖ of the last solve, used by the residual-floor rule of the Newton iteration: 0 for SparseLU, the measured relative residual for MKLPardiso, its relative_tolerance for GMRES;
• newton_operator_exact(cache): whether the solver's Newton operator is the exact Jacobian when the system has terms with an approximate stored block (jacobian_kind; false for solvers on the stored matrix alone);
• matrix_unchanged!(cache, unchanged), optional: told before every solve whether the matrix has the values of the last solve, so that a factorization of it can be kept (a system linear in the state, NewtonSolver). Without a method the solver treats every matrix as new.

MKLPardiso(; threads=1, refinement_steps=0)

Monolithic sparse direct solve by Intel MKL PARDISO, the solver of DREAM's LINEAR_SOLVER_MKL, through the package Pardiso.jl. Its solves come with the extension Runes2PMKLPardisoExt, loaded by using Pardiso; without it, a NewtonWorkspace with this solver is refused with an ArgumentError. It takes the systems that SparseLU takes, the field-particle term (low_rank_jacobian) by the same Woodbury update. Its factorization has another pivot order than UMFPACK's, so the converged states differ from those of SparseLU within the Newton tolerance, and the Newton iterations of a step can differ.

PARDISO factorizes the real unsymmetric matrix with the defaults of its pardisoinit: a nested-dissection ordering (METIS), a weighted matching with scaling, and pivots below $10^{-13}$ times the infinity norm of the scaled matrix perturbed. The ordering and the matching need the values of the matrix: they are computed at the first solve, from its matrix, and kept for every later matrix of the same pattern until a solve fails (below); their time counts in the factorization time of that solve (SolverStatistics). A run resumed from a checkpoint computes them from its own first matrix, so it repeats the uninterrupted run as with another pivot order, not bit for bit: within the Newton tolerance, and within round-off where the Newton iterations of the steps are the same. The numeric factorization is made in every solve, except when the matrix is that of the last solve (a system linear in the state, NewtonSolver). PARDISO reads the compressed columns of the matrix as the compressed rows of its transpose and solves the transposed system, so no matrix is converted. Its integers have 32 bits, so the factors hold at most $2^{31} - 1$ nonzeros; a larger factorization ends the run with an error that names the limit. The factors are in memory that MKL allocates and Julia's garbage collector does not count: they are released when the workspace is collected, so repeated runs in one process can hold several factorizations until GC.gc().

• threads: the threads of PARDISO, MKL's own; Julia's threads and the BLAS threads do not apply. The count is a setting of MKL for the whole process, set before every solve, so one workspace solves one task at a time, and two workspaces with different counts do not solve at the same time. 1 is DREAM's value (fvm/Solvers/MIMKL.cpp), and the results of a run are then reproducible, which MKL does not promise with more threads.
• refinement_steps: the largest number of steps of iterative refinement in each solve (PARDISO's iparm(8)). 0 is a choice of this package (see the docs page Choices of this package): on the disruption case ten steps take longer than the factorization, and the solve without them reaches a relative residual of about $10^{-9}$. With 0, PARDISO still refines twice after a perturbed pivot. DREAM sets 10 (fvm/Solvers/MIMKL.cpp).

Every solve measures its relative residual $\|b - A x\|/\|b\|$ with one product by the matrix, before the Woodbury update, with refinement too, since the refinement stops at its cap and promises no residual. That residual is the bound that the residual-floor rule of the NewtonSolver compares with its floor_linear_tolerance: a perturbed pivot that spoils the solve within the bound below then turns both waivers of the residual test off, where SparseLU reports 0.

PARDISO perturbs a small pivot instead of failing, so a singular or nearly singular matrix gives a finite solution with a large residual. A solve whose relative residual exceeds MAX_RELATIVE_RESIDUAL of the extension, $10^{-6}$ (a choice of this package, see the docs page Choices of this package), or is not a number, fails with :singular_matrix, as a zero pivot does with SparseLU: the Newton iteration ends at that solve, and the time stepper retries the step shorter. MKL's errors -4 (a zero pivot, or a failed numeric factorization or iterative refinement) and -7 (a singular diagonal) also return :singular_matrix; its other errors throw. A solve that returns :singular_matrix drops the ordering and the matching with the factorization, and the next factorization computes them again from its own matrix, since those of another matrix of the same pattern can spoil its solves.

SparseLU()

Monolithic sparse direct solve by LU factorization (UMFPACK, through SparseArrays), the DREAM-like baseline. The symbolic analysis of the sparsity pattern is done once and reused for every later matrix with the same pattern. Needs every implicit term in the stored Jacobian, except matrix-free terms with a low-rank Jacobian (low_rank_jacobian, the field-particle term), which a Woodbury update of the factorization solves exactly. A system with other implicit matrix-free terms (the implicit knock-on term) is refused; those need the Newton-Krylov solver GMRES.

GMRES(; preconditioner=FluidSchurLU(), relative_tolerance=1e-10, max_iterations=200,
      restart_length=50, preconditioner_refresh=:every_iteration, refresh_step_change=0.5,
      refresh_count_growth=2.0, operator=:exact)

Newton-Krylov linear solver: restarted GMRES (Krylov.jl) on the scaled Newton operator W (M/dt - J) S, where J is the stored Jacobian plus the implicit matrix-free terms, right preconditioned by preconditioner: FluidSchurLU by default, which holds the couplings between the radii, or a RadiusBlockLU or RadiusBlockILU of the stored matrix. With right preconditioning the residual GMRES minimizes is the true residual of the scaled system; a solve succeeds when ‖b - A x‖ ≤ relative_tolerance ‖b‖ within max_iterations iterations, and otherwise the Newton iteration fails with :linear_solve_failed.

operator sets what J is for terms whose stored block is an approximation of their Jacobian (jacobian_kind): :exact adds their jacobian_correction_product! to every product, so Newton is exact and the approximate block only preconditions; :stored uses the stored matrix as it is, a quasi-Newton iteration (DREAM's AD_INTERP_JACOBIAN_LINEAR and UPWIND). Without such terms the two are the same.

preconditioner_refresh sets when the preconditioner is factorized again:

• :every_iteration: every Newton iteration, from the current matrix. For a system that is linear in the state the matrix of the later iterations of a step is that of the first, and its factors are kept, so the preconditioner is factorized once per step (NewtonSolver, reuse_linear_matrix);
• :every_step: in the first Newton iteration of each time step;
• :on_demand: when the time step differs from that of the factorization by more than the relative refresh_step_change, or when the last solve needed more than refresh_count_growth times the Krylov iterations of the first solve after the factorization. The Newton matrix W (M/dt - J) S changes with the step through the mass term and through the row weights W of a kinetic unknown, which follow its entries at the start of the step (NewtonSolver); the factors of a nearby step precondition it well, and the count rule catches coefficients and entries that drift in time and the factors of a longer step applied at a shorter one. A halved step is a change of exactly 0.5 and keeps the factors with the default; a doubled step refreshes. The stored matrix of a system with prescribed coefficients, a constant time step and no kinetic unknown is the same in every step, so one factorization then serves the whole run. Under a step controller that grows the step at its cap of 2 per step, the default refresh_step_change factorizes at every such step, and :on_demand pays off where the step changes by less.

A solve that fails with a preconditioner from an earlier iteration is repeated once with a fresh one.

low_rank_update adds the low-rank Jacobians of implicit matrix-free terms (low_rank_jacobian, the field-particle term) to the blocks they lie in by a Woodbury update: k block solves per factorization and O(n k) work per application for rank k. Without it the preconditioner leaves those terms out, as it leaves out the knock-on term.

On hot-tail grids of 270 × 15 cells (fully kinetic) and of 125 × 15 cells (superthermal), with radii coupled by radial transport:

• RadiusBlockLU leaves out the radial coupling, and the Krylov count grows with the number of radii. restart_length has to exceed the count: below it the iteration stagnates and the step ends with :linear_solve_failed.
• FluidSchurLU holds the radial coupling by line solves and the fluid unknowns exactly. Its count grows more slowly with the number of radii. The default restart length is 50 (a choice of this package; see the docs page Choices of this package).
• relative_tolerance stays at or below the floor_linear_tolerance of the NewtonSolver (1e-8). With 1e-6 the Newton iteration cannot accept a residual at its round-off floor and ends with :max_iterations; the Newton workspace of a larger tolerance warns. A tolerance of 1e-8 takes fewer Krylov iterations per Newton iteration than 1e-10, and more Newton iterations per step.
• The memory is that of the block factors and of the Krylov basis, restart_length vectors.

RadiusBlockILU(; drop_tolerance=1e-4)

BlockPreconditioner that computes an incomplete LU factorization of each block (Crout ILU of IncompleteLU.jl), which drops entries of the factors below drop_tolerance in magnitude. Couplings between blocks (radial transport, scalar unknowns) and the matrix-free terms are left out.

RadiusBlockLU()

BlockPreconditioner that factorizes each block exactly (UMFPACK, symbolic analysis reused). Couplings between blocks (radial transport, scalar unknowns) and the matrix-free terms are left out.

BlockPreconditioner

Preconditioners for GMRES: the scaled stored Newton matrix restricted to the blocks of radius_block_range (all unknowns at one radius) and scalar_block_range, each block factorized on its own. The subtypes are RadiusBlockLU, RadiusBlockILU and FluidSchurLU.

FluidSchurLU(; radial_lines=:auto)

The default preconditioner of GMRES, for systems whose radii are coupled through the fluid unknowns (Ampère's law, the plasma current and the other scalars) and through the radial transport of the distributions. RadiusBlockLU leaves both couplings out, and its Krylov count grows with the number of radii.

• The kinetic unknowns of each radius form one block, factorized exactly (UMFPACK), with the Woodbury update of low-rank terms as in RadiusBlockLU.
• The fluid unknowns of all radii (every radial and scalar unknown) form one block. With the kinetic blocks eliminated, its Schur complement S = A_ff − Σ_r A_fk,r A_kk,r⁻¹ A_kf,r is a dense matrix of the size of the fluid unknowns (about 11 per radius and 3 scalars in a self-consistent run), factorized exactly. The columns A_kf,r (the electric field, the magnitude of the particle source) are few per radius. This stage is the exact inverse of the stored matrix without the kinetic couplings between radii.
• The radial lines are the tridiagonal systems along the radii of each kinetic cell: the diagonal of the stored matrix and the couplings of the cell to the same cell at the two neighbouring radii. The application is symmetric, lines, Schur stage, lines, each stage on the residual of the ones before: one block solve, two line solves and two products with the stored matrix and the low-rank terms of the blocks. The lines need one pitch grid for all radii (the same restriction as the four-point schemes in r). Without stored kinetic couplings between radii (no radial transport) they are left out. radial_lines = :auto (default) takes them where every kinetic unknown has one pitch grid at all radii and leaves them out otherwise, with a warning: the fluid coupling stays exact, and the radial transport of the distributions is left to the Krylov iteration. true refuses a grid whose pitch faces differ between radii, false never takes them.

The preconditioner is a fixed linear map for each linear solve. With preconditioner_refresh = :on_demand the factors can be older than the stored matrix, whose current values the residuals between the stages take: still a fixed linear map for each solve, and a preconditioner of the current matrix. A kinetic row must not couple to a radial unknown of another radius; GMRES refuses such a system.

On the self-consistent system with the values of reference/dream/cases/self_consistent_field.py (the hot-tail and the runaway grid, the electric field as an unknown, minor radius 0.22 m), at 6 to 24 radii and optionally with radial diffusion: without radial transport the Schur stage alone is exact and GMRES takes 1 to 2 iterations at 6 to 24 radii, against 8 to 22 with RadiusBlockLU. With radial diffusion at a coupling typical of a disruption run (D Δt/Δr² ≈ 5) the symmetric application takes 5 iterations at 12 and 24 radii, against 49 and 55.

POLISHED_SCHEMES

The schemes of the advection whose kinetic unknowns polish_negatives = :auto of NewtonSolver polishes, a choice of this package (see the docs page Choices of this package): :upwind and :exponential_fitting, whose step is an M-matrix as long as the central cross diffusion keeps the off-diagonal couplings non-negative (docs page Positivity), and the limiters :tcdf, :ospre and :van_albada, which stayed non-negative in the tests without such a guarantee.

POLISH_MISSES

The number of polished Newton solves in a row whose extra iteration leaves negative entries of an unknown, after which polish_negatives = :auto of NewtonSolver stops polishing it; a choice of this package (see the docs page Choices of this package). An unconverged limiter is cleaned by one extra iteration, so a few misses in a row show a scheme whose converged solution is negative, and five bound the cost of finding out.

NewtonResult

Outcome of newton_solve!: status (:converged, or the reason of failure: :max_iterations, :singular_matrix, :non_finite_residual, :non_finite_jacobian (a non-finite entry of the stored Jacobian at a finite residual), :non_finite_update, :domain_error, :line_search_failed, :linear_solve_failed: an iterative linear solve did not reach its tolerance, :below_lower_limit: an entry at its lower limit would be pushed below it, or the converged state lies below it), the number of iterations, the scaled residual norm before the first iteration and at the end (which may exceed residual_tolerance on a stiff step, see NewtonSolver), the scaled norm max |δ|/scale of the last update, the estimated convergence_order (NaN when fewer than three updates are above round-off), the smallest step length used, the Krylov iterations and preconditioner factorizations summed over the iterations (0 for direct solvers), residual_floor, true when the step converged with the residual above residual_tolerance: every row above it was within its round-off bound, or the residual had stopped decreasing (or, for a linear step, it was not checked against it; see NewtonSolver), and linear_step, true when the step was accepted as a linear system, on the update test alone (at the second iteration, see NewtonSolver), and polished, true when the iteration took one more iteration for negative entries of a kinetic unknown (polish_negatives of NewtonSolver).

NewtonSolver(; linear_solver=SparseLU(), max_iterations=30, max_reduction=0.9,
             residual_tolerance=1e-8, floor_linear_tolerance=1e-8, line_search=false,
             sufficient_decrease=1e-4, max_backtracks=20, linear_acceptance=true,
             reuse_linear_matrix=true, polish_negatives=:auto)

Settings of the Newton iteration of a backward-Euler step.

• Convergence needs both a small update and a small residual, since an iteration can stall with a vanishing update at a large residual: every entry of the full Newton update δ satisfies |δ| ≤ absolute + relative |u| with the unknown's Newton Tolerance, and the residual, measured as the change it causes in one step relative to the unknown's scale (dt |F|/scale on differential rows, |F|/scale on algebraic rows), is at most residual_tolerance everywhere. On the differential rows of a kinetic unknown the scale is the larger of the row scale and the entry whose time derivative the row holds, at the start of the step (|u_old|): the collision fluxes of the bulk cancel to a net term orders of magnitude smaller, and the round-off of that cancellation grows with the entry, which exceeds the scale of the unknown by orders of magnitude in the bulk. The tail, below the scale, keeps the row scale. In a bulk cell the residual test is then looser by |u_old|/scale in absolute terms, as the update test is relative to |u| there. A bulk that grows by orders of magnitude within one step can still meet the round-off floor of a stiff step, the exception of the next item. An infinite absolute tolerance excludes an unknown from both tests; the merit function of the line search keeps all rows. The unknown of a Constraint is a rate: the system fixes the change dt u that it causes over the step, and round-off of the constrained quantity divided by dt is the accuracy it can reach. Its update test is |δ| ≤ absolute + relative |u| + relative scale/dt, with the scale of its equation, which bounds the error of the change over the step as the tolerances bound the other unknowns (such as the particle source at steps below 1e-6 s, where the round-off of the density divided by dt exceeds a tolerance on the source itself).

• A stiff step is an exception to the residual test. When the terms of a row cancel to many digits, round-off holds its residual above residual_tolerance whatever the iterate: in a bulk cell of a kinetic unknown the fluxes of the Fokker-Planck stencil through the cell exceed their net term, which the particle source balances, by many orders of magnitude, and more so for a long step. The waiver of the residual floor then accepts an iterate whose update test passes when one of two tests holds:

  • The round-off bound. Every row whose scaled residual exceeds residual_tolerance is within its round-off bound, |F_i| ≤ c eps Σ_j |(∂F/∂u)_ij u_j|, with ∂F/∂u = M/dt - J from the stored Jacobian J of the iteration. The terms of the row carry rounding errors of about eps times their sizes, so no iterate brings the row lower; |(∂F/∂u)_ij u_j| estimates the size of a term by its linearization, exact for a term linear in u_j and off by a factor of order one or more for a strongly nonlinear one, at the level of eps either way. The explicit terms, which are not in J, and the start state of the transient are left out, so the bound errs low and the stall below then serves. The factor c = 1 is a choice of this package (see the docs page Choices of this package). The step is accepted at the first iteration that reaches the floor. The bound needs every implicit term in the stored Jacobian, exactly: it is not used with a matrix-free term or an approximate stored block (jacobian_kind).
  • The stall. The residual stops decreasing (it grew, or fell by less than 10%, in the last iteration). It serves where the bound is not available or does not hold; the residual at the floor wanders, so the stall can take a few iterations more than the bound.

  Such a step is flagged (residual_floor in NewtonResult, accepted_at_residual_floor in SolverStatistics): the residual is above residual_tolerance and only the update test bounds the error. The same waiver would also accept a linearly converging iteration whose small update understates its error (an inexact Jacobian or linear solve). An inexact linear solve with relative residual η contracts the error by about η per iteration, so both of its tests apply only when the last linear solve was accurate to floor_linear_tolerance (always true for SparseLU; for MKLPardiso its measured relative residual; for GMRES its relative_tolerance, and the workspace warns when that is above floor_linear_tolerance), and never when the Newton matrix is itself approximate (a term with jacobian_kind(term) === :approximate solved on the stored matrix): such a quasi-Newton iteration converges linearly by design, its residual falls by a constant factor per iteration that the stall would mistake for the floor, and only the residual criterion holds it.

• A system that is linear in the state (every implicit term is_linear) with an exact Newton matrix and a linear solve accurate to floor_linear_tolerance is accepted as soon as the update test passes, without the residual test (linear_acceptance, flagged as linear_step in NewtonResult): the first iteration solves the step, the second confirms it with an update at round-off, and the residual is then at its round-off floor whatever it is, so the iterations that would only establish the floor are saved. The second iteration is kept, since a solve bounds the residual and not the error: accepting the first solve alone loses its refinement on ill-conditioned steps. A step shortened by a lower limit or the line search is not accepted this way.

• The stored Jacobian of such a system does not depend on the Newton iterate, only on the coefficients, which are functions of the time, and on the state at the start of the step (the column of a lagged knock-on term, KnockOnTerm with :explicit_primary). With reuse_linear_matrix it is assembled once per step and kept, with the Newton matrix and its factorization, for the later iterations of the step and for a step repeated to the same end time from the same start. The key is the end time and the length of the step, which identify the start state for every sequence of steps integrate produces: the step doubling of a step controller solves to one end time from two start states with two step lengths. An update of the setup by the time stepper discards them (invalidate_jacobian!). reuse_linear_matrix = false assembles and factorizes in every iteration, as DREAM does.

• Lower limits (LowerLimit, declared on the unknowns) shorten the whole step so that no limited entry loses more than the fraction max_reduction of its distance to its limit in one iteration. Only steps towards the limit are shortened, never growth. max_reduction = 0.9 is DREAM's value (MaximalPhysicalStepLength). The temperature unknown of a PlasmaFromUnknowns has the floor of that source as its limit, so an update that would cross the floor is damped and a step whose solution lies below the floor is refused. A temperature that only fluid terms read, with no declared limit, takes the default floor of 0.5 eV (fluid_temperature_floors). DREAM damps towards zero and has no floor.

• polish_negatives takes one more iteration when an iterate that passes both tests has negative entries in a polished kinetic unknown, and accepts the next iterate that passes them; at most once per step, and never on a linear step. A scheme whose converged discrete solution stays non-negative leaves negative entries near zero in the tail as long as Newton has not converged it there; the residual test does not see them, since their size is far below the tolerance of the unknown, and the extra iteration removes them. A scheme whose converged solution has negative entries keeps them and pays the iteration.

  • :auto (the default, a choice of this package; see the docs page Choices of this package) polishes a kinetic unknown whose equation has Fokker-Planck terms that all use, in p and ξ, a scheme of POLISHED_SCHEMES. When the extra iteration leaves negative entries of an unknown in POLISH_MISSES polished solves in a row, its schemes do not keep that run non-negative: :auto stops polishing the unknown for the rest of the run (a resumed segment starts again) and writes a warning to the log of the run, or to the logger of Julia outside integrate.
  • true polishes every kinetic unknown, whatever its scheme, and never stops.
  • false never polishes, as DREAM's nonlinear solver.

• line_search = true adds backtracking on ½‖F‖² of the scaled residual (Armijo condition with sufficient_decrease, at most max_backtracks reductions), starting from the step length the lower limits allow. Off by default; DREAM's nonlinear solver (src/Solver/SolverNonLinear.cpp) has no line search (see the docs page Comparison with DREAM).

The linear systems are solved in scaled form: rows weighted as the residual above, columns multiplied by the unknowns' scales. The default linear_solver = SparseLU() (UMFPACK, SparseLU) is a choice of this package (see the docs page Choices of this package); MKLPardiso and GMRES are the others.

NewtonWorkspace(system, solver::NewtonSolver)

Buffers, tolerances, lower limits and the linear-solver cache for Newton iterations on system, reused over all time steps. The per-iteration history of the last solve is in the fields residual_norms, update_norms and step_lengths. exact_jacobian records whether the Newton matrix is the exact Jacobian (no term with an approximate stored block, or a linear solver whose operator restores it, GMRES(operator=:exact)), and linear whether the system is linear in the state and the matrix exact, so the first iteration solves a step and the second confirms it. For such a system assembled_time and assembled_step are the end time [s] and the length [s] of the step whose stored Jacobian and Newton matrix are kept (NaN: none kept, see invalidate_jacobian!); jacobian_assemblies counts the assemblies of the stored Jacobian since the workspace was made.

invalidate_jacobian!(workspace::NewtonWorkspace)

Discard the stored Jacobian and the Newton matrix that a system linear in the state keeps between the iterations of a step (NewtonSolver, reuse_linear_matrix). To be called when coefficients of the system were changed by hand at an unchanged time. The time stepper calls it after an update that changed the setup.

newton_solve!(u, u_old, workspace::NewtonWorkspace, t, dt; reference=u_old)

Solve the backward-Euler step from u_old to time t [s] with step dt [s], starting from the guess in u and overwriting it. reference replaces u_old in the time derivative alone (residual!): BDF2 passes its combination of the two last states with its effective step. Returns a NewtonResult; on failure u holds the last iterate. Never throws for a failed iteration: a singular matrix, a non-finite residual, Jacobian or update, a DomainError from a term (evaluation outside the physical domain), a linear solve that did not converge, a failed line search, an entry below its lower limit at convergence and the iteration limit are reported through the status, for the time stepper to handle.