Time stepping

FixedSteps(; max_halvings=10)

The stepper of DREAM's constant time stepper: each interval between two time points of the run is one step. When a step fails, it is retried with half the time step, up to max_halvings times; the shorter steps continue to the end of the interval, whose end point is always hit exactly, and the next interval starts again with its full length. A step that still fails throws TimeStepFailure.

RunOptions(; solver=NewtonSolver(), stepper=FixedSteps(), save_times=nothing,
           updates=(), checkpoint=nothing, resume=false, stop=nothing, log=nothing,
           description=Dict())

The options of a run of integrate, which takes them as keywords or as one record. save_times = nothing saves at the time points. The others are described with integrate.

SolverStatistics

Solver statistics of a run, one entry per attempted step (accepted or not), with its own time array:

• time, time_step: end time and length of the attempted step;
• accepted, status: whether the step was accepted, and the NewtonResult status; a step that converged and was not accepted was rejected by the stepper on its error estimate, or was a half step of a step doubling;
• kind: :step for a step of the run (accepted, failed or rejected), :half_step for the two half steps of a step doubling (StepDoubling);
• error: the error estimate of the step relative to its tolerance (accepted at most 1), NaN where the stepper made none;
• newton_iterations, krylov_iterations, preconditioner_updates (factorizations of the Krylov preconditioner);
• initial_residual_norm, residual_norm, update_norm: scaled norms before the first and after the last Newton iteration (see NewtonSolver);
• convergence_order: estimated order of the Newton iteration (about 2 with a correct Jacobian; about 1 flags an inexact Jacobian, with which the iteration still converges);
• smallest_step_length: the shortest Newton step length used (1 when neither the lower limits nor the line search shortened a step);
• accepted_at_residual_floor: the step converged with its residual above residual_tolerance, at the round-off floor of a stiff step: within the round-off bound of its rows, or where the residual had stopped decreasing (see NewtonSolver); such a step costs no iteration beyond the floor when the bound applies, and a few when only the stall does;
• accepted_as_linear: the step was accepted as a linear system, on the update test alone at the second Newton iteration (see NewtonSolver);
• polished_for_negatives: the step converged after one more Newton iteration for negative entries of a kinetic unknown (polish_negatives of NewtonSolver);
• wall_time: seconds spent in the step, of which jacobian_time in the assembly of the stored Jacobian blocks, residual_time in residual evaluations (line search included), linear_solve_time in the linear solves and, within it, factorization_time in the LU or preconditioner factorizations (the rest of the solve time is the Krylov iteration, with the matrix-free terms' products).

The history of every Newton iteration is in iteration_attempt (index of the attempted step it belongs to), iteration_residual_norm, iteration_update_norm and iteration_step_length.

StepController(; initial_step=nothing, relative_tolerance=1e-3, max_growth=2.0, min_shrink=0.2,
               safety=0.9, gains=nothing, max_failures=5, failure_factor=0.2,
               max_rejections=20, tail_cut=1.0, tail_floor=1e-12, order=1)

The step controller of Runes2P: an error estimate of the backward-Euler step without an extra solve, tolerances per unknown with a floor, and a PI filter on the step length with a cap on its growth. It is the adaptive stepper of this package; DREAM's StepDoubling is the compatibility option.

The error estimate is the difference between the solution of the step and its predictor, the linear extrapolation through the two last accepted states. Two backward-Euler states satisfy u_n - u_{n-1} = h_old f(u_n), so the extrapolation u_n + (h / h_old) (u_n - u_{n-1}) is the forward-Euler predictor u_n + h f(u_n) with the derivative of the last step, whatever h_old is, and the implicit solution differs from it by h (f(u_{n+1}) - f(u_n)) = h² u'', twice the local error h² u'' / 2 of the step. The estimate is half the difference and costs one vector operation. The first step of a run has no history: it is solved as one step doubling (StepDoubling, twice the cost of that step alone), whose difference u_full - u_half is half the local error.

The derivation holds for an unknown whose own row holds its time derivative. The estimate is taken on every entry of every unknown, and an unknown of an Algebraic or a Constraint equation, or one whose row holds the time derivative of another unknown (TimeDerivative(of)), does not satisfy u_n - u_{n-1} = h_old f(u_n). On such an unknown the predictor is the linear extrapolation of its own two last values, and half the difference h (h + h_old) u''/4 (h² u''/2 at constant steps) measures the curvature of its path in time: the local error that backward Euler would make on a differential unknown with that path. It is not the error of its equation, which the Newton solve meets to its tolerance, nor the error that the differential unknowns carry into it. The doubling of a first step measures that carried error on an algebraic unknown: the difference of the values its equation gives at the end states of the full and the half steps. The unknown of a constraint is the mean rate over a step, and its doubling difference is the change of that rate between the full step and the second half step, of first order in the step (h S'/4 for a rate S, and the estimate twice that). With order = 2 the estimate of such an unknown measures the third derivative of its path. The unknown counts in the largest error and can set the step; a Tolerance with all parts zero in time_step_tolerances leaves it out.

The error of an unknown is the root mean square over its entries of the local error over absolute + relative |u|, with the Tolerance the system holds for the time-step error (time_step_tolerances of EquationSystem); an unknown without one gets relative_tolerance and the absolute tolerance relative_tolerance times its scale, a floor that keeps entries near zero (the tail of a distribution, a density before its source acts) from holding the step. A radial unknown that declares groups of its components (component_groups, the charge states of the ion species, also of a single species) gets no absolute part and the group part relative_tolerance in its place: each entry is measured against relative_tolerance times the sum over its group, the density of its species, so a trace impurity is controlled on its own scale. An unknown with all tolerances zero is left out. The error of the step is the largest over the unknowns; the step is accepted when it is at most 1.

The next step follows the PI step-size control of K. Gustafsson, M. Lundh and G. Söderlind, BIT 28, 270 (1988) and K. Gustafsson, ACM Trans. Math. Softw. 17, 533 (1991); see E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, 2nd ed. (Springer, 1996), section IV.2: h_new = h · safety · error^(-gains[1]) · error_old^(gains[2]), with the factor kept between min_shrink and max_growth. The default gains are the elementary 0.7 and 0.4 over the order of the local error, (0.35, 0.2) for backward Euler, whose local error is of second order in the step (DREAM's stabiliser has them over the order 1). The two gains are the fields gain_error and gain_memory of the record. The memory error_old is floored at 1e-4, as in Hairer's codes, so that a step after one with an error near zero is not shrunk by the limit. An interval update that reports a jump of the setup restarts the estimate: the next step takes it from a doubling, as the first step of a run. A change that is continuous in time, such as that of a RechesterRosenbluthUpdate whose level has no forced_time_points (any function other than a PiecewiseConstant), keeps the history. A rejected step is repeated with h · max(min_shrink, safety · error^(-1/(k+1))), k the order of the step, and keeps the memory of the filter; more than max_rejections rejections in a row end the run with TimeStepFailure (a settings error, such as a min_shrink near 1 with a step whose error does not fall). A step whose solve fails is retried with failure_factor times the length tried, at most max_failures times in a row.

The first step has the length initial_step (nothing: the first interval of the time points). Every step ends at or before the next time point and the next save time, so every saved state is a computed one. A step cut short by a time point or a save time is estimated like any other, and the next length follows from the step taken (at most max_growth times it), not from the length that was planned: after a run of close save times the step grows back by the cap, each step under the tolerance.

The cells of the hot tail of a kinetic unknown lie far below the absolute floor of the cell measure (relative × scale, 1e17 m⁻³ at a density of 1e20 m⁻³ and a tolerance of 1e-3, against 1e15 m⁻³ in the tail), so the cells alone leave the tail uncontrolled. With tail_cut (a momentum in units of mₑc, nothing for none) the error of every kinetic unknown also holds its moments above the cut: the density and the current density at every radius (density and current_density with a MomentumRegion), each estimated by the moment of the cell estimate, over absolute + relative |moment| with the absolute tolerance tail_floor times the density of the electrons on all kinetic grids at that radius for the density, and e c times that for the current, so that the floor does not depend on the scale of the unknown and stays a floor for a runaway grid whose electrons all lie above the cut. The error of the unknown is the larger of the root mean square over its cells and that over its 2 nr moments, taken per radius on the momentum grid of that radius. On a two-grid layout the runaway unknown lies above the cut as a whole, so its moments are its density and current.

With order = 2 the steps are those of the variable-step BDF2, u - a u_n + b u_{n-1} = c h f(u) with a = (1+ω)²/(1+2ω), b = ω²/(1+2ω), c = (1+ω)/(1+2ω) and ω the ratio of the step to the one before. The solve is the backward-Euler residual with the reference a u_n - b u_{n-1} in the time derivative and the effective step c h (residual!); the terms see u_n as with backward Euler, and the statistics record h. The estimate is Milne's: the quadratic extrapolation through the last three accepted states, ρ (u - P) with ρ = h(h+h1)/((h+h1+h2)(2h+h1) + h(h+h1)), 2/11 at constant steps, for a local error of third order in the step. The first two steps of a run and the first two after an update that jumps the setup are backward-Euler steps, with the gains scaled to their order. BDF2 is zero-stable for ω < 1 + √2, so max_growth must stay below that. It keeps neither the positivity of backward Euler with an M-matrix operator nor the second order for explicit and lagged terms, which it warns about. The default is order = 1, backward Euler, which DREAM uses.

The estimate assumes that every term is implicit. An explicit or lagged term contributes h h_old u''/2 in place of h² u''/2 to the difference, which the growth cap bounds to a factor of order 2.

The state of the controller (the step length, the memory of the filter, the last accepted state and its step) is in every checkpoint. With SparseLU a resumed run repeats the uninterrupted one bit for bit. With MKLPardiso the ordering and the matching of PARDISO come from the first matrix of the resumed run, so its states differ from those of the uninterrupted run as with another pivot order: within the Newton tolerance, and within round-off where the Newton iterations of the steps are the same. With GMRES(preconditioner_refresh = :on_demand) the factors of the preconditioner are not in the checkpoint, so the Krylov solves after a resume differ within the Newton tolerance. After either, the step lengths may differ from those of the uninterrupted run. A run resumed onto another layout (an unknown added or removed) keeps the step length and starts with a doubling, since the stored state of the predictor does not fit; the count of rejections in a row starts at zero in every segment.

StepDoubling(; initial_step=nothing, relative_tolerance=1e-2, check_every=0, safety=0.95,
             exponent=1.0, max_failures=5, failure_factor=0.2)

DREAM's adaptive time stepper (TYPE_ADAPTIVE, TimeStepperAdaptive.cpp), the compatibility option of the step control. A checked step is solved three times: two steps of half the length and one of the full length from the same state, so the error estimate costs two solves per step of progress. The error of an unknown is the 2-norm of the difference between the full step and the two half steps over the 2-norm of the full step, ‖u_full - u_half‖ / (absolute + relative ‖u_full‖), taken over every entry of the unknown at every radius; the error of the step is the largest over the unknowns. An unknown at zero with no absolute tolerance (a zero scale) stays out of that largest error but fails the step when it moves between the full and the half steps, as DREAM's acceptance dx ≤ epsa + epsr x does (ConvergenceChecker.cpp:193). The step is accepted when no such unknown fails it and its error is at most 1, and its solution is that of the full step, as in DREAM. The next step is safety · h / error^exponent in either case, without a cap on its growth (DREAM: safety 0.95, exponent 1.0; its stabiliser, exponent 0.7 with a memory of the last error, is compiled out). A step that such an unknown fails shows in the record as rejected, also with an error of at most 1. When it fails a step in which every other unknown has no error at all, DREAM's next step is unbounded; this step is retried with half its length instead (a choice of this package; see the docs page Choices of this package). check_every steps between two checked steps are taken at the current length without a check, the documented intent of DREAM's checkevery. DREAM's code takes them from a stale initTime, with the lengths h₁ − h₀, −h₀, −(h₀ + h₁), and checkevery + 1 of them (TimeStepperAdaptive.cpp:119-122, 145-152, 334-338), so no option reproduces it.

The tolerances are the Tolerances the system holds for the time-step error (time_step_tolerances of EquationSystem); an unknown without one gets relative_tolerance and no absolute tolerance (DREAM's default absolute tolerance is zero for every unknown but the runaway density, the runaway distribution f_re and the runaway current, which get 1e-10, 1e-10 and e c 1e-10; ConvergenceChecker::GetDefaultAbsTol, ConvergenceChecker.cpp:110-120). An unknown with both tolerances zero is left out of the check. A step whose solve fails is retried with failure_factor times its length, at most max_failures times in a row (DREAM's STEP_REDUCTION_AT_EXCEPTION and MAX_STEPS_WITH_EXCEPTION); in a step between checks DREAM ends the run instead, Runes2P retries there too. DREAM's step update reads the tolerances of the wrong unknowns when a prescribed unknown precedes others; Runes2P uses one set of tolerances for the acceptance and the next step.

The first step has the length initial_step (nothing: the first interval of the time points). Every step ends at or before the next time point of the run, and a save time must be a time point, so that every saved state is a computed one, as with FixedSteps (DREAM saves computed steps only and has no save times). A step shortened by a time point sets the next length from its own error, as a step to tmax does in DREAM.

TimeSolution

Result of integrate: the layout, the saved times, the saved states (one column per saved time; views of the unknowns through unknown_view), the SolverStatistics, the record of the run's settings (run_settings) and update_states, per interval update the state in force when each saved state was computed (one column per saved time). Of a resumed run it holds the states saved by this segment, those after the time it resumed from.

TimeStepFailure

Thrown by integrate when a step still fails after the retries the stepper allows; holds the time and time_step of the last attempt, its Newton status and the statistics up to the failure.

TimeStepper

Supertype of the steppers of integrate, which set the length of every step: FixedSteps, StepDoubling and StepController. The state a stepper carries between steps (its step length, its counts) is recorded in every checkpoint and restored by a resume.

integrate(system, u0, times; solver=NewtonSolver(), stepper=FixedSteps(), save_times=times,
          description=Dict(), updates=(), checkpoint=nothing, resume=false, stop=nothing,
          log=nothing)
integrate(system, u0, times, options::RunOptions)

Backward-Euler integration of system from the state u0 at times[1] through the time points times (any increasing array, e.g. nonuniform or matching the restarts of another run). The stepper (TimeStepper) sets the length of every step, and every time point is the end of a step: with FixedSteps, the default, each interval between two time points is one step, retried with half the time step when it fails (max_step_halvings = n is a shorthand for stepper = FixedSteps(max_halvings=n)). Each step is solved by newton_solve! with explicit terms at the previous time level. A step that still fails throws TimeStepFailure.

States are saved at save_times (strictly increasing, inside [times[1], times[end]]; a save time within a relative 1e-12 of a time point is taken as that time point, and two save times taken as one time point are refused), and every saved state is a computed one: with FixedSteps and StepDoubling a save time must be a time point (an ArgumentError names one that is not; DREAM saves computed steps only), while StepController ends a step on every save time. time_points builds times and save_times from Legs in one arithmetic. Returns a TimeSolution.

The solution carries the record of the run's settings (run_settings), taken before the first step and before a stored state is restored, so it describes the setup as the script built it and is the same in every segment of a run. It holds the grids, the layout, every term, the solver, the stepper, the time points, the code version with the commit of Runes2P, and description, a dictionary of plain values the caller adds (a run name, the origin of the initial state).

updates are IntervalUpdates, applied at the start of the run and after every accepted step, before the step that follows (the knock-on cutoff that follows the plasma, KnockOnCutoffUpdate). An update that reports a changed setup discards the stored Jacobian.

checkpoint (CheckpointSettings) writes checkpoints after accepted steps and the output of the segment with each of them; it needs HDF5 loaded. A checkpoint holds the state of the stepper, so the resumed run repeats the uninterrupted one.

resume continues a run from a stored state: true takes the last checkpoint of checkpoint when there is one and starts from u0 otherwise, so the same script serves every segment of a run; a file name takes the checkpoint or the last saved state of an output file; a StoredState is what read_checkpoint and read_state return. The stored time must lie inside the time points of the run; a state from an output file, or from a checkpoint without the record of the stepper, must be at one of them, equal to the last bit, so that the stepper can start there. The run continues from the stored time; a script that builds its time points the same way in every segment gets the same numbers. The unknowns are matched by name: those in the stored state replace the values of u0, on the same grids, and an unknown the stored state does not have keeps its value of u0, or takes it from its equation (below). The updates are set back to their stored states. The settings may differ from those of the run that stored the state (tolerances, the solver, later time points, terms); the record names the sections that changed (segment in the settings) and a warning lists them.

An unknown whose equation has initial_value = :from_equation (the energy of the cold electrons, thermal_energy_equation) is set from its row at the start of the run, whatever u0 holds there (fill_initial_values!); a resumed run sets those that the stored state does not have. The record of the segment names them (initial_values_from_equations).

stop (StopCondition) ends the segment after an accepted step before the last time point, when the wall time left does not cover another step or a stop file appears. The segment writes its checkpoint and returns what it has saved; run_status of the solution is then :stopped, and exit_status is the exit status for the script.

log (ProgressLog) reports while the run is going: a header, entries at milestones of the simulated time, after a silence and on demand, and the events of the run when they happen. The warnings of the run are written to it.

read_settings(filename)

The settings group of an output file as a plain record, equal to the record written (run_settings); empty for a file without the group. Defined by the package extension Runes2PHDF5Ext, which loads with HDF5 (using HDF5).

run_settings(system; solver=NewtonSolver(), stepper=FixedSteps(), times=Float64[],
             save_times=times, description=Dict(), updates=(), checkpoint=nothing,
             stop=nothing, log=nothing)
run_settings(system, times, options::RunOptions)

The record of a run, written by write_output as the settings group. Sections:

• grids: the kinetic grids by name, with their faces;
• layout: the unknowns, their order, sizes, scales and lower limits;
• system: the equations with every term and its options, the tolerances;
• updates: the interval updates by name (IntervalUpdate);
• prescribed: the tables of the prescribed functions of time, sampled at the time points (describe_prescribed_plasma), which the terms refer to by name. They are a section of their own so that the hash of system does not change with the time points;
• solver: the Newton and linear solver settings;
• time: the time points, the save times and the stepper (TimeStepper; max_step_halvings = n stands for stepper = FixedSteps(max_halvings=n));
• run: description, what the caller adds (a name, the origin of the initial state, the parameters a hand-built coefficient came from), as plain values;
• code: code_version, with the commit of Runes2P;
• environment: run_environment;
• checkpoint: the checkpoint settings, or the type none;
• stop: the stop condition, or the type none;
• log: the settings of the progress log, or the type none;
• hashes: settings_hash of each of the first eight sections, and setup_hash of the eight together. code, environment, checkpoint, stop and log are outside the hashes, so the same setup has the same hash on another machine and commit;
• complete: false when some object of the run has no record of its own (a term defined outside the package, a prescribed function that fails at a time point).

integrate builds the record at its start and returns it in the solution, with the record of the segment added under segment: its number, its start time, the number of states saved before it, its status (running, stopped or finished) with the time reached, the intervals done and the reason of a stop, and for a resumed segment the file it resumed from, the commit and setup hash of the run that wrote it and the sections whose hash changed.

write_output(filename, solution::TimeSolution; atomic=false)

Write a TimeSolution to an HDF5 file: the root attributes format_version (OUTPUT_FORMAT_VERSION), runes2p_commit, runes2p_version and runes2p_tree_clean, the record of the run under settings (run_settings, read back by read_settings), the describe record of the layout under layout, the saved time, one dataset per unknown under unknowns (scalar (t), radial (components, r, t), kinetic on a shared grid (p, xi, r, t), on per-radius grids a group with one (p, xi, t) dataset per radius named radius_<ir>), and the SolverStatistics under solver_statistics with their own time array, and the states of the interval updates per saved time under update_states. The string attribute unit holds the unit of every entry of unknowns (unknown_unit, empty for an unknown without one), of time ("s") and of the faces of the grids under layout ("m", "m_e c" and "1"). atomic = true writes to a temporary file that is renamed onto filename (write_atomically). read_output reads the file back. Defined by the package extension Runes2PHDF5Ext, which loads with HDF5 (using HDF5).

Sections of the settings that define the setup of a run, covered by setup_hash.

The error estimate per unknown of the last accepted step, for the log; nothing without one.

Record of the state of a stepper, stored in the checkpoint.

stepper_state(stepper, system, times, t, record)

The state of stepper at the time t of the run of system with the time points times: restored from the record of a checkpoint written by the same kind of stepper, fresh otherwise.

Leg(duration; steps, saves)
Leg(duration; dt, save_dt)
Leg(duration; save_dt)

One leg of a run, duration seconds long, in one of three keyword forms:

• steps time steps of equal length with saves saved states, at the save steps of DREAM's nSaveSteps: a running sum of steps / saves in floating point, rounded half away from zero, so that a tie the sum reaches exactly goes to the later step and one it falls just short of to the earlier (Leg(T; steps = 5, saves = 2) saves after the steps 3 and 5);
• steps of length dt [s] with a saved state every save_dt [s], where dt and save_dt divide the duration (to 1e-12 relative) and the saved states are those of the first form with duration / save_dt saves;
• for a StepController alone, no prescribed steps and a saved state every save_dt [s], which divides the duration.

The last time point of a leg is its end, and it is saved. Any other combination of keywords, a positional form or a non-positive value throws an ArgumentError that names the leg. time_points turns a list of legs into the times and save_times of integrate.

time_points(legs; start=0.0)

The times and save_times [s] of a run made of legs (Leg), one after the other from start [s], as a pair of vectors. Both begin with start. A prescribed-step leg adds its time points, with the saved ones among them; a controller leg adds only its end to the time points and its save times in between, which a StepController makes ceilings of its steps. The end of each leg is start plus the durations so far, in one arithmetic for both vectors.