Unknowns and the equation system

A run of Runes2P solves one system of equations for one state vector $u$. The UnknownLayout places the unknowns in $u$. Each unknown has one Equation, whose right-hand side is a sum of EquationTerms. The EquationSystem collects the equations into

\[M \dot u = \mathrm{rhs}(u, t),\]

with a constant mass matrix $M$, checks that the parts fit together, and supplies the residual and the Jacobians to the solvers.

Unknowns

An unknown has a name, a kind, a reference scale, a unit and optionally a LowerLimit.

The momentum grid of a kinetic unknown may differ between radii (see Grids). The components of a radial unknown split into consecutive groups (component_groups). For the ion densities of ion_densities_unknown a group holds the charge states of one species. Lower limits and tolerances can be relative to the sum over a group, so a trace impurity is converged on its own density.

Layout of the state vector

Each radius owns one contiguous block of $u$. Within it come the radial unknowns, then the kinetic unknowns, each in the order of declaration. A kinetic unknown is stored with $p$ fastest, then $\xi$. The scalar unknowns form one block after the radius blocks. UnknownLayout(...; scalars = :first) places that block first, which serves only to check that no code depends on the order.

Code reaches the entries only through the layout, never through offsets of its own:

• unknown_id turns a name into the integer id used in inner loops.
• radial_view, kinetic_view and scalar_view give type-stable views, also of vectors of dual numbers.
• kinetic_array and unknown_values present an unknown as a 3D array or as its values over the radii to the operators.
• unknown_range, radius_block_range and scalar_block_range give index ranges, which the block preconditioners use.

The layout is immutable. copy_unknowns! carries a state by name to a layout with another order or other unknowns, as a resumed run does. unknown_scales gives the reference scale of every entry, which the tolerances and the row scales build on.

Units of the unknowns

An unknown carries its unit as a string, which unknown_unit reads and the output file stores (Unknowns). The constructors and equations of the package give their unknowns the default names below, and a scalar or radial unknown of such a name takes the unit of this table. The column "In the file" gives the unit in the ASCII form of the file. The keyword unit of ScalarUnknown, RadialUnknown and KineticUnknown sets another unit, an ASCII string without control characters. Without that keyword, a KineticUnknown takes "m^-3" whatever its name, and a scalar or radial unknown of a name outside the table has the empty string as its unit.

A keyword sets another name, such as total_energy of total_energy_equation for :total_electron_energy. induction_unknowns and ion_densities_unknown keep the unit of the default name under another name. An unknown declared by the caller under another name takes its unit from the keyword unit.

A distribution and the density of the hot electrons have no default name. The caller names them, as the examples name the distribution :f. The density of the hot electrons is in m⁻³, which the caller's declaration passes as unit = "m^-3". A distribution is in m⁻³ per unit of momentum volume (Notes on the units), hence the default unit of a KineticUnknown.

Equations

An Equation names its unknown $x$, its kind and its terms.

Pairing of rows and unknowns

TimeDerivative(of = y) puts $\dot y$ in the row of $x$, so the row of $y$ can be algebraic and each row still holds one unknown. The induction equations of induction_equations follow the pairing of DREAM. The row of the parallel electric field $E$ holds Faraday's law, $\dot \psi = 2\pi E$, with $\psi$ the poloidal flux divided by the major radius (V s/m), as in DREAM. The row of the flux holds Ampère's law, which is algebraic in $\psi$. The plasma current, the flux at the edge and the flux at the wall are the scalar unknowns of induction_unknowns. The equation of the total current density, the sum of the currents of the run, comes from the caller.

The cold-electron temperature is paired in the same way. Its row holds the time derivative of the thermal energy $W_\mathrm{cold}$, whose own row is the algebraic $0 = \tfrac32 e\, n_\mathrm{cold} T_\mathrm{cold} - W_\mathrm{cold}$, with $T_\mathrm{cold}$ in eV (thermal_energy_equation). In the fully kinetic mode the row holds the derivative of the total energy instead (total_energy_equation, see Energy and current balance).

A Constraint fixes a rate that appears in other equations. The magnitude of a particle source is of this kind: its row is the density that the distribution must keep, and the source stands in the kinetic equation. The source enters only through $\dot f$, and the constraint does not contain it, so the system has differential index 2 in that unknown, which backward Euler solves.

Scales and initial values

The scale of a row is the magnitude of the quantity it balances. By default it is the scale of the unknown whose time derivative the row holds, or of the row's own unknown. Equation(...; scale) sets it where the row balances another quantity: Ampère's law, and every Constraint. residual_scales divides differential rows by the time_scale of the system.

Equation(...; initial_value = :from_equation) sets an algebraic unknown from its row at the start of a run, as DREAM's INITRULE_EVAL_EQUATION does (fill_initial_values!). The fill solves each entry from its own row, so it refuses a row that couples the entries of its own unknown (radial transport, Ampère's law), unless the values solved entry by entry satisfy it. Such an unknown is set in the initial state. The thermal energy of the cold electrons takes its initial value this way. fill_maxwellian! fills a distribution with a Maxwellian from profiles of the density and the temperature, the start of a fully kinetic run in DREAM. fill_runaway_distribution! fills a distribution from a density in one of DREAM's shapes of f_re (RUNAWAY_SHAPES). fill_coronal_equilibrium! puts species of the ion densities in coronal equilibrium at the free-electron density they give, DREAM's init_equil. It leaves the electron densities to the caller, who sets them from the rows of the system or by hand.

Terms

A term adds its contribution to the rows of one unknown, the row unknown, and may read any unknowns (dependencies). rhs! evaluates it as a pure function of the state and the time. A few properties decide how a term enters the system and its Jacobian:

A row that holds every cell of a radius, such as a density moment, would take one color per cell in the colored AD. Such terms assemble their block directly, and the coloring sees the other terms alone. The flux through the momentum interface (InterfaceParticleFluxTerm, InterfaceEnergyFluxTerm) reads the cells at the face of both grids and the field in one row, and assembles that row directly for the same reason. A term whose rows read all charge states of a radius builds its block by one dual-number pass over a reduced layout (Runes2P.add_local_jacobian!, Runes2P.local_layout).

The knock-on term with an implicit primary (KnockOnTerm) is matrix-free, since its block is dense per radius. With a lagged primary the knock-on term assembles a sparse block directly.

A term defined outside the package adds methods to these functions under their qualified names, Runes2P.rhs!(out, term::MyTerm, ...).

The system

EquationSystem sorts the terms into implicit sparse, implicit matrix-free and explicit terms. mass_matrix has an entry 1 in each row with a time derivative, in the column of the unknown whose derivative the row holds. It is diagonal when every such row differentiates its own unknown.

Residual and Jacobian

A backward Euler step solves

\[F(u) = M \frac{u - u_\mathrm{ref}}{h} - \mathrm{rhs}(u, t) = 0\]

with residual!. The explicit terms are evaluated at the start state, and prepare_step! sees the start state before every evaluation. For backward Euler $h = \Delta t$ and $u_\mathrm{ref}$ is the start state. For the variable-step BDF2 of StepController $h = c\,\Delta t$ and $u_\mathrm{ref}$ is a combination of the two last states, with $c$ set by the ratio of the last two steps. Newton's method (NewtonSolver) uses the Jacobian $M/h - \partial \mathrm{rhs}/\partial u$ (see Time stepping and solvers).

The stored part of $\partial \mathrm{rhs}/\partial u$ (jacobian!) has three sources:

• The terms with automatic assembly add their entries by colored forward-mode AD (jacobian_colors counts the passes).
• The terms with direct assembly add their own entries.
• The terms whose coefficients follow unknowns add their plasma columns (add_plasma_columns!).

jacobian_operator adds the matrix-free terms and the corrections of approximate blocks, so the operator is the exact Jacobian. The structure (jacobian_sparsity, jacobian_prototype) is fixed when the system is built. The construction evaluates the sparse terms at initial_state and initial_time. A term that is defined only within the run's interval or at a physical state therefore builds without error.

rhs!(du, u, system, t) is the ODE form for a general integrator. There a lagged part sees the state itself, so ode_function passes no Jacobian for a system with a lagged part (Runes2P.lags_in_ode_form).

Tolerances

A Tolerance passes an entry when

\[|\delta| \le a + r|x| + g \sum_{k \in \text{group}} |x_k|,\]

with $g = 0$ without a group part. tolerance gives the default values. A radial unknown with several groups, such as the ion densities, takes a group part for the Newton iteration, so a trace species is converged on its own density. The system holds separate tolerances for Newton and for the error control of the time steps.

Construction checks

The constructor refuses a system in these cases:

• An unknown has no equation, or more than one.
• A term reads a name that is not in the layout, or fails its own check against the layout.
• An algebraic row has no implicit term that reads its unknown.
• A constraint has no scale, or no other equation reads its unknown implicitly.
• A TimeDerivative(of) pairs unknowns of different kind or sizes, or nothing determines the row's unknown.
• The time derivative of one unknown appears in two rows.
• A tolerance with a group part belongs to an unknown other than a radial one.

The algebraic block is not checked for numerical regularity.

Coefficients that follow the unknowns

Terms read the plasma through one interface, a PlasmaSource evaluated by plasma_values. A PrescribedPlasma is prescribed in time. A PlasmaFromUnknowns takes the cold-electron temperature and optionally the density and the ion charge states from radial unknowns, and the rest from a prescribed plasma. A source that reads the ions reads the density as well. is_state_dependent tells the two apart.

Caches and rebuilds

The values of a PlasmaFromUnknowns carry no derivative parts, so the coefficients of a term hold plain numbers. A cache of coefficients keeps the plasma values per radius it was built from in a PlasmaKey. changed_plasma! compares them with update_key! and returns the plasma only when the values of some radius differ. A FokkerPlanckTerm then rebuilds its coefficients at all radii (refresh_coefficients!). A prescribed source is compared once per time, so the evaluations of a term at one time share one comparison.

Plasma columns of the Jacobian

The derivatives with respect to the plasma unknowns enter the stored Jacobian as columns of their own (has_plasma_columns). dual_plasma seeds the temperature as partial 1, the density as partial 2 and the charge states after them. The coefficients at a radius depend only on the plasma of that radius, so one evaluation gives the derivatives at every radius. Runes2P.plasma_columns maps the partials to entries of the state, and plasma_column_sparsity! lists them in the pattern. add_plasma_columns! applies the dual coefficients to $f$ with the scheme's own face weights, so the plasma columns are exact also when the stored block in $f$ is approximate. The dual coefficients are cached under a PlasmaKey of their own and rebuilt only when the plasma values of some radius change. The term and the terms that read its coefficients (the boundary flux, the interface, the flux through it) share that one rebuild per state.

Lagged parts

The keyword lagged of FokkerPlanckTerm names the part of the plasma taken at the start of the step:

With lagged = :all the plasma is fixed for the step, so the block in $f$ is the term's whole derivative, exact unless jacobian selects an approximation. The default :ions is a choice of this package (see Choices of this package). A term that lags the ions is of first order under BDF2. The kinetic ionization and a corrected conductivity read the charge states of the iterate, so within a step the ions enter at two states, a difference of first order in the step.

ConductivityModel(...; kinetic = :discrete_operator) takes a prescribed plasma. Its correction comes from the steady state of the discrete operator, and following the unknowns would cost a factorization of a bordered operator per radius at every rebuild.

Temperature floor

A PlasmaFromUnknowns raises a DomainError for a temperature at or below its floor. The Newton iteration reports it as a failed solve, and the stepper answers with a shorter step. The floor is also the lower limit of the temperature unknown in the Newton iteration (temperature_floors). A temperature that only fluid terms read takes the default floor DEFAULT_TEMPERATURE_FLOOR unless it declares a limit of its own (fluid_temperature_floors). The value of the floor is a choice of this package (see Choices of this package). How the Newton iteration keeps an unknown above its limit is on the page Positivity and flux limiters.

Pairing checks

Terms that describe one physical process from two sides must agree on their plasma, temperature, treatment, region, ion unknown, composition or tables. A mismatch gives a result that is wrong without an error, so EquationSystem refuses it at construction.

Plasma and treatment

Shared coefficients

Hot region and kinetic ionization

Ions

Rows of the total energy

A row that holds the time derivative of $W_\mathrm{tot} = W_\mathrm{cold} + W_\mathrm{hot}$ charges every unlisted change of $W_\mathrm{hot}$ to the cold electrons. For each distribution in $W_\mathrm{hot}$ the system therefore requires in that row:

• the sinks of its collisions where they have any (KineticEnergyLossTerm with part = :sinks), over the region of the energy moment.
• the energy that leaves through $p_\mathrm{max}$, as MomentumBoundaryFluxTerm with weight = :energy or as InterfaceEnergyFluxTerm of the interface to the runaway grid.

Each of these terms takes a negative coefficient.

The system refuses radial transport of that distribution, RadialEdgeLossTerm, KnockOnTerm and RosenbluthPutvinskiSourceTerm on it, since no term of the row carries their energy. The energy balance itself is on the page Energy and current balance.