Equation terms

EquationTerm

Supertype of the terms of an Equation. A term implements

• dependencies(term): names of the unknowns it reads;
• coupling(term): :sparse (its Jacobian block is stored; sparsity! lists the structural nonzeros) or :matrix_free (applied through jacobian_vector_product!, e.g. the knock-on term, whose block is dense per radius);
• rhs!(out, term, u, layout, row, t): add the term to out in the rows of the unknown row, as a pure function of the state u and the time t. :sparse terms must accept state vectors of dual numbers, since their Jacobian comes from AD;
• sparsity!(entries, term, layout, row) for :sparse terms: push (row index, column index) pairs of the structural nonzeros;
• jacobian_vector_product!(out, term, u, v, layout, row, t) for :matrix_free terms: add (∂ rhs_term / ∂u) v to out;

A :sparse term whose row holds many unknowns (a moment of a distribution, the loss through the edge) declares jacobian_assembly(term) === :direct and adds its Jacobian entries itself with add_jacobian!, so that the coloring of the AD Jacobian does not see its rows.

A term may also define treatment(term) (:implicit by default; :explicit terms are evaluated at the previous time level and contribute nothing to the Jacobian), prepare_step! (a part of an implicit term lagged to the previous time level), validate(term, layout, row) (checks against the layout, run when the system is built), is_linear(term) (false by default: whether the term is linear in the state, up to a constant, which lets a linear system accept a Newton step on the update test alone) and, for :sparse terms whose stored block is an approximation of their Jacobian (a quasi-Newton matrix), jacobian_kind(term) === :approximate with jacobian_correction_product! supplying the difference to the exact product.

A term defined outside the package adds its methods under the qualified names (Runes2P.dependencies(term::MyTerm) = ..., function Runes2P.rhs!(out, term::MyTerm, ...)). The functions that callers use are exported. validate, whose name is common to many packages, add_local_jacobian!, lags_in_ode_form and the plasma columns of a term (has_plasma_columns, plasma_column_sparsity!, add_plasma_columns!), which a term implements and the system calls, are declared public and not exported.

LayoutKineticArray

3D array view f[i, j, ir] of a kinetic unknown with the same cell counts at all radii in a state vector, so that operators written for 3D arrays apply to the state directly.

LayoutRadialValues

The values of a radial unknown with one component over the radii, as a vector view of a state vector: values[ir] is the value at radius index ir. Operators take the electric field of the state in this form. unknown_values builds it.

LowRankBlock(rows, columns, left, right)

Part left * right' of a term's Jacobian, with rows rows and columns columns of the state vector (contiguous ranges). left is length(rows) × k, right is length(columns) × k, k the rank.

add_jacobian!(J, term, u, layout, row, t)

Add the entries ∂ rhs_term / ∂u of a term with jacobian_assembly(term) === :direct to the stored matrix J at the state u and the time t. The entries lie in the pattern of sparsity!. A method may allocate work arrays per call, such as the dual numbers of ParticleSourceTerm, InterfaceParticleFluxTerm and InterfaceEnergyFluxTerm or the row of KineticEnergyLossTerm.

add_local_jacobian!(J, term, u, layout, row, t)

Add the Jacobian of a term whose rows at a radius read only the radial unknowns of that radius (its dependencies) to J, by one evaluation of its rhs! with dual numbers: partial k is seeded on the k-th component of the dependencies at every radius at once, so that the partials of a row at radius ir are its derivatives with respect to the unknowns at ir. A term with a row that reads all charge states of a radius (the radiated power, the conductivity, the energy exchange with the ions) assembles its block this way, so that the coloring of the automatic Jacobian does not need a color per charge state. The tag is fixed, so the code compiles once per number of local columns.

The evaluation runs on a state of the dependencies and the row alone, in a layout of these radial unknowns on the same radial grid (local_layout), so that its cost does not grow with the kinetic unknowns of the system.

coupling(term)

:sparse (default) when the Jacobian block of the EquationTerm term is stored, with its structural nonzeros listed by sparsity!. :matrix_free when it is applied through jacobian_vector_product!, as for the knock-on term, whose block is dense per radius.

dependencies(term)

Names of the unknowns that the EquationTerm term reads, as a collection of Symbols. Every term defines it.

has_low_rank_jacobian(term)

Whether low_rank_jacobian returns the whole Jacobian of term (false by default). A direct solver takes such matrix-free terms by a Woodbury update.

jacobian_assembly(term)

How the stored Jacobian block of a :sparse term is filled: :automatic (default) by colored forward-mode AD of rhs!, or :direct by the term's own add_jacobian!.

The colored AD needs one pass of the sparse terms per color, and two columns that share a row need different colors. A row that holds every cell of a radius (the density or current moment of a distribution, the loss through the edge) would therefore take one color per cell. Such terms assemble their block directly, and the coloring sees the other terms alone. A direct term is tested against the AD of its rhs!.

jacobian_correction_product!(out, term, u, v, layout, row, t)

For a :sparse term with jacobian_kind(term) === :approximate: add (J_exact - J_stored) v to out, the difference between the exact Jacobian-vector product of the term's residual and the product of its stored block, both at the state u. Nothing for exact terms.

jacobian_kind(term)

:exact (default) when the stored Jacobian block of a :sparse term is the derivative of its residual, :approximate when it is a cheaper matrix (a frozen limiter, another scheme's Jacobian, FokkerPlanckTerm). With an approximate block the Newton iteration on the stored matrix converges linearly (a quasi-Newton method), which NewtonSolver accounts for; GMRES(operator=:exact) restores the exact Newton operator through jacobian_correction_product!.

jacobian_vector_product!(out, term, u, v, layout, row, t)

Add the product (∂ rhs_term / ∂u) v of a :matrix_free EquationTerm (see coupling) to out.

kinetic_array(u, layout, id)

LayoutKineticArray of kinetic unknown id in state u; the unknown's grid must have the same cell counts at all radii.

low_rank_jacobian(term, u, layout, row, t)

For a :matrix_free term whose Jacobian at the state u is a sum of low-rank blocks: the vector of LowRankBlocks. nothing (the default) for any other term. The block preconditioners of GMRES and the direct solvers SparseLU and MKLPardiso add the blocks by a Woodbury update, so the term needs no stored Jacobian.

prepare_step!(term, u_old, layout, row, t)

Called with the state at the start of the step before every evaluation of the residual (residual!), and with the state itself in the ODE form (rhs!). fill_initial_values! calls it for the terms of the rows it fills, with the state and the time it is given, before the fill. The construction of an EquationSystem calls it with the initial state for a KineticEnergyLossTerm of a source with lagged = :all that was never prepared. A term with a part lagged to the previous time level inside an implicit term (the knock-on term with treatment = :explicit_primary) evaluates that part here. Nothing by default.

A term with a lagged part, evaluated on its own outside a system, takes that part from the state of its last preparation, whatever state is passed. Before any preparation, the knock-on term with a lagged primary and the KineticEnergyLossTerm of a source with lagged = :all throw an ArgumentError; a FokkerPlanckTerm with lagged = :all and a FieldParticleTerm with a plasma of unknowns keep the coefficients of the operator they were built with, and a Fokker-Planck term with lagged ions takes the ions of the state.

rhs!(out, term, u, layout, row, t)

Add the EquationTerm term to out in the rows of the unknown with id row, as a pure function of the state u and the time t. A :sparse term (see coupling) must accept state vectors of dual numbers, since its Jacobian comes from automatic differentiation.

sparsity!(entries, term, layout, row)

Push the (row index, column index) pairs of the structural nonzeros of the Jacobian block of a :sparse EquationTerm (see coupling) to entries.

treatment(term)

:implicit (default) or :explicit. An explicit EquationTerm is evaluated at the previous time level and contributes nothing to the Jacobian.

validate(term, layout, row[, t])

Check term against the layout layout of its system, in the equation of the unknown row, and throw an ArgumentError when it does not fit; run when an EquationSystem is built, with its initial time t for terms that evaluate a function of time. The default accepts every term.

local_layout(layout, ids)

The layout of the radial unknowns ids of layout alone, on its radial grid and in the order of ids, without their lower limits, whose references may be missing. With add_local_jacobian! a term evaluates on the state of this layout.

LinearTerm(unknown; coefficient=1)

coefficient times the unknown unknown, added entry by entry to the row unknown, which must have the same kind and sizes (e.g. -1 times an unknown in its own algebraic equation, or a decay rate). coefficient is a number, or one number per radius for a radial or kinetic unknown: a conductivity in Ohm's law 0 = σ E − j. The vector is referenced, so an update of the run may change its values between the intervals. The term takes the unit of its row, so coefficient has the unit of the row over that of unknown.

PrescribedValueTerm(values)

A prescribed value added to the row of a radial unknown with one component or of a scalar unknown: a number, one number per radius, or a function of the time [s] that returns either. The values take the unit of the row. For a radial unknown the function may also return a callable of r, such as an interpolator: it is evaluated on the cell centres once per time, and a profile on points of its own that is constant in time is t -> profile. With LinearTerm(name; coefficient = -1) in an algebraic equation it holds the unknown at the prescribed value, 0 = value(t) - unknown, which makes a prescribed profile an unknown of the system: a prescribed electric field E(r, t) that the kinetic terms read from the state.

CurrentMomentTerm(kinetic, grid; coefficient=1, pitch_rule=:linear, region=nothing)

coefficient times the current density [A/m² for f in m⁻³] of the kinetic unknown kinetic on grid (current_density; pitch_rule = :cell_center is DREAM's j_hot moment), added to a radial row unknown with one component (the ohmic or hot-electron current defined by an algebraic equation). The pitch weights are built per radius when the term is created, so different momentum grids at different radii work; grid must be the unknown's grid (the same radii and momentum faces, checked when the system is built).

With region, a MomentumRegion, the moment covers that part of the momentum grid, with the limit of the region at the time of the evaluation: the current of the hot electrons above the threshold of the hot region.

The term keeps the positions of its row in the stored Jacobian between calls, so its add_jacobian! is not safe to call concurrently.

DensityMomentTerm(kinetic; coefficient=1)
DensityMomentTerm(kinetic, grid; coefficient=1, region=nothing)

coefficient times the density n(r) = (1/V') Σ 𝒱' Δp Δξ f [m⁻³ for f in m⁻³] of the kinetic unknown kinetic, added to a radial row unknown with one component (e.g. a hot-electron density defined by an algebraic equation). Works with different momentum grids at different radii.

With region, a MomentumRegion, the moment covers that part of the momentum grid: the hot electrons above the threshold of the hot region. The limit of the region at the time of the evaluation counts, and grid is the grid of the unknown. A prescribed limit makes the term linear in the distribution. A limit that follows the temperature unknown (ThermalMomentum with the unknown's name) makes it nonlinear, and its Jacobian holds the column of the temperature through the threshold.

The term keeps the positions of its row in the stored Jacobian between calls, so its add_jacobian! is not safe to call concurrently.

RegionEnvelopes(region, grid)

The envelopes of a region on a kinetic grid, per radius, kept by a term. A limit that depends on time is evaluated once per time, and the envelopes are rebuilt when the limit of a radius changes. The evaluations of a term share them, so a term with a region allocates nothing between the changes.

FokkerPlanckTerm(unknown, operator::AdvectionDiffusion; jacobian=:exact, electric_field=nothing,
                 plasma=nothing, settings=nothing, synchrotron=nothing, lagged=:ions)

Advection-diffusion (Fokker-Planck) operator applied to the kinetic unknown unknown, added to the equation of the row unknown (normally the same). Stored sparse, with the 3×3 (p, ξ) stencil and the second neighbours of the four-point schemes. The corners of the 3×3 block are stored when the operator declares cross diffusion, and the radial neighbours when it declares radial transport: the cell of the same pitch at each neighbouring radius, or with different pitch grids the cells that overlap it (stencil_pattern). Needs the same cell counts and p faces at all radii.

electric_field names the radial unknown of the electric field [V/m] when the field is an unknown of the system. The advection of every face is then that of the operator, which holds the collisions, plus the field of the radius times the coefficients of a unit field (set_unit_field_coefficients!), inside the scheme. The operator must not hold the field of a prescribed plasma as well (fokker_planck_operator(...; field_in_coefficients = false)). The term is bilinear in the distribution and the field, and its stored block holds the column of the field of each radius, which is DREAM's Jacobian of its ElectricFieldTerm.

plasma is the source of the plasma that the collision coefficients follow: a PlasmaSource, a PlasmaState or a function of time that returns one (PrescribedPlasma). The coefficients of operator are rebuilt from it with settings and synchrotron (set_fokker_planck_coefficients!), and the scheme and the boundary conditions of operator are kept. Without them the term takes those that built operator (fokker_planck_operator), and CollisionSettings() with the synchrotron loss for an operator built otherwise (on by default, a choice of this package; see the docs page Choices of this package). The rebuild happens when the plasma values of any radius differ from those of the last rebuild (update_key!), not at every new time, and rebuilds all radii. Without plasma the coefficients of operator are fixed. prescribed_fokker_planck_term builds such a term from a function of time. The rebuild writes no radial transport and no cross diffusion, so a term with plasma clears both declarations of operator (AdvectionDiffusion), and its stored block couples neither the radii nor the corners of the 3×3 block.

With a plasma that reads unknowns (PlasmaFromUnknowns) the term depends on them and is not linear. lagged names the part of that plasma taken at the start of the step; the rest is taken at the state of the Newton iterate, as DREAM does, and the stored Jacobian holds its columns (add_plasma_columns!) next to the block in f. With the ion densities lagged and the temperature and the density implicit, the stored Jacobian has the columns of the temperature and the density and none of the charge states, and the residual is consistent within the step. DREAM's disruption example drops these columns from its Jacobian alone (enableIonJacobian(False), examples/Disruptions/generate.py:243).

• lagged = :all: the whole plasma at the start of the step, set in prepare_step!, with the prescribed part of the plasma at the time of the evaluation. The term is linear in f within a step and its block is exact. In the ODE form (rhs! of the system, the SciML interface) prepare_step! sees the state itself, so the term is implicit there. A term evaluated on its own keeps the coefficients of its last preparation, or those of operator when it was never prepared.
• lagged = :ions (the default, a choice of this package; see the docs page Choices of this package): the ion densities of the plasma when it reads them, no group otherwise. A term that lags them is of first order under BDF2, and ode_function passes no Jacobian for a system with it (lags_in_ode_form). The kinetic ionization and a corrected conductivity read the charge states of the iterate, so within a step the ions enter the collisions of f and those terms at two states, a difference of first order in the step. A term evaluated on its own takes the ion densities of its last prepare_step!, or those of the state when it was never prepared.
• lagged = (): nothing lagged; the columns of the charge states are part of the plasma columns.
• lagged = (:ion_densities,): the ion densities by name, as :ions for that unknown.

jacobian selects what the stored Jacobian block holds for a flux-limiter scheme (DREAM's ad_jac setting of AdvectionInterpolation). The residual is always the operator's own, so every mode converges to the same solution:

• :exact: the Jacobian of the residual by colored forward-mode AD, including the limiter's derivative ψ′ (DREAM AD_INTERP_JACOBIAN_FULL, which carries ψ′ only for the limiters DREAM counts as smooth);
• :frozen_limiter: the limiter weights frozen at the iterate and differentiated as constants (apply!(...; limiter_state), DREAM AD_INTERP_JACOBIAN_LINEAR);
• a linear scheme name (:upwind, :central, :exponential_fitting, :quick, :upwind_second_order): the Jacobian of the operator with that scheme in p and ξ (with_scheme; :upwind is DREAM AD_INTERP_JACOBIAN_UPWIND), with that scheme's stencil (3-point for the two-point schemes).

For a scheme linear in f the frozen limiter is the exact Jacobian, and another linear scheme's is still a different matrix (upwind's for a central residual, the approximation that rescues some DREAM runs at high momentum resolution with central interpolation; see the docs page Comparison with DREAM). With an approximate block the term declares jacobian_kind(term) === :approximate: the Newton iteration on the stored matrix is a quasi-Newton one that converges linearly, or, with GMRES(operator=:exact), the stored block only preconditions the exact Newton operator (jacobian_correction_product!, one dual-number pass of the operator per product).

add_plasma_columns!(J, term, u, layout, row, t)

Add the plasma columns of a term with has_plasma_columns to J at the state u and the time t [s]. A FokkerPlanckTerm builds its coefficients from the plasma with dual numbers in the temperature and the density of every radius (dual_plasma) and applies them to f: the partials are the derivatives of the residual through the scheme's own face weights, the Péclet number of exponential fitting included. The scheme is the operator's own, as is the residual, whatever jacobian selects for the block in f. The coefficients in dual numbers are rebuilt when the plasma values of any radius differ from those of their last rebuild, as the coefficients of the residual (refresh_coefficients!), so the terms that read them share one rebuild per state.

has_plasma_columns(term)

Whether the stored Jacobian holds columns of the plasma unknowns of a term: the derivative of the term with respect to the temperature and the density of each radius, which the colored AD of the term does not see, since its coefficients hold numbers. Such a term lists them with plasma_column_sparsity! and adds them with add_plasma_columns!. True for a FokkerPlanckTerm whose plasma reads unknowns and whose coefficients are implicit. The terms that read the coefficients of such a term (MomentumBoundaryFluxTerm, MomentumInterfaceTerm, InterfaceEnergyFluxTerm, InterfaceParticleFluxTerm) pass its plasma columns on.

lags_in_ode_form(term)

Whether a term takes a part of its residual at the state of prepare_step!, the start of the step: a knock-on term with a lagged primary, a Fokker-Planck term with lagged coefficients or lagged ions, the energy loss of such a term, a field-particle term whose plasma reads unknowns. In the ODE form (rhs! of the system) that state is the state itself, so the stored Jacobian, which treats the lagged part as fixed, is not the Jacobian of the ODE, and ode_function passes none.

plasma_column_sparsity!(entries, term, layout, row)

Push the entries of the plasma columns of a term with has_plasma_columns: every cell of the row unknown at a radius against each plasma unknown of that radius.

prescribed_fokker_planck_term(unknown, operator::AdvectionDiffusion, plasma_at;
                              settings=nothing, synchrotron=nothing, jacobian=:exact,
                              electric_field=nothing)

A FokkerPlanckTerm whose collision coefficients follow a plasma prescribed in time (such as a cooling plasma with T(t) prescribed): plasma_at(t) returns the PlasmaState at the time t [s]. The same as FokkerPlanckTerm(unknown, operator; plasma = plasma_at, settings, synchrotron, jacobian, electric_field), which this function returns. The coefficients are rebuilt when the plasma values change (refresh_coefficients!).

With electric_field, the name of the unknown of the field, the field of the prescribed plasma is left out of the coefficients, and the term takes the field of the state.

has_plasma_source(term)

Whether the coefficients of a term follow a plasma source (refresh_coefficients!).

refresh_coefficients!(term, u, layout, t)

Bring the coefficients of a term up to date with its plasma at the state u and the time t, and return the operator that holds them. A FokkerPlanckTerm with a plasma source rebuilds its operator when the plasma values of any radius changed since the last rebuild. A term without a plasma source returns its operator unchanged, and so does a term with lagged = :all, which prepare_step! refreshes. u and layout may be nothing for a prescribed plasma.

KnockOnTerm(primary, operator::KnockOnOperator; target_density, treatment=:explicit_primary)

Knock-on source [m⁻³ s⁻¹], the rate of the distribution of the kinetic row unknown, produced by the kinetic unknown primary on that unknown (operator maps the primary grid to the row unknown's grid, per radius, see apply_at_radius!). target_density is the target (total electron) density [m⁻³]: the name of a radial unknown, or a vector with one value per radius. Works with different momentum grids at different radii.

treatment:

• :explicit_primary (default, DREAM's Boltzmann operator MollerBoltzmannOperator): the primary distribution at the previous time level, the target density at the new one. The term is then the target density times a source fixed for the step, evaluated once per start state of a step (prepare_step!) and again after a change of the cutoff (set_cutoff!), and its Jacobian is the column of the target density at each radius. With a prescribed target density it is the same as :explicit.
• :explicit: the whole term at the previous time level, target density included.
• :implicit: the whole term at the new time level, matrix-free, since it couples all cells of a radius. Every Krylov iteration of GMRES applies the operator once, and the direct solvers (SparseLU, MKLPardiso) do not take it.

The explicit forms are the cheaper step; the implicit form is the more accurate one per step. The explicit default follows DREAM.

MomentumInterfaceTerm(lower, upper, lower_term, upper_operator;
                      pitch_mapping=:overlap, scheme=<scheme of the lower operator>)

Particle flux through the momentum face shared by the kinetic unknowns lower and upper: the last p face of the grid of lower (the hot-tail grid) is the first p face of the grid of upper (the runaway grid) at every radius. The pitch grids of the two may differ. The same term goes into both equations. In the equation of lower it removes the flux from the last momentum cells, in the equation of upper it adds the flux to the first momentum cells, and both are computed from the same numbers, so no particle is lost or created at the interface.

lower_term is the Fokker-Planck term of lower (FokkerPlanckTerm or prescribed_fokker_planck_term), or its operator. A prescribed term is passed as the term: a system refuses the bare operator of a prescribed term of the system, whose coefficients belong to the time of its last rebuild. The flux takes the coefficients Aᵖ and Dᵖᵖ of that operator at its last p face, as DREAM's PXiExternalKineticKinetic does, and the two operators add no flux of their own through the face: the lower one has p_upper_boundary = :interface, upper_operator has p_lower_boundary = :interface. When lower_term takes the electric field as an unknown, the advection at the face holds the field of the state, and the term depends on the field as well. The cross term Dᵖᵠ at the face is not included (none in DREAM either, and zero for the collision operator). When the coefficients of lower_term follow the plasma unknowns (has_plasma_columns), so does the flux: its columns of the temperature, the density, and the ion densities unless the lower term lags them, come from the flux of the term's coefficients in dual numbers (add_plasma_columns!), as DREAM's PXiExternalKineticKinetic adds the derivatives of its coefficients (AddPartialJacobianContributions).

pitch_mapping:

• :overlap (Runes2P): one flux per overlap cell of the two pitch grids (PitchOverlap), from the value of the lower cell and the value of the upper cell the overlap cell lies in,

  Γ_s = flux(Aᵖ_j, Dᵖᵖ_j, f_lower[end, j], f_upper[1, J]; scheme),

  with Γ_s [m⁻³ s⁻¹] the momentum flux Γᵖ at the face, before its weight 𝒱', and flux the two-point flux of scheme between the two cell centres (:upwind, :central or :exponential_fitting, by default the p scheme of the lower operator, and :upwind when that is a four-point scheme). The couplings are those of the scheme between two neighbouring cells, so exponential fitting and upwind keep their sign properties across the interface.

• :dream: the flux of DREAM, one per pitch cell of the lower grid. The upper value is interpolated linearly in the pitch to the centre of the lower cell, the advection is first-order upwind whatever the scheme of the interior, the diffusion is the difference of the two values over the distance of the centres, and the flux is given to the upper cells by their overlap with the lower cell. Ported as written, including the rule for a lower centre beyond the last upper centre, which takes the second-to-last upper cell (the nearest one at the other end). scheme does not apply. Where the pitch grids differ, the interpolation couples an upper cell to its neighbours in the pitch with negative weights in the diffusive part, so this mapping does not keep the sign of f by construction.

Both conserve particles to round-off. Linear in the state unless the coefficients follow the field or the plasma unknowns; stored sparse. interface_flux gives the flux through the face as a rate per volume.

interface_flux(term::MomentumInterfaceTerm, u, layout, ir, t=0)

Particle flux per unit volume [m⁻³ s⁻¹] through the momentum interface at radius index ir, positive from the lower to the upper grid, for the state u at the time t [s]: the rate at which the upper grid gains particles from the lower one, in the normalization of density and momentum_flux. With the interface above the critical momentum this is the runaway rate.

LossRateTerm(unknown, rate)

Prescribed loss -ν f of the kinetic unknown unknown with the rate ν [1/s] of each radius (rate, one value per radius), added to the row unknown, which has the same grid. The exponential loss of a region with open field lines, as DREAM's loss outside the last closed flux surface. The vector is referenced: an update of the run may change its values between the intervals.

RadialEdgeLossTerm(kinetic, operator; coefficient=1)

coefficient times the rate [m⁻¹ s⁻¹], the electrons per second per metre of the cylinder, divided by 2π, at which the kinetic unknown kinetic loses particles through the edge of the plasma by the radial transport of operator (radial_particle_flux at the last radial face), added to a scalar row unknown. With TimeDerivative() the row unknown is the content lost so far [m⁻¹], in the normalization of total_particles, integrated with the same time levels as the distribution, so that the content of the distribution plus the lost particles stays constant to round-off when nothing else changes the content. The term reads the coefficients of operator as they are, also when an update changes them during the run.

MaxwellianSourceTerm(kinetic, grid, temperature; target_density, rate,
                     bulk_momentum=Inf)

Source with the shape of a non-relativistic Maxwellian exp(-p²/2Θ), Θ = T/mₑc², at temperature [eV] (DREAM's ParticleSourceTerm uses the relativistic one), whose magnitude drives the density of the kinetic unknown kinetic back to target_density [m⁻³] (a number, or one value per radius) at the rate rate [1/s]:

S(p, ξ) = rate (n_target - n(f)) f_M(p) / n_M,

with f_M the Maxwellian on grid and n_M its density, so the source adds rate (n_target - n(f)) particles per second and per unit volume and nothing else. kinetic and the row unknown both live on grid. It keeps a thermal bulk from draining into the runaway region, which is what the published Dreicer rates assume; with rate = 1/Δt a step restores the particles the previous step lost.

n(f) counts the cells below bulk_momentum (in units of mₑc), and n_target is the density of the bulk there. The bound matters: with the whole grid counted, particles moving from the bulk into the runaway tail leave the total unchanged, the source sees no deficit, and the bulk still drains. Set it below the critical momentum and above the thermal momentum.

The term is :explicit (evaluated at the previous time level), so the density moment it reads adds nothing to the Jacobian and the sparsity stays local. This is a device for comparisons with the Dreicer rates, not the self-consistent particle source: DREAM's ParticleSourceTerm fixes the magnitude from the cold-electron density through quasi-neutrality (ParticleSourceTerm).

MomentumBoundaryFluxTerm(kinetic, source; weight=:density, coefficient=1)

coefficient times the flux of the kinetic unknown kinetic through the upper momentum boundary of its grid, per volume, added to the row of a radial unknown with one component: the electrons that leave the grid at p_max. source holds the coefficients of the flux: the Fokker-Planck term of kinetic (FokkerPlanckTerm or prescribed_fokker_planck_term) or its operator. The flux is that of the boundary condition of the operator at p_max, so the term and the operator count the same particles. With a term that takes the electric field as an unknown, the flux holds the field of the state.

• weight = :density: the particles per volume and second [m⁻³ s⁻¹], (1/V') Σ_j Δξ_j 𝒱' Γᵖ. With TimeDerivative() the row unknown is the density of the electrons above the grid, or a fluid runaway density fed by the hot-tail grid.
• weight = :current: the flux weighted with e c ξ, the rate [A m⁻² s⁻¹] of the current density that the electrons carry when they leave, at the speed of light. With TimeDerivative() the row unknown is the current density of the electrons above the grid.
• weight = :energy: the flux weighted with the kinetic energy mₑc²(γ − 1) of the last momentum cell [W/m³], the energy that the discrete energy moment of the grid loses with the electrons; the energy at the face itself differs by about Δp/(2p). With coefficient = −1 in the row of the temperature it is the energy that leaves the total energy of the electrons on the hot-tail grid (total_energy_equation).

The entries of the row are assembled directly, from the derivative of the flux of each pitch cell with respect to the cells at the boundary. With a Fokker-Planck term whose plasma reads unknowns implicitly (has_plasma_columns) the flux follows them as well: its columns of the temperature, the density and the charge states come from the flux of the term's coefficients in dual numbers, through the weights of the scheme of the boundary (add_plasma_columns!).

ParticleSourceTerm(source, grid, temperature; shape=:relativistic_maxwellian)

Particle source of a kinetic unknown on grid with the magnitude of the radial unknown source [1/(m³ s)]: S(r) g(p), with the shape g normalized on the grid to one particle per volume, so that the term adds S particles per volume and second and nothing else (DREAM's ParticleSourceTerm). The magnitude is an unknown whose equation is the constraint it enforces, e.g. 0 = n(f) − n_cold − n_hot with the kind Constraint: the source moves electrons between the cold population and the distribution as electrons run away or the free-electron density changes.

shape:

• :relativistic_maxwellian: exp(-(γ − 1)/Θ) with Θ = T/mₑc² (DREAM's default);
• :maxwellian: the non-relativistic exp(-p²/2Θ) of MaxwellianSourceTerm;
• :lowest_cell: the momentum cell at p = 0 alone (DREAM's delta shape). It does not depend on the temperature, and the stored Jacobian holds the first momentum cell of every pitch cell in the column of the magnitude and no column of the temperature. A temperature that is an unknown stays among the term's dependencies, and the term is not linear (is_linear), although its column of the temperature is zero.

temperature is the cold-electron temperature in eV at each radius, or a function of time that returns the temperatures or a PlasmaState. The shape follows it in time and is rebuilt when the temperature of a radius changes. It is prescribed data: the term is linear, and its Jacobian entries are assembled directly.

temperature may also be the name of the radial unknown of the temperature. The term is then not linear, and its Jacobian holds the column of the temperature, S ∂g/∂T, with the derivative of the normalization on the grid, which DREAM's column leaves out (ParticleSourceTerm.cpp). The momentum integral of that column vanishes, as the source adds S particles at any temperature.

source_shape(term::ParticleSourceTerm, ir, t=0)

Shape of the source at the radius index ir and the time t [s]: the particles per volume, second and phase-space volume that a magnitude of 1/(m³ s) adds to each cell. A term whose temperature is an unknown needs the state, source_shape(term, ir, u, layout, t).

ExternalAvalancheTerm(plasma; p_cutoff, runaway_density=:runaway_density,
                      total_density=:total_electron_density, electric_field=:electric_field,
                      temperature=:cold_electron_temperature, density=:cold_electron_density,
                      ions=nothing, settings=CollisionSettings(), exponent=-2)

DREAM's ExternalAvalancheTerm (src/Equations/Fluid/ExternalAvalancheTerm.cpp, placed by src/Settings/RunawaySourceTerms.cpp:36-58) in the equation of the fluid runaway density n_re above a hot-tail grid without a runaway grid: dn_re/dt = Γ n_re [m⁻³ s⁻¹] with the rate Γ [s⁻¹]

\[Γ = \left(Γ_\mathrm{ext}^{a} + Γ_\mathrm{fluid}^{a}\right)^{1/a} \text{ for } Γ_\mathrm{fluid} > 0, \qquad Γ = Γ_\mathrm{fluid} \text{ otherwise},\]

a = exponent = −2 (DREAM's value), the Rosenbluth-Putvinski rate Γ_ext [s⁻¹] of the knock-ons above p_cutoff (units of mₑc), Γ_ext = n_tot rosenbluth_putvinski_rate(p_cutoff), and DREAM's fluid growth rate [s⁻¹] with the effective critical field, fluid_avalanche_growth_rate. The match makes the external rate, positive at every field, follow the fluid rate below the effective critical field, where it decays. DREAM takes p_cutoff as the upper momentum face of the hot-tail grid, p_max (RunawaySourceTerms.cpp:44-46), not the cutoff of the kinetic source (RosenbluthPutvinskiSourceTerm), and so must the case. One p_cutoff serves every radius, as in DREAM ("assume same momentum grid at all radii", RunawaySourceTerms.cpp:41): a grid whose p_max differs between radii needs a term per radius.

The plasma of the rate is that of plasma (a PlasmaState or a function of time: the magnetic field, the composition) with the temperature temperature, the cold density density and, with ions, the charge states of that radial unknown, as the other fluid terms read them; settings gives the screening, the Coulomb logarithm and the tables of DREAM's RunawayFluid (runaway_fluid_settings). The Jacobian is exact: the columns of n_re, n_tot, the field, the temperature, the density and the charge states by forward differentiation through the solutions of the effective critical field and of p⋆ (one Newton correction in dual numbers each). DREAM's own Jacobian takes the columns of the field and of n_tot alone, with the approximation Γ ∝ |E| − E_ceff (RunawayFluid.cpp:881-907). The solution of the effective critical field takes about 2.4 ms per radius on one core (the D and Ar of DREAM's disruption example; 2.5 ms with the derivatives). The rate of a radius is evaluated once per state: the rows that hold the term (the runaway density and the particle source of a hot-tail grid with fluid runaways) and the residual and the Jacobian at one state share the evaluation, and the Jacobian takes every partial in one pass of forward differentiation.

RosenbluthPutvinskiSourceTerm(kinetic, grid; runaway_density=:runaway_density,
                              total_density=:total_electron_density,
                              electric_field=:electric_field, p_cutoff=0.01)

DREAM's kinetic Rosenbluth-Putvinski avalanche source (AVALANCHE_MODE_KINETIC, src/Equations/Kinetic/AvalancheSourceRP.cpp, placed by src/Settings/Equations/f_hot.cpp:116-127) in the equation of the kinetic unknown kinetic on grid: the knock-on electrons that a fluid runaway density n_re above the grid creates on the electrons of the total density n_tot (free and bound), at every momentum cell above p_cutoff (units of mₑc) and at the pitch ξ* = ±√((γ − 1)/(γ + 1)) that the momentum of a knock-on from a primary at the speed of light fixes, in the cell that holds it:

\[S_{ij} = n_\mathrm{re} n_\mathrm{tot}\, r_0^2 c\, \frac{1}{Δp_i}\left(\frac{1}{γ_{i−1/2}^* − 1} − \frac{1}{γ_{i+1/2} − 1}\right) \frac{2π V'}{\mathcal{V}'_{ij} Δξ_j},\]

with γ*_{i−1/2} at the larger of the lower face and p_cutoff and Δp_i the whole cell. The source S_ij [m⁻³ s⁻¹] is the rate of the distribution [m⁻³] in the cell (i, j). Summed over the grid it is rosenbluth_putvinski_rate(p_cutoff, p_max) times n_re n_tot exactly. The pitch is that of DREAM in a cylinder (B_max = B_min, fvm/Grid/FluxSurfaceAverager.avalancheDelta.cpp:114-140): ξ* takes the sign of the electric field of each radius (RP_SOURCE_PITCH_ADAPTIVE), or + without electric_field, and a ξ* on a cell face gives no cell (DREAM's strict inequalities). The term is bilinear in n_re and n_tot; the sign of the field enters without a derivative, as in DREAM. The Jacobian is assembled directly.

The source creates electrons, which the particle source of the grid takes out of the cold electrons: DREAM's row of the particle source, S = dn_free/dt − internal − external, holds the internal number with the factor −1 (TotalElectronDensityFromKineticAvalanche, f_hot.cpp:244-255, 362-374). Written as dn_free/dt = S + internal + external (TimeDerivative(:free_electron_density)), the term is ProductTerm(total_density, runaway_density; coefficient = +rosenbluth_putvinski_rate(p_cutoff, p_max)). The value p_cutoff = 0.01 is that of DREAM's disruption example (examples/Disruptions/generate.py:329).

rosenbluth_putvinski_rate(p_lower, p_upper=Inf)

The rate of the Rosenbluth-Putvinski avalanche source between the momenta p_lower and p_upper (units of mₑc) per runaway and per electron of the total density, in m³/s:

\[\frac{1}{n_\mathrm{re} n_\mathrm{tot}} \int_{p_\mathrm{lower}}^{p_\mathrm{upper}} S\,\mathrm{d}^3p = 2π r_0^2 c \left(\frac{1}{γ_\mathrm{lower} − 1} − \frac{1}{γ_\mathrm{upper} − 1}\right),\]

the knock-ons of a runaway at the speed of light on electrons at rest: the Møller cross section in the limit of an infinitely fast primary, 2π r₀²/(γ − 1)² per unit γ (M. N. Rosenbluth and S. V. Putvinski, Nucl. Fusion 37, 1355 (1997), doi:10.1088/0029-5515/37/10/I03, eq. (4)). DREAM's AvalancheSourceRP::EvaluateNormalizedTotalKnockOnNumber (src/Equations/Kinetic/AvalancheSourceRP.cpp:96-115), with its prefactor e⁴/((4πε₀ mₑ c)² c) = r₀² c. Infinite at p_lower = 0.

AvalancheRateCache

The last evaluation of the matched rate of ExternalAvalancheTerm at one radius, shared by the rows that hold the term and by the residual and the Jacobian at one state: the local values x and the time t of the entry, the value of the rate times n_re, and its derivatives by x once the Jacobian needed them. One entry per radius, so that a loop over radii in threads writes separate entries. The entries take a state in Float64; a state in another number type is evaluated without the cache.

THERMAL_ENERGY_FACTOR

3/2 e [J/eV]: the thermal energy density W = 3/2 e n T [J/m³] of electrons with the density n [m⁻³] and the temperature T [eV] (DREAM's ElectronHeatTerm).

EnergyMomentTerm(kinetic, grid; coefficient=1, region=nothing)

coefficient times the kinetic energy density W(r) = (1/V') Σ 𝒱' Δp Δξ mₑc²(γ − 1) f [J/m³ for f in m⁻³] of the kinetic unknown kinetic on grid, added to a radial row unknown with one component, with γ at the momentum cell centres (DREAM's KineticEnergyFromDistributionFunction). With region, a MomentumRegion, the moment covers that part of the grid, as DensityMomentTerm: the energy of the hot electrons above the threshold, whose limit may follow the temperature unknown. Linear in the distribution; the column of the temperature is that of the envelope.

The term keeps the positions of its row in the stored Jacobian between calls, so its add_jacobian! is not safe to call concurrently.

ProductTerm(first, second; coefficient=1)

coefficient times the product of the radial unknowns first and second, each with one component, at every radius, added to the row of a radial unknown with one component (DREAM's DiagonalQuadraticTerm). Not linear in the state. The term takes the unit of its row, so coefficient has the unit of the row over the product of the units of first and second.

RunawayFluidEnergyTransferTerm(plasma; runaway_density=:runaway_density,
                               temperature=:cold_electron_temperature,
                               density=:cold_electron_density,
                               coulomb_logarithm=:energy_dependent, ions=nothing)

DREAM's heating of the cold electrons by the friction of fluid runaways when there is no runaway grid (include/DREAM/Equations/Fluid/CollisionalEnergyTransferREFluidTerm.hpp, the momentum 20 at line 26, used in src/Settings/Equations/T_cold.cpp), in W/m³ for the row of the temperature:

\[e c E_c n_\mathrm{re} = 4π r_0^2 m_e c^3 \ln Λ_{ee}(p = 20)\, n_\mathrm{cold} n_\mathrm{re},\]

the friction at the speed of light on runaways of momentum 20 mₑc. The Coulomb logarithm of electrons takes the temperature unknown and the free-electron density of the prescribed ions in plasma (a PlasmaState or a function of time), or with ions of that radial unknown of the charge states; coulomb_logarithm is the model of coulomb_logarithm. DREAM evaluates the logarithm with the model of its collision settings, whose code default is the thermal one, independent of p; the default :energy_dependent here is the setting of DREAM's disruption example (examples/Disruptions/generate.py:145). The term is the DREAM-compatible mode of the energy equation.

set_thermal_energy!(u, layout; energy=:cold_thermal_energy,
                    temperature=:cold_electron_temperature, density=:cold_electron_density)

Set the energy of the cold electrons in the state u from its temperature and density, W = 3/2 e n T [J/m³] at every radius. DREAM's order of the initial values: the temperature from the input, the density from its row, and then the energy (T_cold.cpp, n_cold.cpp). The density must already hold its value.

thermal_energy_equation(; energy=:cold_thermal_energy, temperature=:cold_electron_temperature,
                        density=:cold_electron_density)

The algebraic equation of the energy of the cold electrons, 0 = 3/2 e n T − W with W in J/m³ and T in eV (DREAM's W_cold). The row of the temperature holds the time derivative of the energy, Equation(temperature, TimeDerivative(energy), ...), with the heating and loss terms in W/m³, as in DREAM, where the energy is the unknown of the time derivative and the temperature that of its own row.

The energy takes its initial value from the row (initial_value = :from_equation): integrate sets it from the temperature and the density of the initial state, and the value of the energy in that state is not read.

total_energy_equation(kinetic, grid; region, total_energy=:total_electron_energy,
                      energy=:cold_thermal_energy)

The algebraic equation of the total energy W_tot [J/m³] of the electrons in the fully kinetic mode:

\[W_\mathrm{tot} = W_\mathrm{cold} + \sum_{p > p_0} m_e c^2 (γ − 1) f\, ΔV,\]

the energy energy of the cold electrons and that of the kinetic unknown kinetic on grid in the hot region region (EnergyMomentTerm), the region of the hot density. The row of the temperature then holds the time derivative of total_energy (TimeDerivative(total_energy)) with the channels that leave or enter the electrons: the work of the field on all of them (E j_ohm + E j_hot in the fully kinetic mode, where j_ohm = j(f) − j_hot + E(σ − σ_num) already leaves out the hot current, so the two products count each electron once), the radiated power of the ions, the heat transport and the exchange with the ions. The collisional transfer from the hot electrons to the cold ones, the electrons that cross the threshold, and the threshold that moves with the temperature are then exact for any time stepper, since the time derivative acts on the energy of the hot region itself. DREAM's fully kinetic mode has none of them (T_cold.cpp:313-321).

Every other change of the energy of the hot region is charged to the cold electrons, since W_cold = W_tot − W_hot. The row must therefore list each term of the equation of kinetic that moves energy out of the electrons: KineticEnergyLossTerm with part = :sinks for the bremsstrahlung, the synchrotron loss and the friction on bound electrons, MomentumBoundaryFluxTerm(...; weight = :energy) for an outflow at p_max, and InterfaceEnergyFluxTerm for the runaway grid. The system refuses a row that misses one of them, and it refuses radial transport of kinetic and knock-on secondaries on its grid, which have no such term: their energy would come out of the cold electrons of the radius where it is moved. The sinks take no flux through p_max, whose energy is the boundary term's alone. The sinks term is required by the settings of the collisions (synchrotron, bremsstrahlung, partial or no screening on), not by its value: at B = 0 a synchrotron loss of zero still needs the term, so that the row stays closed when the field of the plasma changes.

total_energy takes its initial value from the row, after the energy of the cold electrons.

RadiatedPowerTerm(plasma; radiation=:recombination, temperature=:cold_electron_temperature,
                  density=:cold_electron_density, ions=nothing,
                  kinetic_ionization=nothing, extrapolation=:dream,
                  bremsstrahlung_corrections=:gould,
                  directory=adas_cache_directory(), initial_time=0.0)

The power that the cold electrons lose to the ions, −P [W/m³] in the row of their temperature (which holds the time derivative of their energy), at every radius:

\[P = n_e \sum_{i, Z_0} n_{i, Z_0} \left(L_{i, Z_0} + B_{i, Z_0}\right)\]

with the coefficients L and B [W m³] of each charge state, the electron density n_e and temperature from the unknowns density and temperature, and the densities of the charge states n_{i,Z_0} from the ions of plasma, a PlasmaState or a function of time that returns one, or from the radial unknown ions of the charge-state densities of the species of plasma (ion_densities_unknown). The species are named by their element ("Ar", "D"), whose ADF11 tables are read from directory when the term is built (read_adas). A plasma that is a function of time gives these species at the time initial_time [s]. With ions the term assembles its Jacobian with the columns of the charge states (add_local_jacobian!). With kinetic_ionization, the KineticIonizationTerm of the split of IonRateTerm, the binding energy of the ionization is that of the cold electrons below its threshold, at the fluid coefficient of fluid_ionization_coefficient.

• radiation = :recombination (the default, a choice of this package; see the docs page Choices of this package): the line radiation L = PLT + PRB with the recombination radiation and bremsstrahlung, the binding energy that ionization takes, I(Z_0) SCD, and the one recombination returns, −I(Z_0 − 1) ACD, with the ionization energies I of NIST (ionization_energy). DREAM's mode with recombination takes I(Z_0) in the gain from its second radial cell on.
• radiation = :dream: DREAM's default, L = PLT, the ionization loss alone, and the bremsstrahlung of Gould with a relativistic correction, C √T (Σ n Z_0² + (a Σ n Z_0² + κ n_e) T/mₑc²) n_e, C = 1.69e-38 W m³ eV^{-1/2} (GOULD_BREMSSTRAHLUNG_PREFACTOR). With bremsstrahlung_corrections = :gould (default) a = 11/24 and κ = 5(44 − 3π²)/(24√2), Gould's values after his erratum and by his own eqs. (18) and (22) (GOULD_ION_CORRECTION, GOULD_ELECTRON_CORRECTION); with :dream DREAM's 19/24 and 5(44 − 3π²)/(8√2). DREAM's radiation therefore needs both radiation = :dream and bremsstrahlung_corrections = :dream; the second has no effect with :recombination. The ion sum is that of the charge states the term reads, n_e the unknown.

The coefficients are those of ADAS (H. P. Summers, The ADAS User Manual, version 2.6 (2004), http://www.adas.ac.uk, with the data of OpenADAS, https://open.adas.ac.uk), DREAM's src/Equations/Fluid/RadiatedPowerTerm.cpp. extrapolation goes to adas_coefficient; its default :dream takes DREAM's extrapolation below the tables (a choice of this package; see the docs page Choices of this package). The derivatives with respect to the temperature and the density are exact through ForwardDiff. radiated_power gives the radiation and the binding energy apart.

ohmic_heating_term(; electric_field=:electric_field,
                   current_density=:ohmic_current_density)

The ohmic heating j_ohm E [W/m³] of the cold electrons, for the row of their temperature (DREAM's OhmicHeatingTerm, whose geometric weight is 1 in a cylinder): a ProductTerm of the field electric_field [V/m] and the ohmic current density current_density [A/m²]. The heating takes the treatment of the equation that sets the current (the conductivity of a ConductivityTerm), so the two agree by construction.

radiated_power(term::RadiatedPowerTerm, u, layout, t)

The two parts of the power of RadiatedPowerTerm at the state u and the time t [s], per radius [W/m³]: radiation, which leaves the plasma (line radiation, recombination radiation and bremsstrahlung), and binding, the binding energy that the ions take on ionization, net of what recombination returns, which stays in the plasma.

DREAM_GOULD_ELECTRON_CORRECTION

DREAM's electron-electron bremsstrahlung relative to the electron-ion one, 5(44 − 3π²)/(8√2) ≈ 6.36 times T/mₑc², eq. (41) of Gould 1980 as printed (src/Equations/Fluid/RadiatedPowerTerm.cpp:52, "e-e brems correction"); GOULD_ELECTRON_CORRECTION holds the value of eqs. (18) and (22).

DREAM_GOULD_ION_CORRECTION

DREAM's relativistic correction of the electron-ion bremsstrahlung, the 19/24 of Gould 1980 before its erratum (src/Equations/Fluid/RadiatedPowerTerm.cpp:51, "relativistic-maxwellian correction"); GOULD_ION_CORRECTION holds the corrected value.

GOULD_BREMSSTRAHLUNG_PREFACTOR

The prefactor of the bremsstrahlung power of a Maxwellian plasma in the mode :dream of RadiatedPowerTerm, (32/3) α r₀² c √(2 mₑc² e/π) [W m³ eV^-1/2], the 1.69e-38 of the NRL Plasma Formulary: R. J. Gould, "Thermal bremsstrahlung from high-temperature plasmas", Astrophys. J. 238, 1026 (1980), ADS 1980ApJ...238.1026G (erratum Astrophys. J. 243, 677 (1981)), as DREAM's src/Equations/Fluid/RadiatedPowerTerm.cpp (lines 47-52 and 174-189, the 1.69e-38 named at line 179). Gould's expansion holds for 10⁷ to 10⁹ K, about 1 to 100 keV. Neither Runes2P nor DREAM has its first-order Born correction, (ln 2 − 1/2) π (Σ Z³ n/Σ Z² n) (π Ry/kT)^{1/2} (eq. (41) with the erratum's π), which at 10 eV is 1.25 times the leading term for charge 1: at the temperatures of a disruption the formula is outside its range, where the line radiation of the impurities dominates the losses.

GOULD_ELECTRON_CORRECTION

The electron-electron bremsstrahlung of a Maxwellian plasma relative to the electron-ion one of charge 1, 5(44 − 3π²)/(24√2) ≈ 2.12 times T/mₑc²: the ratio of eqs. (22) and (18) of R. J. Gould, Astrophys. J. 238, 1026 (1980), the Born approximation of the electron-electron and the electron-ion power, (20/9)(44 − 3π²) π^{-1/2} α³Λ² m^{-1/2} n_e² (kT)^{3/2} over (32/3)(2/π)^{1/2} α³Λ² m c³ n_e Σ Z² n_Z τ^{1/2}; also eq. (58) of H. Xie, "Bremsstrahlung radiation power in fusion plasmas revisited", arXiv:2404.11540 (2024). Gould's summary, eq. (41), which says it combines eqs. (18) and (22), prints 5(44 − 3π²)/(8√2) ≈ 6.36, three times larger, and the erratum (Astrophys. J. 243, 677 (1981)) does not correct it; DREAM takes that value (DREAM_GOULD_ELECTRON_CORRECTION). The factor rests on the consistency of the paper's own equations and on Xie, not on an erratum.

GOULD_ION_CORRECTION

The first-order relativistic correction 11/24 of the electron-ion bremsstrahlung of a Maxwellian plasma, the coefficient of T/mₑc² in the power of GOULD_BREMSSTRAHLUNG_PREFACTOR: R. J. Gould, Astrophys. J. 238, 1026 (1980), eqs. (18) and (41), as corrected in the erratum, Astrophys. J. 243, 677 (1981) ("In equations (18) and (41), 19/24 should be replaced by 11/24", found by E. Haug). It is the slope in T of the Jüttner average of the cross section 4BN of H. W. Koch and J. W. Motz, Rev. Mod. Phys. 31, 920 (1959), bremsstrahlung_cross_section, which gives 11/24 to 1e-5 (test_independent_references.jl). DREAM's src/Equations/Fluid/RadiatedPowerTerm.cpp:51 keeps the 19/24 of the paper (DREAM_GOULD_ION_CORRECTION).

HeatTransportTerm(; diffusion=nothing, diffusion_at=nothing, delta_b_over_b=nothing,
                  delta_b_over_b_at=nothing, major_radius=nothing, safety_factor=1.0,
                  temperature=:cold_electron_temperature, density=:cold_electron_density,
                  edge=:zero_flux, coefficient=1)

Radial heat diffusion of the cold electrons, (1/V') ∂_r (V' χ ∂_r T) [W/m³] with T in eV, for the row of their temperature. The flux through the face between two cells is V' χ (T_k − T_{k−1})/Δr with the distance Δr of the two centres, and χ at the face takes the density and the temperature interpolated linearly between the two cells (DREAM's deltaRadialFlux). No heat crosses the first face. One of two models:

• a diffusivity D [m²/s], χ = 3/2 e D n (DREAM's src/Equations/Fluid/HeatTransportDiffusion.cpp:85-87): diffusion a radial profile (radial_values on the faces), or diffusion_at a function of the time [s] that returns one;
• the heat diffusivity of Rechester and Rosenbluth in stochastic field lines (A. B. Rechester and M. N. Rosenbluth, Phys. Rev. Lett. 40, 38 (1978), doi:10.1103/PhysRevLett.40.38; DREAM's src/Equations/Fluid/HeatTransportRechesterRosenbluth.cpp:80 for 3√(2π) c e and :105 for (1 − 5Θ/8), cylinder),

  \[χ = 3\sqrt{2π}\, e\, c\, q R_0 (δB/B)^2\, n \sqrt{Θ} (1 − 5Θ/8),\qquad Θ = T/m_ec^2,\]

  the energy moment of the kinetic diffusion π q R₀ |v_∥| (δB/B)² of add_rechester_rosenbluth_diffusion! on a Maxwellian of uniform density. The part of that moment in ∂_r n is dropped, as in DREAM, whose derivation assumes ∂n_cold/∂r = 0 because the cold electrons are not transported (doc/notes/heat-transport.tex:236-245). delta_b_over_b is a profile, or delta_b_over_b_at a function of the time [s] such as the PiecewiseConstant transport history of a RechesterRosenbluthUpdate. The same object in both keeps the kinetic and the heat transport on one level. The major radius R₀ [m] and the safety factor q are inputs, as for the kinetic operator. The default q = 1 is DREAM's, which fixes it (HeatTransportRechesterRosenbluth.cpp:98-100, marked as a TODO there).

A PiecewiseConstant level takes, in the step that ends at a switch time, the value before the switch (level_over_step).

Its switch times must be time points of the run (forced_time_points), so that a step does not straddle a change.

edge is the condition at r = a:

• :zero_flux (default, DREAM's BC_CONSERVATIVE, the setting of DREAM's disruption examples);
• :absorbing: T = 0 on the edge face, over the distance from the last centre to the face;
• :absorbing_dream: DREAM's BC_F_0, T = 0 at the distance of the last two centres;
• :constant_gradient_dream: DREAM's BC_DF_CONST, the gradient of the last two cells continued one cell outside (two radii at least).

heat_flux_through_edge(term::HeatTransportTerm, u, layout, t)

The heat that leaves through the edge of the plasma at the time t [s], −V' χ ∂_r T at r = a [W/m]: zero with edge = :zero_flux. With the spatial Jacobian V' = r of a cylinder, it is the heat per second through the edge per metre of the cylinder divided by 2π. The volume integral Σ V' Δr of the term, in the same unit, is minus this flux.

IonEnergyExchangeTerm(plasma; side, ion_energy=:ion_thermal_energy,
                      temperature=:cold_electron_temperature, density=:cold_electron_density,
                      ion_ion=true, ion_coulomb_logarithm=:nrl, masses=nothing, mass_table=:standard,
                      ions=nothing, initial_time=0.0)

The collisional exchange of thermal energy between Maxwellian species (DREAM's src/Equations/Fluid/MaxwellianCollisionalEnergyTransferTerm.cpp, the prefactor at line 49), in W/m³: the power that species a receives from species b,

\[Q_{ab} = \sqrt{3/π}\,\frac{(e^2/ε_0)^2 \sqrt{m_a m_b}}{4π} \ln Λ\, \sqrt{N_a N_b}\, (nZ^2)_a (nZ^2)_b\, \frac{N_a W_b − N_b W_a}{(m_b N_b W_a + m_a N_a W_b)^{3/2}},\]

the two-temperature relaxation of Spitzer (L. Spitzer, Physics of Fully Ionized Gases, 2nd ed., Interscience, 1962), Q_{ab} = (3/2) N_a ν_ε (T_b − T_a) with the energy relaxation rate of the NRL Plasma Formulary (A. S. Richardson, NRL/PU/6770–19-652, Naval Research Laboratory, 2019). For the electrons N = (nZ²) = n_e and W = 3/2 e n_e T_e from the unknowns density and temperature; for an ion species N is its total density, (nZ²) = Σ n_{Z₀} Z₀² (neutrals hold energy and exchange none), and W the component of the unknown ion_energy (J/m³, one component per species of plasma, in its order). The logarithm of electrons and ions is the thermal one, as DREAM's. Among the ions it is

• ion_coulomb_logarithm = :nrl (default): the mixed ion-ion logarithm of the NRL Plasma Formulary (A. S. Richardson, NRL/PU/6770–19-652, Naval Research Laboratory, 2019, p. 34), λ_ab = 23 − ln[Z_a Z_b (μ_a + μ_b)/(μ_a T_b + μ_b T_a) (Σ_γ n_γ Z_γ²/T_γ)^{1/2}] with the densities in cm⁻³, the temperatures in eV and μ = m/m_p. Formula (c) of the formulary sums over the two colliding ions; its definition of the screening length on the same page sums over every species whose thermal speed exceeds the relative speed of the pair, ū² < kT_γ/m_γ, which holds for the electrons unless the ions are hotter by the mass ratio, so the sum here takes every ion species and the electrons (n_e/T_e). Each ion species takes its own temperature T = W/(3/2 e N) and the electrons theirs, at least the 0.5 eV of DEFAULT_TEMPERATURE_FLOOR, and for Z its mean charge (nZ²)/(nZ), (nZ) = Σ n_{Z₀} Z₀. The formula holds for λ ≫ 1 (its corrections are O(1/λ)). With Ti = Te, deuterium at 10²⁰ m⁻³, argon at 10¹⁹ m⁻³ in the charge state that dominates in collisional ionization equilibrium in the limit of zero density (P. Bryans, E. Landi and D. W. Savin, Astrophys. J. 691, 1540 (2009), doi:10.1088/0004-637X/691/2/1540, Table 2) and quasi-neutral electrons, it is 8.7 for D with Ar³⁺ at 10 eV and 9.7 for D with Ar⁸⁺ at 54 eV. At the density of this plasma the state of 10 eV lies about one charge state higher: coronal_equilibrium with the ADAS rates at ne = 1.3 × 10²⁰ m⁻³ has Ar⁴⁺ dominant, and D with Ar⁴⁺ gives 8.3 (ne = 1.4 × 10²⁰ m⁻³); at 54 eV it keeps Ar⁸⁺. Ar¹⁶⁺ left at 10 eV by a fast cooling of the electrons, out of equilibrium, gives 5.9 with D (ne = 2.6 × 10²⁰ m⁻³). Ar¹⁶⁺ with Ar¹⁶⁺ gives 3.1 when the argon is held in two species, together at 10¹⁹ m⁻³: a species is not paired with itself, so a plasma of one argon species has the pair D–Ar alone. The formula falls below zero for two species of a highly charged ion, cold and dense: at these densities and with Ti = T_e, Ar¹⁶⁺ with Ar¹⁶⁺ crosses zero at 1.26 eV and gives −1.4 at the floor of 0.5 eV, where D with Ar¹⁶⁺ still gives 1.4. It is cut at 0 so that the exchange keeps its sign, a guard against cold iterates of the solver and cold prescribed input.

• ion_coulomb_logarithm = :dream: DREAM's 17.3 − ½ ln(n_free/10²⁰) + 3/2 ln(T_e/1 keV) with the electron temperature and the free-electron density, the same for every pair (src/Equations/CoulombLogarithm.cpp:107-113, with no source named there; BootstrapCurrent.cpp:181 calls the form Wesson's, which takes the ion temperature). With Te in place of Ti it is too small where the ions are hotter than the electrons: 10.4 at Te = 10 eV and 10²⁰ m⁻³ against 14.9 of the NRL form for two deuterium species at Ti = 1 keV (the exchange among them 1.4 times too slow after the thermal quench). Without the charges it is too large for a pair with a multiply charged ion: 10.3 at Ti = Te = 10 eV and n_free = 1.3 × 10²⁰ m⁻³ against 8.7 of the NRL form for deuterium with Ar³⁺ in the plasma above (the exchange about 1.2 times too fast).

• side = :ions, in the row of ion_energy: each component receives Q_{i,e} and, with ion_ion, Σ_j Q_{ij}.

• side = :electrons, in the row of the temperature (which holds the time derivative of the electron energy): −Σ_i Q_{i,e}.

A system needs both sides on the same unknowns, so the exchange conserves energy; a system with one side alone is refused. The masses are those of ion_mass with mass_table (:standard, the standard atomic weights and the nuclei of the hydrogen isotopes, or :dream, DREAM's rule). masses, when given, replaces them for all species: one mass in kg per species of plasma, in its order. A mass outside 1e-27 kg to 1e-24 kg is refused, which catches a mass given in atomic mass units or in grams. The densities of the ions come from the prescribed plasma, a PlasmaState or a function of time, or with ions from that radial unknown of the charge states, whose columns the term then assembles directly; a species of zero density is skipped. ion_ion = true adds the exchange among the ion species, as DREAM does whenever the ion energies are evolved (src/Settings/Equations/T_i.cpp:109). A species of zero total density exchanges nothing and gets no derivative there: the exchange goes as the square root of the density, whose derivative at zero is singular. A Newton iterate with a negative ion energy makes the denominator non-positive, where the exchange is zero with a zero derivative: declare the unknown non-negative, RadialUnknown(:ion_thermal_energy; components, nonnegative = true).

ion_temperatures(u, layout, plasma; ion_energy=:ion_thermal_energy, ions=nothing, t=0.0)

The temperature [eV] of every ion species at every radius (species × radii), T = W/(3/2 e N), from the state u and the densities of plasma at the time t [s]; NaN where a species has no density. With ions the densities N are those of that unknown of u, as in set_ion_thermal_energy!.

set_ion_thermal_energy!(u, layout, plasma, temperature; ion_energy=:ion_thermal_energy, ions=nothing, t=0.0)

Set the thermal energy of every ion species in the state u, W = 3/2 e N T with the total density N of the species in plasma at the time t [s] and the temperature temperature [eV]: a number, one value per species, or a matrix of species × radii (DREAM's initial W_i from T_i, T_i.cpp). With ions, the radial unknown of the charge states, N is that of the state u, as the energy exchange reads it: after deposit_neutrals! the prescribed densities no longer hold.

MAXIMUM_EQUILIBRIUM_ITERATIONS

The largest number of iterations of the free-electron density in fill_coronal_equilibrium!: 50, DREAM's MAX_ITER (src/Settings/IonEquilibrium.cpp:82, DREAM v26.5).

IonDiffusionTerm(composition; diffusion=nothing, diffusion_at=nothing, ions=:ion_densities,
                 states=:charged, edge=:zero_flux)

Radial diffusion of the densities of the ion charge states of composition, (1/V') ∂_r (V' D ∂_r n_j) [m⁻³ s⁻¹] in the rows of the ion densities ions (DREAM's IonChargedPrescribedDiffusionTerm). The diffusivity D [m²/s] is diffusion, a radial profile on the faces (radial_values), or diffusion_at, a function of time that returns one, such as a PiecewiseConstant transport history: its switch times must be time points of the run (forced_time_points). The flux through the face between two cells is V' D (n_k − n_{k−1})/Δr with the distance of the centres, as for the heat (HeatTransportTerm); nothing crosses the first face.

• states = :charged (default) diffuses the charge states Z0 ≥ 1, DREAM's charged diffusion; :all the neutrals too.
• edge = :zero_flux keeps the ions in the plasma; :absorbing sets the density to zero on the edge face, over the distance from the last centre to the face.

The term has the sign convention of every other row: a positive rate adds particles (DREAM's ion diffusion needs a negative coefficient for the sign of its ion transient term, src/Equations/Fluid/IonTransientTerm.cpp:78).

IonRateTerm(composition; ions=:ion_densities, temperature=:cold_electron_temperature,
            density=:cold_electron_density, kinetic_ionization=nothing,
            extrapolation=:dream, directory=adas_cache_directory())

The ionization and recombination of the charge states of composition (IonComposition) by the cold electrons, in the rows of the ion densities ions (DREAM's IonRateEquation):

\[\frac{dn_j}{dt} = n_e \left(I_{j-1} n_{j-1} - I_j n_j + R_{j+1} n_{j+1} - R_j n_j\right)\]

for each species and charge state j = 0, …, Z, with the ionization coefficients I (ADF11 SCD) and the recombination coefficients R (ACD) of OpenADAS (H. P. Summers, The ADAS User Manual, version 2.6 (2004), http://www.adas.ac.uk, with the data of OpenADAS, https://open.adas.ac.uk) at the electron density n_e and temperature of the radial unknowns density and temperature (adas_coefficient with extrapolation), as DREAM's src/Equations/Fluid/IonRateEquation.cpp. The default extrapolation = :dream takes DREAM's extrapolation below the tables (a choice of this package; see the docs page Choices of this package). The species are named by their element, whose tables are read from directory when the term is built (read_adas).

The sum over the charge states of a species is zero, so the term moves particles between the states and conserves each species. Its Jacobian comes from the colored AD of its residual, with the derivatives of the coefficients through ForwardDiff in place of DREAM's one-sided differences, and without DREAM's floors on the density (1 m⁻³) and the coefficients (1e-40) in the columns of the temperature and the density.

kinetic_ionization, a KineticIonizationTerm of the same system, splits the ionization at the threshold of its region: the electrons of f above the threshold ionize through the cross sections of that term, and the cold electrons below it through the fluid coefficient of fluid_ionization_coefficient, the ADAS rate less the part of the Maxwellian above the threshold. No electron is counted twice. The threshold is a multiple of the thermal momentum of at least the factor of MINIMUM_SPLIT_THRESHOLD_FACTORS for the form of its cross sections. Without kinetic_ionization the fluid rate is the whole ADAS rate, which next to a kinetic ionization term is DREAM's choice (ions.cpp:280): the Maxwellian tail above the threshold ionizes twice.

WeightedSumTerm(unknown, weights)

Σ_k w_k x_k over the components x_k of the radial unknown unknown at each radius, added to the row of a radial unknown with one component: the free-electron density n_free = Σ Z0 n_ij of the ion densities with the weights of free_electron_weights. The term is linear, and its Jacobian entries are assembled directly, since its row holds every component of a radius.

coronal_equilibrium(term::IonRateTerm, s, free_electron_density, temperature;
                    species_density=1)

The densities of the charge states of species s in coronal equilibrium at the free-electron density free_electron_density [m⁻³] and the temperature temperature [eV], with the total density of the species species_density: the steady state of the rates of IonRateTerm, where the flow between neighbouring states vanishes, n_{j+1}/n_j = I_j/R_{j+1}. The ratios are multiplied in log space, so states far from the dominant one underflow to 0 and the large ratios at the extremes of the temperature do not overflow (DREAM multiplies the ratios directly, src/Settings/IonEquilibrium.cpp:100-137).

deposit_neutrals!(u, layout, composition, deposits...; ions=:ion_densities)

Add neutral atoms to the unknown ions of the state u at the start of a run (DREAM's IONS_DYNAMIC_NEUTRAL): each deposit species => density adds density [m⁻³], a number or one value per radius, to the neutral state of the species of that name. The rate equations ionize them from the first step on. DREAM's reference case deposits neutral deuterium and argon at 1 eV, uniform in radius (examples/Disruptions/generate.py); a deposit of deuterium into the species of the background deuterium gives the same densities as DREAM's separate species, since the rates are linear in the densities.

fill_coronal_equilibrium!(u, layout, term::IonRateTerm, temperature;
                          species=term.composition.names)
fill_coronal_equilibrium!(u, layout, composition, temperature; species=composition.names,
                          ions=:ion_densities, extrapolation=:dream,
                          directory=adas_cache_directory())

Put the species species (a name, as a String or a Symbol, or a collection of names) of the composition of term, an IonRateTerm, in coronal equilibrium in the state u, in the radial unknown of the ion densities of term, at the temperature temperature [eV], a profile in any form of radial_values on the cell centres. The total density of each of these species is that of u and is kept; the other species keep their charge states. The rates are those of term with its extrapolation; its kinetic_ionization is not used, as DREAM's init_equil takes the fluid rates alone. The second method builds the term from composition (IonComposition), ions, extrapolation and directory, which reads the ADF11 tables of every species. Returns u.

The equilibrium depends on the free-electron density, which it changes, so the free-electron density is iterated to a fixed point at each radius as DREAM's init_equil does (EvaluateIonEquilibrium, src/Settings/IonEquilibrium.cpp:63-155, DREAM v26.5): it starts at the free electrons of the other species plus the total density of each species in equilibrium, the charge states of coronal_equilibrium follow at that density, and the free-electron density is summed from them, until it changes by at most eps() of itself or after MAXIMUM_EQUILIBRIUM_ITERATIONS iterations. A radius where that density starts at or below zero, where the species in equilibrium have no atoms and the other species no free electrons, is left as it is: no rates are defined there.

The free-electron and cold-electron densities of the state are not changed: the caller sets them from the new charge states, as the rows of the system give them (population_equations, fill_initial_values!) or by hand.

Throws an ArgumentError for a species that is not in the composition or is named twice, for a radial unknown of the ion densities whose components are not those of the composition, for a total density of a species in equilibrium that is negative or not finite, or for 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.

free_electron_weights(composition::IonComposition)

The charge Z0 of every component of the ion densities: the weights of WeightedSumTerm for the free-electron density Σ Z0 n_ij.

ion_densities(u, layout, composition; ions=:ion_densities)

The IonSpecies of composition with the charge-state densities of the unknown ions in the state u.

ion_densities_unknown(composition; ions=:ion_densities, scale=1e19, limit_tolerance=1e-10)

The radial unknown ions of the densities [m⁻³] of the charge states of composition (IonComposition), with one group of components per species (component_groups). Its lower limit keeps every charge state at or above −limit_tolerance times the density of its species (LowerLimit(reference = :group)): a state that starts at 0 is not refused for a round-off below 0, and the limit scales with each species, so a trace impurity has its own. DREAM gives the charge states no limit of their own: it damps the Newton update to keep Σ Z0² n per radius above a tenth of its value (src/Solver/SolverNonLinear.cpp:489-508), and the total density N_i of each species likewise where N_i is an unknown (src/Solver/SolverNonLinear.cpp:455-456, 476-487). The default Newton tolerance of such an unknown is relative to the density of each species as well (tolerance).

ion_rate_matrix(term::IonRateTerm, s, free_electron_density, temperature)

The matrix A of the rates of species s at the free-electron density free_electron_density [m⁻³] and the temperature temperature [eV], dn/dt = A n for the densities n of its charge states (IonRateTerm). Its columns sum to zero.

population_equations(composition; ions=:ion_densities, density=:cold_electron_density,
                     hot_density=nothing, runaway_density=nothing)

The closure of the cold-electron density, 0 = Σ Z0 n_ij − n_hot − n_re − n_cold: the free electrons of the ion densities ions (WeightedSumTerm) less those of the distributions, given as the radial unknowns hot_density and runaway_density of their densities (DREAM's n_cold), in the row of the cold density density. Without them the cold density is the free-electron density. The cold density is the target of the e-e collisions, read by a PlasmaFromUnknowns with the same density.

set_ion_densities!(u, layout, species; ions=:ion_densities)

Set the unknown ions of the state u to the charge-state densities of species, a list of IonSpecies or a PlasmaState, whose species and order define the composition of the unknown.

MINIMUM_SPLIT_THRESHOLD_FACTORS

The smallest factor of the thermal momentum at which the threshold of a KineticIonizationTerm may split the ionization with a fluid term (the kinetic_ionization keyword of IonRateTerm and RadiatedPowerTerm), per form of its cross sections (a choice of this package; see the docs page Choices of this package): 7 for CHIANTI's (:chianti, the default), and 5, the hot-region threshold of DREAM's disruption example (examples/Disruptions/generate.py:234), for DREAM's (:dream) and for the fits of DREAM's method (:bethe).

Runes2P's default hot region takes 7 for every form (ThermalMomentum).

KineticIonizationTerm(composition, kinetic, grid; ions=:ion_densities, region=nothing,
                      ionization_cross_section_form=:chianti)

The ionization of the charge states of composition (IonComposition) by the electrons of the kinetic unknown kinetic on grid, in the rows of the ion densities ions (DREAM's IonKineticIonizationTerm):

\[\frac{dn_j}{dt} = K_{j-1} n_{j-1} - K_j n_j, \qquad K_j = \frac{1}{V'} \sum_{\text{cells}} e(p)\, v σ_j(p)\, 𝒱' Δp Δξ\, f,\]

with the ionization rate K_j [s⁻¹] of the charge state j, the cross sections ionization_cross_section of the form ionization_cross_section_form (:chianti, the default, CHIANTI's parameters continued relativistically, with DREAM's form for the elements CHIANTI's table lacks; :dream, DREAM's form and fits; or :bethe, DREAM's form with its relativistic term divided by β², refitted), tabulated at the cell centres of the momentum grid of each radius (DREAM tabulates the grid of radius 0 alone), and the envelope e(p) of the MomentumRegion region: the hot electrons above a threshold that may follow the temperature unknown, or the whole grid without a region. Each species is conserved.

Paired with an IonRateTerm through its kinetic_ionization keyword, the electrons below the threshold ionize at the fluid coefficient of fluid_ionization_coefficient, so a Maxwellian f ionizes at the ADAS rate for every threshold. The term has no energy of its own: the friction on bound electrons of the Fokker-Planck term already takes energy from f for the inelastic collisions. With bound_friction = :ionization_split of CollisionSettings and its ionization_threshold at the factor of this term's region, that friction holds the energy I_j + ⟨E_s⟩_j of each ionization the term counts, up to the discretization. Without a region the two agree to the discretization; at a threshold the friction takes a smooth step of its own on the faces (ionization_threshold_width) while this term takes the envelope of its region in the cells, so the two differ by the electrons within that width of the threshold.

The Jacobian is assembled directly: the columns of f (every cell of a radius, the block DREAM's reference case drops), of the charge states, and of the temperature through the threshold of the region.

SecondaryHeatingTerm(kinetic_ionization::KineticIonizationTerm;
                     widths=SECONDARY_SPECTRUM_WIDTHS,
                     temperature=:cold_electron_temperature, cut=:none, spectrum=:beb)

The heating [W/m³] of the cold electrons by the secondary electrons of the kinetic ionization kinetic_ionization, for the row of their temperature:

\[Q = e \sum_{s,j} n_{s,j} \frac{1}{V'} \sum_{\text{cells}} e(p)\, v σ_j(p)\, ⟨E_s⟩_j(E)\, 𝒱' Δp Δξ\, f,\]

with the rates of kinetic_ionization, the envelope e(p) of its region (1 without one) and the mean energy ⟨E_s⟩_j of the secondaries of the spectrum spectrum and the cut cut below. A kinetically ionized electron joins the cold electrons through the closure of the cold density; without this term it is born at zero energy.

• cut = :none (the default): the mean over every secondary, up to (E − I)/2. All secondaries join the cold electrons, also those born above the hot threshold, so this is the moment of the spectrum over the secondaries that the cold electrons receive, and the one that the friction with bound_friction = :ionization_split takes from f (CollisionSettings). A KnockOnTerm on the bound electrons also puts the secondaries above its cutoff into f, so next to it their energy is counted twice. Their share of the mean is large: with the spectrum of Opal et al. at 1.3 MeV, the secondaries above 1 keV carry about 59 % of it.
• cut = :threshold: the mean below the hot threshold E_c, conditional on lying there, for every secondary.

spectrum = :beb (the default, a choice of this package; see the docs page Choices of this package) takes the binary-encounter-Bethe spectrum of Kim and Rudd with the ionization energy of each charge state as its binding energy (beb_mean_secondary_energy), which needs no width; spectrum = :opal the spectrum of Opal et al. with the widths widths [eV], the neutral element's for every charge state. The friction of the Fokker-Planck term on bound electrons takes the energy of the ionization from f already, so in the total energy of the electrons the term moves the energy of the secondaries from that sink to the cold electrons. That holds where the friction takes the energy. With the default screening window of the collisions (screening_window = :excitation_energy) the friction on bound electrons is zero from the ionization potential up to the mean excitation energy of the state (about 190 eV for Ar⁰, more for the higher states): there f ionizes without paying, and this term adds the energy of the secondaries on top. screening_window = :dream narrows that band, and bound_friction = :ionization_split of CollisionSettings closes it: f then pays the ionization potential and the energy of the secondaries for every ionization above the threshold, and this term returns the energy of the secondaries to the cold electrons. DREAM has no such term: kinetic_ionization alone is DREAM's model. The region of kinetic_ionization sets E_c; its threshold follows the temperature unknown temperature or is fixed.

fluid_ionization_coefficient(adas, fraction, coefficient)

The ionization coefficient [m³/s] of the cold electrons below the threshold of a kinetic ionization: (I − A_>)/(1 − f_>), the ADAS coefficient I (adas, in m³/s) of the whole Maxwellian less its part above the threshold, A_> (coefficient, in m³/s), per electron below it, with the fraction f_> of the Maxwellian above the threshold (threshold_split). With the closure n_cold = n_M (1 − f_>) and the kinetic rate n_M A_> of a Maxwellian f of density n_M, the two add up to n_M I for every threshold the grid resolves.

The coefficient has no floor on its value. Its denominator is bounded by the grid: threshold_split takes 1 − f_> as the weight of the first momentum cell where the threshold lies in that cell, and the sum then holds only to within that part. A_> comes from the cross sections of the kinetic ionization and I from ADAS. The Maxwellian rate coefficient of DREAM's fits exceeds the ADAS coefficient by up to a factor of about 1.6 for a state of argon, that of CHIANTI's cross sections by up to 2.6 for a state of neon or argon (Ne²⁺ at 3 to 20 eV). Where the threshold lies below the ionization energy of a state, almost all of its Maxwellian ionization comes from above the threshold, so A_> > I and the coefficient is negative. On the grid of the tests (n_e = 5e19 m⁻³, 5 to 50 eV), with DREAM's cross sections, it is negative for the dominant argon states at threshold factors 2 and 3, by 3 to 18 % of I, and for no state with a coronal fraction above 1e-17 at factors 5 and 7. From 1 eV to 10 keV, for D, Ne and Ar, it is negative for no state with a coronal fraction above 1e-3 with DREAM's cross sections at factors 5 to 7. With CHIANTI's it is negative at factor 5 for three such states of a cold plasma: Ar⁺ at 1 eV (−0.19 I, coronal fraction 0.087), Ar²⁺ at 1.6 eV (−0.50 I, 1.2e-3) and Ne²⁺ at 2.5 eV (−2.8 I, 0.025), and for none at factors 6 and 7. With a Maxwellian f the sum with the kinetic rate is n_M I whatever the sign. The split refuses thresholds below the factor of MINIMUM_SPLIT_THRESHOLD_FACTORS of its form.

ionization_energy_flows(kinetic_ionization::KineticIonizationTerm, u, layout, t;
                        widths=SECONDARY_SPECTRUM_WIDTHS, spectrum=:beb)

The energy [W/m³ per radius] that the kinetic ionization kinetic_ionization takes from f at the state u and the time t [s], in two parts:

• potential: the ionization potentials, e Σ I_j K_j n_j, stored in the ions and released at recombination (the gain of RadiatedPowerTerm(...; radiation = :recombination));
• secondaries: the kinetic energy of the secondary electrons, which SecondaryHeatingTerm with its default cut = :none gives to the cold electrons; with spectrum = :opal and the widths widths [eV], nothing when an element of the composition has no width of the secondary spectrum (SECONDARY_SPECTRUM_WIDTHS), since potential needs none.

Both are part of the sink of f through the friction on bound electrons (KineticEnergyLossTerm with part = :sinks), so an energy diagnostic that books that sink as lost and the recombination radiation as radiated counts potential twice: it is the share of the sink that goes to the potential energy of the ions.

kinetic_ionization_rates(term::KineticIonizationTerm, u, layout, t, ir)

The rate coefficients K [1/s] of all charge states at radius ir, the moments of f of the state u at the time t [s] with v σ above the threshold of the region, one per component of the ion densities.

threshold_split(term::KineticIonizationTerm, ir, u, layout, t)

The parts of a Maxwellian at the cold temperature of radius ir above the threshold of the region of term, on its momentum grid, for fluid_ionization_coefficient: the fraction f_> of its density and the rate coefficient A_> [m³/s] of each charge state, Σ e(p) v σ(p) M w / Σ M w with the envelope e, the cell weights w of a density moment and the relativistic Maxwellian M = exp(−(γ − 1)/Θ) at the cell centres. The sums run over the cells of the grid, the quadrature of the term itself. The temperature is that of the region's unknown in the state u, in its number type.

A threshold p₀ in the first momentum cell, below its upper face, is not resolved by the grid: the part of the Maxwellian below it then comes from the envelope within one cell, not from the shape of the Maxwellian, and tends to 0 as the threshold nears the lower face, while the fluid coefficient divides by it. There 1 − f_> is taken as the weight of that cell, M_1 w_1 / Σ M w (a choice of this package; see the docs page Choices of this package): at the threshold factors the split takes, 5 thermal momenta or more, a threshold in the first cell leaves nearly the whole Maxwellian in that cell, below the threshold. A warning (once per session) names the radius and the temperature; a finer first momentum cell, or a threshold of more thermal momenta, resolves it.