Comparison with DREAM

This page compares Runes2P with DREAM, the code whose design it follows (M. Hoppe, O. Embréus and T. Fülöp, Comput. Phys. Commun. 268, 108098 (2021), doi:10.1016/j.cpc.2021.108098). The comparison is with DREAM release v26.5 and its branch boltzmann-eqterm at commit 93d27e1, the builds of the Reference data. The file and line citations of DREAM's source refer to commit 93d27e1 unless they name another version.

Features of DREAM

The tables follow the groups of DREAM's settings (py/DREAM/Settings). The column "DREAM" names the setting, class or unknown. The column "Runes2P" gives the state, the names that implement the feature and, for a partial feature, the part that is in the code. A cell that ends with "(note)" has a note under its table, which gives the detail.

Runes2P is cylindrical: a straight magnetic field, constant on each flux surface, so no trapping and no flux-surface averages (cylindrical_grid). Every toroidal feature of DREAM is therefore not implemented.

Run modes

Runes2P has no switch of the run mode. A run is an EquationSystem of the terms the caller chooses, so a mode is a choice of unknowns and terms.

• Superthermal. ConductivityTerm is for Ohm's law of the cold electrons. FrictionEnergyTransferTerm is for their heating by $f_\mathrm{hot}$. The disruption case in this mode is compared with DREAM end to end (The disruption case).

Grids

• Major radius. A cylinder has none, as in DREAM's cylindrical grid. major_radius is only an input of the Rechester–Rosenbluth coefficient (add_rechester_rosenbluth_diffusion!, RechesterRosenbluthUpdate, HeatTransportTerm).
• Inner radius. A RadialGrid that starts above the axis has no flux through its first face.

Collisions

Electric field and current

Temperature and energy

• Collisional heating by $f_\mathrm{hot}$ and $f_\mathrm{re}$. For the hot tail, Runes2P takes DREAM's moment. For the runaway grid, Runes2P's default is free_electron_energy_transfer_term, the energy change of the electron-electron operator (KineticEnergyLossTerm). DREAM's moment is the option for the runaway grid (Defaults that differ). The moment is compared with DREAM's saved heating of both grids (test_collisional_heating_dream.jl). It is also compared in a run of the runaway grid with the cold temperature self-consistent against DREAM's run (test_runaway_heating_dream.jl).

Ions and ionization

• Start in coronal equilibrium. fill_coronal_equilibrium! is DREAM's fixed point of the free-electron density over the species started in equilibrium. It takes the charge states of coronal_equilibrium. It is compared with DREAM's examples/ionequil (test_initial_state_dream.jl).

Runaway generation

• Fluid runaway density. Above a hot-tail grid without a runaway grid, it is a radial unknown fed by ExternalAvalancheTerm and the flux through $p_\mathrm{max}$ (MomentumBoundaryFluxTerm). With a runaway grid, it is fed by the particle flux through the interface (InterfaceParticleFluxTerm).
• Fluid avalanche. avalanche_growth_rate is the Rosenbluth–Putvinski rate. The fluid rate enters $n_\mathrm{re}$ only above a hot-tail grid, through ExternalAvalancheTerm.

Transport

• Prescribed transport coefficients of $f$. The caller sets them on the faces of the grid. There are no coefficients given in time, radius and momentum and interpolated to the grid.

Kinetic options

• Particle source. The explicit source has no packaged form. The disruption case builds it from an unknown of the source with its own row (test/example_disruption.jl).
• Initial distribution. A run takes any state array of the caller on the grid of the run. DREAM's Maxwellian is built from profiles of $n$ and $T$.
• Initial runaway distributions. fill_runaway_distribution! gives the four shapes in the first momentum cell. They are compared with DREAM's (test_initial_state_dream.jl).

Solver and time stepping

• Residual check. The residual is a convergence criterion except for the linear steps that linear_acceptance of NewtonSolver accepts on the update. The steps at the round-off floor of the residual are the other exception.

Output

The disruption case

The page The disruption case compares the two codes in whole runs of DREAM's disruption case of examples/Disruptions/generate.py. It covers the reduced size of the tests and the whole quenches of scenarios 1 and 4 at full size.

What Runes2P has that DREAM does not

• Exponential fitting, the default scheme of fokker_planck_operator (Positivity and flux limiters).
• The flux limiter van Albada, scheme = :van_albada of AdvectionDiffusion.
• A smooth momentum grid, tanh_p_faces.
• A pitch grid per radius, with the same cell counts and momentum faces at all radii, and conservative radial fluxes between different pitch grids (PitchOverlap, radial_face_overlap).
• A conservative flux between the hot-tail and the runaway grid by overlap cells, pitch_mapping = :overlap of MomentumInterfaceTerm.
• The momentum-conserving field-particle term of the electron-electron collisions, FieldParticleTerm.
• The conductivity of the discrete kinetic operator, kinetic_conductivity, as the correction of the ohmic current (kinetic = :discrete_operator of ConductivityModel).
• The total energy of the electrons in the fully kinetic mode, total_energy_equation, which closes the energy of the cold electrons.
• A smooth window of partial screening above the excitation energy, screening_window of CollisionSettings.
• The friction of the ionization on $f$, bound_friction = :ionization_split of CollisionSettings.
• The heating of the cold electrons by the secondaries of the kinetic ionization, SecondaryHeatingTerm.
• The split of the ionization at the threshold of the kinetic ionization, with the fluid coefficient of the cold electrons below it (fluid_ionization_coefficient, kinetic_ionization of IonRateTerm).
• CHIANTI's cross sections of the kinetic ionization continued by the relativistic Bethe form, and DREAM's form with the factor $1/β^2$ (form = :chianti, :bethe of ionization_cross_section).
• The knock-on pitch integrated over each pair of momentum cells and averaged over the primary pitch cell (MollerPitchKernel).
• An implicit knock-on term, treatment = :implicit of KnockOnTerm.
• A knock-on cutoff per radius that follows the plasma during a run, KnockOnCutoffUpdate.
• The extrapolated flux at $p_\mathrm{max}$ by DREAM's stated definition, p_upper_boundary = :extrapolated_flux_consistent.
• Absorbing edges with $f = 0$ or $T = 0$ on the edge face: r_boundary = :absorbing of AdvectionDiffusion, edge = :absorbing of HeatTransportTerm.
• A safety factor and a major radius in the Rechester–Rosenbluth coefficient (add_rechester_rosenbluth_diffusion!).
• Lower limits per unknown, relative to a reference unknown or to the density of a species (LowerLimit), and the temperature floor of PlasmaFromUnknowns.
• Tolerances of the charge states relative to the density of their species, group of Tolerance.
• A step controller with an error estimate from a predictor, and variable-step BDF2 (StepController).
• A backtracking line search of the Newton iteration, line_search of NewtonSolver.
• Preconditioners of GMRES for systems coupled through the fluid unknowns: FluidSchurLU, RadiusBlockLU, RadiusBlockILU.
• Exact Jacobians by automatic differentiation, the columns of the plasma unknowns included (jacobian!, jacobian_correction_product!).
• A schedule of legs with prescribed steps or a controller (Leg, time_points).
• Checkpoints and segments of a run (CheckpointSettings, read_output). With SparseLU a resumed run repeats the uninterrupted one bit for bit.
• The system as an ODE function of SciML (ode_function).

Departures from DREAM

The sections below list the defaults that differ from DREAM's and the behaviour without an option for DREAM's. The last three give the details of the extrapolated flux at the highest momentum, the approximate Jacobians and the line search.

Defaults that differ

Each row names the default of Runes2P and the option that gives DREAM's behaviour. The reasons of the defaults that are choices of this package are on Choices of this package.

The sources of the values of Runes2P in the table are these.

Ion sizes of partial screening. DREAM's lookup has a list of nuclear charges with 17 argon entries for 18 charge states. This shifts the sizes from Ar¹⁷⁺ on.

Flux through $p = 0$ on a runaway grid from $p = 0$. Runes2P takes the operator's default, p_zero_boundary = :outflow of fokker_planck_operator, as DREAM's hot-tail grid. The friction carries electrons out by $\lim_{p \to 0} p^3 ν_s$ (p3_slowing_down_frequency_at_zero). With it the runaway grid of the case heating of reference/dream/cases/initial_state.py loses 1.55e-4 of the density of $f_\mathrm{re}$ in the first step, where DREAM's stays constant to round-off.

The caller sets the option p_zero_boundary = :zero_flux for a runaway grid from $p = 0$, as both compat values of runaway_heating_example do (test_runaway_heating_dream.jl). In DREAM's source:

• $f_\mathrm{re}$ has no internal boundary condition (src/Settings/Equations/f_re.cpp:87).
• $f_\mathrm{re}$ keeps the Dirichlet lower boundary of its advection in $p$ (fvm/Equation/AdvectionTerm.cpp:36).
• Unlike $f_\mathrm{re}$, $f_\mathrm{hot}$ has an internal boundary condition (src/Settings/Equations/f_hot.cpp:96). With it, the lower boundary of its advection in $p$ is mirrored (AD_BC_MIRRORED) in place of Dirichlet (src/Settings/Equations/distribution.cpp:234-237).

Heating of the cold electrons by a runaway grid. Runes2P's default is the energy change of the electron-electron operator of the grid, with a linear scheme of the operator. DREAM's moment, FrictionEnergyTransferTerm, also counts the energy that the friction carries out of the grid through $p = 0$ where that flux is applied (p_zero_boundary = :outflow). The test case runaway_heating_example (test/example_cases.jl) takes the option at compat = dream. It is compared with DREAM's run in test_runaway_heating_dream.jl.

Relativistic corrections of Gould's bremsstrahlung. The coefficients of $T/m_ec^2$ are these.

Behaviour without an option for DREAM's

In these points Runes2P differs from DREAM, and no option gives DREAM's behaviour.

• Newton's method converges on both the update and the scaled residual, except for the linear steps that linear_acceptance accepts on the update (NewtonSolver). A stiff step at the round-off floor of its residual is accepted on the update too, within the round-off bound of its rows or where the residual stopped decreasing. DREAM converges on the update. With setCheckResidual on, its default, DREAM warns when the residual is not converged (src/Solver/SolverNonLinear.cpp:305-310, py/DREAM/Settings/Solver.py:28).
• A failed Newton step is reported as a status of NewtonResult. The time stepper retries it with a shorter step.
• Each charge state has a lower limit relative to the density of its species (ion_densities_unknown). In DREAM the charge states have no limit of their own (src/Solver/SolverNonLinear.cpp:441-508). DREAM's MaximalPhysicalStepLength damps these quantities (src/Solver/SolverNonLinear.cpp:441-508):
  • $T_\mathrm{cold}$, $n_\mathrm{tot}$ and $n_\mathrm{cold}$.
  • $W_\mathrm{cold}$, $W_i$ and the species densities $N_i$, where they are unknowns.
  • $\sum Z_0^2 n$ per radius.
• A lower limit shortens only steps towards the limit, never growth (LowerLimit).
• A doubled step can fail on an unknown at zero with no absolute tolerance while every other unknown has no error. Such a step is retried with half its length (StepDoubling). DREAM's next step is unbounded there (src/TimeStepper/TimeStepperAdaptive.cpp:458-482).
• The ion rows take the sign convention of every other row, so IonDiffusionTerm takes a positive coefficient (DREAM: src/Equations/Fluid/IonTransientTerm.cpp:78).
• coronal_equilibrium multiplies the ratios of the rates in log space, where DREAM multiplies them directly (src/Settings/IonEquilibrium.cpp:100-137).
• The derivatives of the ADAS coefficients are exact, without DREAM's floors in the Jacobian (adas_coefficient, IonRateTerm).
• The Jacobian of ExternalAvalancheTerm holds every column of the rate, where DREAM's holds those of the field and of $n_\mathrm{tot}$.
• The temperature column of ParticleSourceTerm holds the derivative of the normalization of its shape, which DREAM's column omits (src/Equations/Kinetic/ParticleSourceTerm.cpp).
• The recombination gain of radiation = :recombination takes the ionization energy of the recombined state at every radius. DREAM's takes that of the recombining state from its second radial cell on (src/Equations/Fluid/RadiatedPowerTerm.cpp). test/test_adas_cache.jl checks the radiated power in coronal equilibrium, which detects this mismatch.
• RechesterRosenbluthUpdate takes the level of the perturbation at the start of each step. DREAM interpolates its table linearly in time and takes the value at the end of the step.
• StepDoubling reads the tolerances by unknown, for the acceptance and for the next step. DREAM's step update reads them by position in its list of non-trivial unknowns (GetErrorScale, src/ConvergenceChecker.cpp).
• check_every of StepDoubling takes the steps between two checked steps at the current length, the documented intent of DREAM's checkevery. DREAM's code takes them from a stale start time, with other lengths and one step more (src/TimeStepper/TimeStepperAdaptive.cpp).
• The :smooth envelope of a MomentumRegion whose limit lies at or above the last cell centre takes the width of the last cell. DREAM takes the limit minus the last centre (FindThresholdStep, fvm/Equation/MomentQuantity.cpp). That gives zero width and a NaN at that centre, and a jump of the envelope as the limit crosses it.

Extrapolated flux at the highest momentum

DREAM's BC_DPHI_CONST is defined by a constant $\mathrm{d}Φ/\mathrm{d}p$ at $p_\mathrm{max}$, with $Φ = \mathcal{V}' Γ^p$, the phase-space Jacobian $\mathcal{V}'$ times the momentum flux $Γ^p$. Its code (fvm/Equation/BoundaryConditions/PXiExternalLoss.cpp) shares the $\mathcal{V}'$ of the last interior face between the two terms of the extrapolation. It multiplies both terms by the advection coefficient $A^p$ of that face. The interpolation and the diffusion of the second term are those of the face before it. p_upper_boundary = :extrapolated_flux of AdvectionDiffusion ports this code literally, and the comparisons with DREAM use it. :extrapolated_flux_consistent extrapolates with the $\mathcal{V}'$ and the $Γ^p$ of each face, as the definition states. test/test_pmax_boundary.jl checks both. test/test_pmax_boundary_dream.jl compares the literal port with DREAM.

Approximate Jacobians

The residual of a FokkerPlanckTerm is always that of its scheme. jacobian selects the stored block alone:

• :exact differentiates the limiter.
• :frozen_limiter is DREAM's AD_INTERP_JACOBIAN_LINEAR.
• A linear scheme name gives that scheme's Jacobian (with_scheme).

jacobian = :upwind with a central residual is DREAM's AD_INTERP_JACOBIAN_UPWIND, which makes some of DREAM's runs at high momentum resolution converge (DREAM issue #126). A term with an approximate block declares jacobian_kind(term) === :approximate. Newton's method on the stored matrix is then a quasi-Newton iteration. With GMRES(operator = :exact) the stored block only preconditions the exact Newton operator (jacobian_correction_product!). test/test_jacobian_strategies.jl checks the modes.

Line search

DREAM's nonlinear solver (src/Solver/SolverNonLinear.cpp) damps the Newton step. It has no line search. With line_search = true, NewtonSolver backtracks on $\tfrac{1}{2}\lVert F \rVert^2$ of the scaled residual $F$ with the Armijo condition. It starts from the step length that the lower limits allow. The default is off, DREAM's behaviour. test/test_time_stepping.jl checks it.