Getting started

This page covers the installation, a first run with the examples, the units of the inputs, the settings of a run and its output files.

Installation

Runes2P needs Julia 1.13. It is not registered in the General registry and is added from its repository:

using Pkg
Pkg.add(url = "https://github.com/osrep/Runes2P.jl")

The HDF5 output is a package extension. It loads with using HDF5 next to using Runes2P and defines write_output and read_output. Checkpoints need it loaded.

The examples run from a checkout of the repository. Their environment, examples/Project.toml, takes Runes2P from the parent directory and adds HDF5:

git clone https://github.com/osrep/Runes2P.jl
cd Runes2P.jl
julia --project=examples -e 'using Pkg; Pkg.instantiate()'   # once
julia --project=examples examples/relaxation.jl
julia --project=examples examples/cooling.jl scheme=van_albada solver=gmres
julia --project=examples examples/avalanche.jl treatment=implicit solver=gmres

The scripts are described in The examples, and the arguments such as scheme= in Switches.

The examples

Each script solves one kinetic equation for the electron distribution $f(p, \xi, r)$ and runs without arguments. Its header lists its default values, the sources of these values and further command lines.

relaxation.jl

• The initial state is a Maxwell-Jüttner distribution at the plasma temperature, or at T0=.
• The plasma is deuterium, and the plasma and the field stay fixed.
• The operator is that of fokker_planck_operator.
• p_star= adds the fraction of the density above that momentum at every save.
• The plasma values are choices of this package (Choices of this package).
• By default the field lies between the critical and the Dreicer field (Choices of this package).

cooling.jl

• The temperature decays exponentially from the initial temperature T= to T_final= over t_end=.
• prescribed_fokker_planck_term rebuilds the collision coefficients from the temperature.
• The field stays at e_over_ec= times the critical field at the initial temperature.
• With nr= above one, radial diffusion along stochastic field lines couples the radii (add_rechester_rosenbluth_diffusion!). Its strength is dbb= ($\delta B/B$), and $f = 0$ at the wall. The diffusion is that of A. B. Rechester and M. N. Rosenbluth, Phys. Rev. Lett. 40, 38 (1978), doi:10.1103/PhysRevLett.40.38.
• The plasma is that of reference/dream/cases/cooling_positivity.py, cooled over 2 ms where that case takes 1 ms.
• The radial values are those of DREAM's SPI example.
• The grid and the time step are a choice of this package.
• With the default scheme, check the loss through p_max against a finer np= (Exponential fitting).

avalanche.jl

• The seed is a Maxwell-Jüttner distribution of temperature 100 keV.
• The knock-on collisions are those of KnockOnOperator and KnockOnTerm.
• The setup is that of examples/avalanche/generateCylindrical.py of DREAM's boltzmann-eqterm branch.
• It takes that example's values, among them superthermal collisions, the energy-dependent Coulomb logarithm and DREAM's ion-size table.
• The interpolation scheme is the default of this package, where the example uses TCDF (scheme=tcdf).

The switch treatment= of KnockOnTerm sets the time level of the knock-on term.

implicit applies the term matrix-free in each Krylov iteration of GMRES. The direct solvers do not take it, so with treatment=implicit the solver defaults to GMRES, and solver=gmres in the command above only repeats that default.

What every script prints

A run is one segment, or a chain of segments when it is split over several allocations (Checkpoints, segments and the stop). Every script prints the segment and its status, the solver statistics (SolverStatistics) and the negative entries of the saved states. It then prints the density per radius at every save. In a GMRES run the statistics include the Krylov iterations and factorizations. Before the timed run a script integrates over the first three time points, so that the timing leaves out the compilation. relaxation.jl and cooling.jl end with exit_status as the exit status of the process.

The parts of a run

The three scripts build a run in the same order.

1. A MomentumGrid from momentum faces (uniform_p_faces, biuniform_p_faces, tanh_p_faces) and pitch faces (uniform_xi_faces, biuniform_theta_xi_faces), combined with a RadialGrid by cylindrical_grid (Grids).
2. A PlasmaState with its IonSpecies, and the CollisionSettings of the collision frequencies (Collision models).
3. The kinetic operator, an AdvectionDiffusion with its interpolation schemes and boundaries (Positivity and flux limiters).
4. An UnknownLayout with one KineticUnknown :f, and an EquationSystem of one Equation: a TimeDerivative and the terms, among them FokkerPlanckTerm or prescribed_fokker_planck_term (Unknowns and the equation system).
5. A NewtonSolver with its linear solver, and integrate over the time points, which returns a TimeSolution (Time stepping and solvers).

The shared helpers of the scripts are in examples/common.jl: the switches, a Maxwell-Jüttner initial state normalized with density, the report and the output.

Units

The inputs of a run are in SI units, except temperatures and particle energies in eV and momenta in units of $m_e c$, with $m_e$ the electron mass and $c$ the speed of light. One setting is in minutes, guard_minutes of ProgressLog. The code checks almost no unit, so a value in another unit usually gives a wrong result without an error. The table lists the units by group, for the inputs and for the results of the functions it names.

The units of the unknowns that the package names by default are on the page Unknowns and the equation system. The output files store the unit of every unknown, of the times and of the grid faces (Unknowns). The settings record stores no unit with its values (Records).

Notes on the units

These notes explain the distribution, the ion masses, the thermal energy densities, the radial totals and the atomic data that the package reads.

The distribution $f$ is a density per unit of momentum volume, with the momentum $p$ in units of $m_e c$. So $f$ is in m⁻³. Its density moment at a radius is $n = \int 2\pi p^2 f \, \mathrm{d}p \, \mathrm{d}\xi$ (density). The tests against DREAM read DREAM's f_hot and f_re without a conversion factor and compare them with $f$ value by value.

IonEnergyExchangeTerm takes the masses of the species from ion_mass. Its keyword masses replaces them all, with one mass per species in kg.

A thermal energy density is $W = \tfrac32 e\, n T$ in J/m³, with $e$ the elementary charge (The energy of the cold electrons). The factor $e$ converts the temperature $T$ from eV to J (THERMAL_ENERGY_FACTOR).

radial_particle_content gives $n V' \Delta r$ at each radius, in m⁻¹. Its terms are the density $n$ in m⁻³, the spatial Jacobian $V' = r$ of a cylinder in m and the radial cell width $\Delta r$ in m. Its sum over the radii is total_particles. The number of electrons of $f$ per metre of the cylinder is $2\pi$ times total_particles.

The package converts the atomic data that come in other units. read_adf11 converts the densities of an ADAS file from cm⁻³ to m⁻³ and its coefficients from cm³ to m³ (ADAS rate coefficients). The fits of CHIANTI (CHIANTI_IONIZATION_FITS) keep the units of their files. These are eV for the shell energies and 10⁻¹⁴ cm² eV² for the scaled cross sections (ChiantiDirectShell), and Ry for the excitation energies (ChiantiAutoionization). chianti_cross_section gives the cross sections computed from them in m².

Switches

A script takes name=value arguments, and every switch has a default. The switches of common.jl are shared, and a script's header names the defaults it changes.

Notes on the grids and schemes

These notes belong to the rows grid=, scheme=, jacobian= and pmax_boundary= of the table above. Momenta are in units of $m_e c$ (Units).

grid=biuniform_270 is a fully kinetic hot-tail grid and a choice of this package (Choices of this package). It ignores the other grid switches.

The schemes up to upwind_second_order are linear in $f$. The rest are flux limiters, which are nonlinear. The default scheme, exponential fitting, is a choice of this package. Its properties and the other schemes are on the page Positivity and flux limiters.

DREAM's default is AD_INTERP_CENTRED (DREAM's py/DREAM/Settings/Equations/DistributionFunction.py). scheme=central matches it (Defaults that differ). DREAM's other choices include the flux limiters (Comparison with DREAM).

The stored Jacobian, or stored matrix, is the Jacobian of the implicit terms that have a stored Jacobian, which the solver assembles (The step). SparseLU factorizes the Newton matrix built from it, and the preconditioners of GMRES factorize the stored matrix in blocks (Linear solvers and preconditioners).

• jacobian=exact takes the exact Jacobian by automatic differentiation.
• jacobian=frozen_limiter freezes the limiter at the iterate (DREAM AD_INTERP_JACOBIAN_LINEAR).
• jacobian=upwind, central and exponential_fitting take the Jacobian of that linear scheme. upwind is DREAM AD_INTERP_JACOBIAN_UPWIND.

extrapolated_flux_consistent extrapolates the flux with the $\mathcal{V}'$ and the flux of each face, as DREAM's definition of BC_DPHI_CONST states (Comparison with DREAM). outflow lets only the advection leave the grid, and AdvectionDiffusion names no DREAM condition for it.

Notes on the solver

These notes belong to the rows preconditioner=, refresh=, restart= and operator= of the table above. The examples set three keywords of GMRES differently from its defaults.

operator= matters only for a term whose block of the stored Jacobian approximates its Jacobian (jacobian_kind), such as a limiter with a jacobian= other than exact. Without such terms the two values are the same. With operator=exact the approximate block only preconditions the exact Newton operator. operator=stored iterates with the stored matrix, a quasi-Newton iteration as in DREAM.

The switches of avalanche.jl

avalanche.jl builds its own grid and plasma, and it takes no checkpoint, stop or log switches. Of the shared switches it reads the schemes, pmax_boundary=, jacobian=, the solver switches, dt=, t_end=, saves= and output=. It reads mode=, screening=, lnlambda= and ion_size_table= with the values of its DREAM example as defaults. Its own switches are these.

Lower limits

nonnegative=true declares a LowerLimit of value 0 and tolerance 0 on $f$. NewtonSolver shortens the Newton update so that no limited entry loses more than the fraction max_reduction of its distance to its limit in one iteration. The default max_reduction = 0.9 is DREAM's value (Lower limits of the Newton step).

A step that would push an entry already at its limit further down fails with :below_lower_limit. FixedSteps then retries the step with half the time step. DREAM skips such entries instead.

At a low temperature on a wide grid, the Maxwell-Jüttner tail underflows to exact zeros, so such a run rejects and halves many steps. Shorter steps do not remove a round-off negative update of an entry at zero. A positive tolerance of the LowerLimit does. No switch of the examples sets the tolerance or max_reduction.

The stop of a segment

A segment stops after an accepted step when the file of stop_file= exists (StopCondition, Checkpoints, segments and the stop). It also stops when the wall-clock time left until its end is less than a margin. That end is the end of wall_budget= or the end of the Slurm job. The margin is the longest accepted step of the segment so far, plus the longest write of a checkpoint with its output, plus slack= seconds. wall_budget= counts its seconds from the start of the timed run, after the warm-up over the first three time points (What every script prints).

In a Slurm job, StopCondition takes SLURM_JOB_END_TIME as the end of the job. A run with checkpoint= stops before the end of the job without further switches. Without checkpoint=, the end of the job counts only when stop_file= or wall_budget= is given. With wall_budget= in a Slurm job, the earlier of the end of wall_budget= and the end of the job counts.

The settings of a run

integrate records the settings of the run with run_settings before the first step and returns them in settings of the TimeSolution. Each object of the package the run is built from has a describe record, and the record collects them by section.

• The sections of the setup are grids, layout, system, updates, prescribed, solver, time and run. Each has a settings_hash, and the record has a hash of them together.
• run holds the description that a script passes to integrate, such as the name of the example and the origin of its initial state.
• The sections code, environment, checkpoint, stop and log are outside the hashes.
• code holds the code_version: the commit of Runes2P when it runs from a git checkout, as the examples do, and "unknown" for an installed package.
• segment holds the number, the start, the status and the reason of a stop of this segment.
• complete is false when some object of the run has no record of its own, such as a term defined outside the package.

The record is taken before a stored state is restored, so restoring does not change it. A resumed segment whose setup differs lists the changed sections in its segment record.

Output files

output=run.h5 writes the solution with write_output: the saved states per unknown, the layout, the solver statistics and the settings record. read_output reads it back as a TimeSolution, and read_settings reads the settings record alone.

checkpoint=name writes, in checkpoint_dir=, the files that CheckpointSettings names.

• name.checkpoint.h5, the last checkpoint, and name.checkpoint.prev.h5, the one before it.
• name.segment<k>.h5 (k with four digits), the output of segment k.

A run that stops writes its checkpoint and exits with exit_status, so that a batch script can start the next segment with resume=true. read_output with the CheckpointSettings of the run joins its segments into one solution:

using Runes2P, HDF5
solution = read_output(CheckpointSettings("cooling"; directory = "run"))
f_end = unknown_view(solution.states[:, end], solution.layout, :f)   # f at the last save

The format of these files is on the page Output files.