Equations and the system

Algebraic()

Kind of Equation 0 = Σ terms that holds its unknown. See also TimeDerivative and Constraint.

Constraint()

Kind of Equation 0 = Σ terms that does not hold its unknown: the unknown is fixed through the equations of other unknowns, which hold it. The magnitude of the particle source is of this kind: its equation is the density that the distribution has to keep, 0 = n(f) − n_cold − n_hot, and the magnitude stands in the kinetic equation (DREAM's S_particle). The system then has the index 2 in that unknown, which backward Euler solves, and the diagonal entry of the row in the Jacobian is zero. A constraint names the scale of the quantity it balances (Equation(...; scale)): its unknown is a rate, whose scale does not fit the row.

Equation(unknown, kind, terms...; scale=nothing)

Equation of the unknown unknown (a name in the layout): TimeDerivative(), TimeDerivative(of) or Algebraic(), with the right-hand side the sum of terms (EquationTerms).

scale is the magnitude of the quantity whose rate or balance the row holds, for the weights of the residual (residual_scales); on the differential rows of a kinetic unknown, Newton takes the larger of it and the entry at the start of the step (NewtonSolver). By default it is the scale of the unknown whose time derivative the row holds, and the scale of the unknown of the equation for an algebraic row. An algebraic row that balances another quantity than its unknown names the scale of that quantity: the row of the poloidal flux holds Ampère's law.

initial_value = :from_equation makes integrate set the unknown of an algebraic equation from its row at the start of a run (DREAM's INITRULE_EVAL_EQUATION), in place of the value in the initial state: fill_initial_values!. The default :given takes the initial state as it is.

TimeDerivative()
TimeDerivative(of)

Kind of Equation du/dt = Σ terms. The other kinds are Algebraic and Constraint. TimeDerivative(of) puts the time derivative of another unknown into the row, d(of)/dt = Σ terms, for an unknown of of the same kind and sizes as the unknown of the equation. Faraday's law is of this form: the equation of the electric field is dψ/dt = 2π E, with the poloidal flux as of, and the equation of the flux is Ampère's law, which is algebraic. Every algebraic row then holds its own unknown. The unknown of an equation with the derivative of another unknown must appear in the terms of that equation, or the equation of of must reach it through algebraic rows: the temperature in DREAM's pairing of the energy of the cold electrons.

RUNAWAY_SHAPES

The shapes of fill_runaway_distribution!: :forward, :xi_positive, :xi_negative and :isotropic, DREAM's INIT_FORWARD, INIT_XI_POSITIVE, INIT_XI_NEGATIVE and INIT_ISOTROPIC of f_re (src/Settings/Equations/f_re.cpp, DREAM v26.5).

EquationSystem(layout, equations; time_scale=1, newton_tolerances, time_step_tolerances,
               initial_time=0, initial_state=nothing)

The system M du/dt = rhs(u, t) for the unknowns of layout, with one Equation per unknown. On construction it checks that every unknown has exactly one equation, that every term's dependencies are unknowns of the layout (and passes the term's own checks), and that every algebraic unknown appears in its own equation (some implicit term depends on it; the numerical regularity of the algebraic block is not checked). The terms are split into three groups:

• implicit :sparse terms, whose Jacobian is stored (jacobian!): by colored AD over the terms with automatic assembly, and by their own entries for the terms with jacobian_assembly(term) === :direct (direct_terms), whose rows the coloring does not see;
• implicit :matrix_free terms, which enter the Jacobian only through jacobian_operator;
• explicit terms, evaluated at the previous time level in residual! and absent from the Jacobian.

Sparse terms whose stored block is an approximation (jacobian_kind(term) === :approximate) are also listed in correction_terms; stored_jacobian_exact is false when there are any, and linear is true when every implicit term is linear in the state (is_linear), both read by the Newton iteration.

time_scale [s] sets the scaling of differential rows (residual_scales); newton_tolerances and time_step_tolerances map unknown names to Tolerances (defaults from tolerance). The time-step tolerances are those of the error control of the adaptive steppers (StepDoubling, StepController).

The construction evaluates the implicit :sparse terms once, at the state initial_state (the zero state by default) and the time initial_time [s], for the structure of the stored Jacobian. A KineticEnergyLossTerm of a source with lagged = :all that was never prepared is prepared there at initial_state (prepare_step!) before its block is assembled; one prepared before keeps its plasma. Set initial_time to the start of the run when a term is defined only on the run's time interval (the plasma_at of a prescribed_fokker_planck_term), and initial_state to the initial state when a term needs a valid state, such as a Fokker-Planck term whose plasma reads the temperature from an unknown (PlasmaFromUnknowns).

Tolerance(absolute, relative; group=0)

Convergence tolerance of one unknown: an entry passes when |δ| ≤ absolute + relative |x|, with absolute in the unknown's physical units. group adds group Σ |x_k| over the components k of the entry's group at its radius (component_groups): the tolerance of a charge state relative to the density of its species, so that a trace impurity is converged on its own scale (DREAM tests the 2-norm of the whole ion block with an absolute tolerance of 0, src/ConvergenceChecker.cpp:110-120, 193). A group part needs a radial unknown: EquationSystem refuses it on a kinetic or scalar unknown.

fill_initial_values!(u, system, t; unknowns=nothing)

Set the unknowns whose Equation has initial_value = :from_equation from their algebraic rows at the time t [s], entry by entry (DREAM's INITRULE_EVAL_EQUATION; the energy of the cold electrons from their temperature and density). A row 0 = c x + r(u) that is linear in its own unknown x gives x = −r/c. The coefficient c comes from the row at x = 0, s, 2s, with s the scale of the unknown, so a row needs no particular form. A row that reads another filled unknown is solved again until every filled row holds, whatever the order of the equations. unknowns limits the fill to those names, the unknowns that a resumed run adds to the stored layout. Returns the names of the filled unknowns.

The coefficient c of an entry is that of its own row with every entry of the unknown moved at once, so a row that reads other entries of its own unknown (a radial transport term, Ampère's law for the poloidal flux) has no value here: its c would be a sum over the entries it reads. Such a row is refused unless the values solved entry by entry satisfy it.

Throws an ArgumentError for a row that does not depend on its unknown or is not linear in it, for a row that reads other entries of its unknown where its coefficient vanishes or it does not hold after the fill, and for filled rows that read each other in a cycle.

fill_maxwellian!(u, layout, unknown, density, temperature; shape=:relativistic_maxwellian)

Fill the kinetic unknown unknown of the state u with a Maxwellian of the density density [m⁻³] and the temperature temperature [eV], profiles in any form of radial_values on the cell centres of the radial grid. Returns u.

The shape is that of ParticleSourceTerm at the cell centres: exp(-(γ − 1)/Θ) with Θ = T/mₑc² for :relativistic_maxwellian, or exp(-p²/2Θ) for :maxwellian. It is scaled at each radius so that its moment density on the grid equals the profile. This is DREAM's initial distribution of a fully kinetic run: the relativistic Maxwellian at the cell centres, rescaled to the density on the grid (ConstructEquation_f_maxwellian in src/Settings/Equations/distribution.cpp, with rescaleMaxwellian true for f_hot in src/Settings/Equations/f_hot.cpp, DREAM v26.5).

Throws an ArgumentError for an unknown that is not kinetic, another shape, a density that is negative or not finite, or a temperature that is not positive or not finite, and a DimensionMismatch for a profile with another number of values than radii; these leave u as it was. A DomainError for a Maxwellian with no density on the momentum grid (a temperature far below the first momentum cell) comes after the radii before it are filled.

fill_runaway_distribution!(u, layout, unknown, density; shape=:forward, electric_field=nothing)

Fill the kinetic unknown unknown of the state u with the electrons of the density density [m⁻³], a profile in any form of radial_values on the cell centres, all in the first momentum cell of its grid, so that their moment density on the grid equals the profile. These are DREAM's starts of f_re from n_re (src/Settings/Equations/f_re.cpp:147-208, DREAM v26.5). A shape of one pitch cell puts the value $n V' / (\mathcal{V}'_{1j} \Delta p_1 \Delta\xi_j)$ in the cell $(1, j)$, with the spatial Jacobian $V'$ of the radius and the phase-space Jacobian $\mathcal{V}'$ of the cell; the isotropic shape puts $n V' / (\Delta\xi_\mathrm{grid} \mathcal{V}'_{1j} \Delta p_1)$ in every pitch cell $j$, with the pitch range $\Delta\xi_\mathrm{grid}$ of the grid, 2 on DREAM's grids from −1 to 1:

• shape = :xi_positive: the last pitch cell, at the largest ξ of the grid (INIT_XI_POSITIVE);
• :xi_negative: the first pitch cell, at the smallest ξ (INIT_XI_NEGATIVE);
• :forward: the last pitch cell where the field electric_field [V/m] (a profile as density, of which only the sign counts) is not negative, the first where it is (INIT_FORWARD, DREAM's default);
• :isotropic: every pitch cell, with f 𝒱' uniform in ξ (INIT_ISOTROPIC).

DREAM's analytic avalanche distribution (INIT_AVALANCHE) is not one of the shapes. The other cells of unknown are set to zero. Returns u.

Throws an ArgumentError for an unknown that is not kinetic, another shape, :forward without a field, a density that is negative or not finite, or a field that is not finite, and a DimensionMismatch for a profile with another number of values than radii; these leave u as it was.

jacobian!(J, system::EquationSystem, u, t)

Values of the stored Jacobian ∂rhs/∂u of the implicit :sparse terms at (u, t), with the time t [s], into J (from jacobian_prototype), by colored forward-mode AD and by the entries of the terms that assemble their own block (jacobian_assembly); for terms with an approximate Jacobian (jacobian_kind) their stored block, not the derivative. Matrix-free and explicit terms are not included.

jacobian_colors(system)

Number of colors of the stored Jacobian's column coloring: the number of directional derivatives (dual-number passes over the sparse terms) one jacobian! costs. The terms that assemble their own block do not count.

jacobian_operator(J, system::EquationSystem, u, t)

∂rhs/∂u at (u, t), with the time t [s], as a linear operator: the stored matrix J (from jacobian! at the same point) plus the products of the implicit matrix-free terms and the corrections of the terms with an approximate stored block (jacobian_correction_product!), so the operator is the exact Jacobian whatever the stored blocks hold. u is referenced, not copied. The products go through a buffer owned by the system, shared by all operators of the same system, so they must not be applied concurrently. The Newton-Krylov solver GMRES builds its own operator.

jacobian_prototype(system)

A new sparse matrix with the structure of jacobian_sparsity, to be filled by jacobian!.

jacobian_sparsity(system)

Structural pattern of the stored Jacobian: the diagonal and the implicit :sparse terms.

mass_matrix(system)

Mass matrix M, constant: 1 in the rows with a time derivative, in the column of the unknown whose derivative the row holds, and nothing in the algebraic rows. A Diagonal when every row holds the derivative of its own unknown, a sparse matrix otherwise (TimeDerivative(of)).

ode_function(system::EquationSystem)

SciML ODEFunction of the system, for OrdinaryDiffEq and NonlinearSolve. Defined by the package extension Runes2PSciMLBaseExt, which loads with SciMLBase (using SciMLBase).

residual!(F, u, u_old, system::EquationSystem, t, dt; reference=u_old)

Backward-Euler residual F = M (u - reference)/dt - rhs(u, t) at the time t [s] that ends the step of length dt [s], with explicit terms at u_old and the terms prepared from u_old (prepare_step!). The reference is the start state for backward Euler; the two-step BDF2 of StepController(order = 2) passes its combination of the two last states and its effective step. M has one entry per row with a time derivative, in the column of the unknown whose derivative the row holds. Allocation free.

residual_scales(system)

Scale of every row of the residual: scale/time_scale for rows with a time derivative and scale for algebraic rows. scale is the scale of the equation when it names one, else the reference scale of the unknown whose time derivative the row holds, else that of the row's unknown. Dividing the residual by these makes all rows of order one. The Newton iteration of NewtonSolver does not use them: it weights differential rows by dt/scale (on a kinetic unknown with the larger of the scale and the entry at the start of the step), since round-off in (u - u_old)/dt at a short step would exceed a tolerance scaled by a fixed time_scale.

rhs!(du, u, system::EquationSystem, t)

Right-hand side of M du/dt = rhs(u, t) with all terms at u and the time t [s] (the ODE form, as passed to a general integrator).

tolerance(system, name; kind=:newton)

Tolerance of the unknown name for Newton convergence (kind = :newton) or time-step error control (kind = :time_step); by default relative = 1e-10 and absolute = 1e-10 times the unknown's reference scale. A radial unknown that declares groups of its components (component_groups, the ion densities, also of a single species) takes group = 1e-10 and absolute = 1e-20 times its scale for the Newton iteration: each species is converged relative to its own density, and the absolute part only covers a species absent at a radius.

fluid_temperature_floors(system)

The temperature unknowns that only the fluid terms read (no PlasmaFromUnknowns gives them a floor), as pairs name => floor with the default floor of PlasmaFromUnknowns, 0.5 eV (DEFAULT_TEMPERATURE_FLOOR). The Newton iteration takes each as the lower limit of an unknown that declares none of its own: at and below zero the ADAS coefficients vanish, so the radiation would switch off, and the Rechester-Rosenbluth transport has no thermal speed. An unknown with a LowerLimit of its own keeps it, so a comparison that follows DREAM below 0.5 eV can set a lower one. SecondaryHeatingTerm is in the list as well: its temperature is that of the cold electrons it heats, though it enters no coefficient of the term.

temperature_floors(system)

The floor [eV] of every temperature unknown that a PlasmaFromUnknowns of the system reads, as pairs name => floor, with the largest floor where several sources read one unknown. The Newton iteration takes each as the lower limit of that unknown.