Post Processing

For each run, an output directory is either created, overwritten, or appended to. This is determined by the command line arguments or by resumption from the input file (see Running Simulations for details).

Within each output directory, GUERNICA writes:

ex9.mesh — a copy of the mesh used for the run
git_diff.patch — the difference between the current code and the recorded git commit
<name>.cfg — a copy of the input file used for the run
run.log — a log file containing everything that would otherwise have been written to the terminal
gf_out/ — a directory containing saved simulation outputs as MFEM grid functions
diag/ — a directory containing end-of-run diagnostic fields (reaction-rate maps, the neutral flux vector, and target/wall flux profiles)
balance.dat — the global particle and energy budget, one row per output step
profile_files/ — the prescribed electron and ion background fields (density, temperature, velocity) projected onto the DG space, saved as .gf files (all-zero fields are skipped), plus copies of any file_rectilinear_2d input tables with a manifest.txt recording where each came from
operator_timing.txt — per-operator wall times, only when operator_timing = true

Output Controls

Output behavior is controlled through the [output] section of the input file.

`display`

display dictates whether anything is output to the log file.

`save_dist`

save_dist dictates whether the full distribution function is written for each intermediate time step, as set by vis_steps.

`save_final_dist`

save_final_dist dictates whether the full distribution function is written at the end of the simulation, independently of save_dist. Set it to true to leave behind a restartable state when the intermediate dumps are not wanted.

Resuming a simulation requires a complete set of distribution slices, so a run that is intended to be continued must set either save_dist = true or save_final_dist = true.

`profile`

profile = true suppresses all writing at intermediate time steps.

Saved outputs are still written at the end of the simulation even if profile = true.

`precision`

precision dictates the numerical precision used when outputs are written.

`track_balance`

track_balance dictates whether the global particle and energy budget is appended to balance.dat at every output step. It allocates no phase-space storage: on a steady run it reuses the residual the solver already computes, and on a transient run it costs one extra right-hand-side evaluation per output into scratch the solver is not using. It follows the same profile rule as every other output.

`operator_timing`

operator_timing writes a table of wall time per operator to run.log and operator_timing.txt at the end of the run. It synchronizes the device around every operator, so it is off by default. See Running Simulations.

Contents of `gf_out/`

Distribution Function

The distribution function is saved as velocity slices with filenames of the form:

ex9-v<iv>-<step>.gf

where:

iv is the flat velocity index
step is the time step index

The name carries identity only. Which velocity a slab holds is recorded in gf_out/velocity_grid.json, written once per run:

{
  "format": "guernica-velocity-grid",
  "version": 1,
  "Nv": 9,
  "grid": { "type": "cartesian", "Nx": 3, "Ny": 3, "Nz": 1 },
  "points": [ [0, -1.0, -1.0, 0.0], [1, 0.0, -1.0, 0.0], ... ]
}

Each entry of points is [iv, vx, vy, vz]. Look a slab up by velocity, not by index arithmetic — the construction of iv is private to the solver and will change when the velocity grid does. scripts/guernica_vgrid.py does this for you:

from guernica_vgrid import VelocityGrid

grid = VelocityGrid.load("runs/my_run/gf_out")
iv = grid.nearest(16.0, 15.0, 7.0)      # or .index_of(...) to require an exact node
path = grid.slab_path(iv, step=2000)

Within each file, the distribution function is stored in dof-inner, element-outer format.

Runs produced before this naming change used ex9-v<ix>_<iy>_<iz>-<step>.gf and have no manifest. They can still be resumed — the solver reads both names — and VelocityGrid.slab_path falls back to the old name when the new one is absent.

Density, Momentum, and Energy

Density, momentum, and energy are saved with filenames of the form:

rho-<step>.gf
rho_ux-<step>.gf
rho_uy-<step>.gf
rho_uz-<step>.gf
E-<step>.gf

Within each file, the values are stored in dof-inner, element-outer format.

Global Moments

The global mass, momentum, energy, and entropy values are written to:

moment.gf

They are saved first in the order:

mass
momentum_x
momentum_y
momentum_z
energy
entropy

and then in time.

moment.gf is only written at the end of a simulation. It is not populated while the simulation is running, so prematurely ending the simulation will result in moment.gf remaining empty.

moment.gf is currently only written when neutral-neutral collisions are on.

Contents of `diag/`

These fields are written once, after the solver has finished, and are intended for post-processing rather than for restarts. Nothing here is written at the vis_steps cadence, and nothing is computed inside the time or pseudo-time loop. The whole directory is skipped when profile = true.

The moment diagnostics reuse a state-sized work vector that the solver has already finished with, so they add no phase-space storage to a run.

Neutral moments and flux vector

File	Meaning
`n_n.gf`	neutral density, `∫ f dv`
`Gamma_n_x.gf`, `Gamma_n_y.gf`, `Gamma_n_z.gf`	neutral particle flux vector `Γ_n = n_n u_n = ∫ v f dv`
`E2.gf`	`∫ \\|v\\|² f dv`, i.e. twice the energy density
`Gamma_n_H1_x.gf`, `Gamma_n_H1_y.gf`	`Γ_n` projected onto a continuous H1 space

Γ_n is the raw first velocity moment, so it is the flux itself and not n_n multiplied by a separately stored u_n. The H1 fields exist only because streamline integrators (ParaView, matplotlib.streamplot) show visible jumps at element faces when handed raw DG data; they carry no extra information.

Reaction rate maps

For each active collision or source channel, the velocity moments of that operator's right-hand side are written using the same quadrature as the solution moments:

File	Meaning
`<channel>_particle.gf`	`Σ_j C(x, v_j) w_j`
`<channel>_mom_x.gf`, `_mom_y.gf`, `_mom_z.gf`	`Σ_j v_j C(x, v_j) w_j`
`<channel>_energy.gf`	`Σ_j ½\\|v_j\\|² C(x, v_j) w_j`

The channel prefixes are Rcx (whichever charge-exchange model is active), Riz (ionization), Rrec (recombination), and Rnn (BGK neutral–neutral).

The sign is the sign of the operator's contribution to df/dt, so a sink is negative. Rcx_particle should be near zero pointwise, since charge exchange conserves neutral particle number locally; it is written as a conservation check.

When ionization is active, nu_iz.gf holds the nodal ionization frequency and Riz_n.gf holds the positive sink rate ν_iz · n_n. Keeping nu_iz means the rate can be reformed against a different density snapshot without rerunning.

Target and wall flux profiles

wall_flux.dat is a plain-text table with one row per boundary-face GLL point:

attr face_id x y nx ny ds Gamma_i_wall alpha_w Gamma_n_emit q_n_emit Gamma_n_reflect q_n_reflect Gamma_n_incident q_n_incident

nx, ny is the outward unit normal and the Gamma_* columns are flux densities along it, so a face-integrated rate is Σ ds · Γ.

Gamma_i_wall — n_i u_i · n̂, the ion flux driving recycling
alpha_w — the wall emission amplitude actually used by the boundary condition
Gamma_n_emit — α_w Σ_{v·n<0} \|v·n\| M_w w_v, the emitted neutral flux
q_n_emit — the same sum with a ½\|v\|² weight
Gamma_n_reflect — Σ_{v·n<0} \|v·n\| c_spec f_spec w_v, the specularly reflected flux
q_n_reflect — the same sum with a ½\|v\|² weight
Gamma_n_incident — Σ_{v·n>0} (v·n) f_n⁻ w_v, using the same interior trace the upwind flux uses
q_n_incident — the same sum with a ½\|v\|² weight

The local recycling coefficient is R(s) = Gamma_n_emit / Gamma_n_incident.

The boundary condition assembles the ghost state as f_ghost = c_spec · f_spec + α_w M_w, so the emitted and reflected columns are the two halves of it and both are gains for the domain. The net wall source in a global budget is therefore

Σ ds · (Gamma_n_emit + Gamma_n_reflect − Gamma_n_incident)

and the same combination of the q_* columns for energy. A purely specular wall has Gamma_n_emit = 0 and Gamma_n_reflect = Gamma_n_incident, so it contributes nothing net, as it must.

Arc length s is deliberately not written: the point coordinates and face IDs are enough to order the profiles and split them per divertor leg in post-processing, and the choice of origin (for example a strike point) belongs in the plotting script. scripts/plot_wall_flux.py does this.

For type = ion_recycling, alpha_w is normalized per face against the face-integrated ion flux, so it is constant across the points of one face while Gamma_i_wall above is the pointwise value. The normalization uses the same surface rule as the advection operator, the collocated GLL points, so the face-integrated emitted flux the solver produces is exactly the one the amplitude was normalized to give:

Σ ds · Gamma_n_emit  ==  recycling_coeff · Σ ds · max(Gamma_i_wall, 0)

per face, which scripts/plot_wall_flux.py reports alongside the profiles.

Contents of `balance.dat`

balance.dat holds the global neutral particle and energy budget, one row per output step, so that the approach to a closed budget can be plotted next to the residual histories. The GLL-DGSEM volume and face terms telescope exactly against a constant test function, so at a discrete steady state

Σ_channels ∫ (moment of that operator's RHS) dV
  + ∮_wall (emitted + reflected − incident) ds
  + plasma-interface leakage
  = 0

once for the particle weight 1 and once for the energy weight ½|v|².

That same telescoping is what makes the running budget cheap. The left-hand side above is nothing other than the constant-test-function projection of the residual,

net = ∫ ( ∫ F(U) dv ) dV,

so net can be read straight off the moment fields the steady solver already forms when it measures the raw residual, without ever summing the individual channels. Summing them would need a phase-space scratch array held for the whole run; this needs two volume integrals over fields that already exist. On a transient march no residual is computed, so F(U) is evaluated once per output into the solver's stage derivative, which is dead between steps — one extra right-hand-side evaluation, still no extra storage.

The columns are:

step time N_total U_total
wall_emit_particle wall_reflect_particle wall_incident_particle leak_particle
net_particle scale_particle rel_particle
wall_emit_energy wall_reflect_energy wall_incident_energy leak_energy
net_energy scale_energy rel_energy

N_total, U_total — the domain integrals ∫ n_n dV and ∫ ½∫|v|²f dv dV
wall_*, leak_* — the boundary half of the budget, with wall_incident_* the positive flux leaving the domain
net_* — the imbalance of the whole budget, which goes to zero as the run converges
scale_* — |wall_emit| + |wall_reflect| + |wall_incident| + |leak|, the boundary throughput
rel_* — net/scale, the closure error relative to that throughput

The wall and leak columns are not a complete decomposition — the volume source and sink channels are not summed here — so they do not add up to net. The full per-channel breakdown is written once at the end of a run into diag/, where the solver's own scratch is dead and can be borrowed, so it also costs no extra storage. scripts/plot_balance.py sums that breakdown, which arrives at net by a route sharing no code with this one; comparing the two is a check on both.

The energy columns carry the ½|v|² weight, matching the diag/<channel>_energy convention rather than the E2 = ∫|v|²f dv field.

The file is flushed every row so a killed run still leaves a plottable history, and is appended to on resume.

Every row is also echoed to run.log as a one-liner giving net and rel for each weight:

  balance: particle net = -3.117e-04 (rel 1.204e-03), energy net = ...

On a steady march that line is printed directly under the residuals of the same iteration, so the closure error and the residual that produced it are read together.

Metadata

normalization.txt records n0, t0, L0, T0, m0 and v0 together with the sign and factor conventions above, so that the rates can be re-dimensionalized in post-processing.

Plotting

Examples of plotting results with Python are provided in the scripts/ directory.

Four scripts read diag/ directly. They share scripts/guernica_diag.py, which parses the mesh and the .gf files, reconstructs the nodal DG fields, builds the GLL diagonal mass matrix so volume integrals use the same quadrature as the solver, and converts to physical units. All of them take a run directory (or its diag/ directory), --save, and --normalized.

Units

The scripts report physical units by default, reading normalization.txt. Pass --normalized to see the solver's own values instead.

Field	Physical unit
`n_n`	m^-3
`Gamma_n_*`	m^-2 s^-1
`E2`	eV m^-3
`nu_iz`	s^-1
`Riz_n`, `<channel>_particle`	m^-3 s^-1
`<channel>_energy`	eV m^-3 s^-1
`<channel>_mom_*`	N m^-3 (force density)
wall `Gamma_*`	m^-2 s^-1
wall `q_n_emit`	eV m^-2 s^-1

Two subtleties are worth knowing before quoting numbers.

The energy factor of two. The neutral is hydrogen (m0 = 1 u) and v0 = sqrt(e · 1 eV / u), so m0 v0^2 = 1 eV exactly. E2 is int |v|^2 f dv with no one-half, i.e. twice the energy density, so it converts with n0 / 2 to eV/m³. As a check, a Maxwellian at rest has E2 = 3 n T and energy density (3/2) n T. The <channel>_energy fields and the wall q_n_* columns already carry the one-half, so they convert with no extra factor: n0 / t0 and n0 v0 respectively.

Integrated quantities are per metre. Configuration space is a 2D (R, Z) plane, so integrating a per-m³ rate over the domain gives m^-1 s^-1, i.e. a rate per metre in the ignored third direction, not a particle count. The solver has no axisymmetric 2*pi*R factor in its advection operator, so multiplying by one afterwards would not be consistent with what was solved.

Script	What it draws
`plot_rate_maps.py`	reaction-rate and moment maps on the (R, Z) mesh
`plot_neutral_flux.py`	`Γ_n` streamlines and quiver over a scalar background
`plot_balance.py`	the global particle balance table and bar chart
`plot_wall_flux.py`	target/wall profiles versus arc length

A fifth, plot_residuals.py, reads run.log rather than diag/ and is described at the end of this page.

python scripts/plot_rate_maps.py runs/<name> --preset rates

plot_rate_maps.py picks a colour normalization per field, since rates span many decades and often change sign: symmetric log for signed fields, log for single-signed ones, --scale linear to override. An all-negative field is drawn as |f| in a different colormap so that bright still means large. On log scales, values at or below the colour floor are clamped to it rather than masked, so exact zeros render as the darkest colour instead of leaving white holes; --floor-quantile moves that floor. --preset selects a group (rates, moments, ionization, momentum), --fields names fields explicitly, and --list shows what a run actually wrote.

python scripts/plot_neutral_flux.py runs/<name> --mode stream

plot_neutral_flux.py reconstructs the DG flux inside every element and interpolates through a triangulation built only from within-element sub-cells, so no triangle bridges a concave boundary or the private-flux gap. Coincident samples on shared element faces are averaged, which is the same arithmetic averaging the Gamma_n_H1_* fields carry. It does not need those files; they remain the route for ParaView. --background selects the scalar underneath (default n_n).

Streamline width follows |Gamma_n|, which typically spans four or more decades. --lw-scale chooses how:

mode	behaviour
`log` (default)	width tracks `log\|Gamma_n\|`; the only mode that separates a field spanning decades
`linear`	width proportional to `\|Gamma_n\|`; almost everything sits at the minimum
`power`	normalized magnitude to the `--lw-gamma` power: below 1 lifts weak flux into view, above 1 pushes contrast to the peaks

--lw-range MIN MAX sets the widths themselves (default 0.3 2.5), and --lw-clip LO HI sets the percentiles used as the endpoints (default 1 99), so one hot cell cannot set the scale. The script prints the flux percentiles with the width each maps to, which is the quickest way to tune it. --stream-color speed colours the lines by |Gamma_n| as well.

python scripts/plot_balance.py runs/<name>

plot_balance.py integrates every channel's particle moment over the domain, adds the wall emission, reflection and absorption from wall_flux.dat, and reports how close the sum is to zero. Since the terms come from unrelated parts of the solver, the closure is a useful check on both a converged run and the diagnostics. The energy budget is reported the same way and closes on the same terms.

python scripts/plot_balance.py runs/<name> --history

--history reads balance.dat instead of diag/ and plots the closure error and the domain inventory against output step. Both panels are dimensionless — the closure error is already a ratio and the inventories are divided by their final value — so the figure needs no normalization file and sits directly alongside a residual history.

python scripts/plot_wall_flux.py runs/<name> --origin 0.9 -1.5

plot_wall_flux.py chains boundary faces through their shared nodes into connected wall segments, builds the arc length along each, and plots the flux profiles and the local recycling coefficient. --origin picks which end of a segment is s = 0.

Residual histories

scripts/plot_residuals.py parses the steady solver's status lines out of run.log and plots the f, rho and E residual norms against iteration:

python scripts/plot_residuals.py runs/<name>

The solver reports the raw steady operator residual F(U) in two flavours per checked iteration, and --stage selects which to draw (all gives one panel each):

stage	meaning
`local` (default)	max over elements of the element-local relative L2 norm, for `rho` and `E`
`global`	relative L2 over the whole phase space, for `f`

There is no element-local norm for the distribution itself, and no global norm is kept for the moments, so each panel carries only the quantities listed above. Convergence against steady_tol is declared on the two local norms together with the global f norm; because these are the unscaled operator residual, they say whether the equations are satisfied rather than whether the solver has stopped taking steps.

--x puts the horizontal axis on iter (default), pseudo_time, walltime (cumulative, summed from the per-report intervals) or record. --tol draws a reference line. Several runs can be passed at once and are overlaid with different dash patterns, which is the direct way to compare runs with different settings. The script also prints first and last values with the reduction factor for each quantity.

Note that the log records a status line only every [time] vis_steps iterations (plus the final one), so that sets the resolution of these curves.

Post Processing

Output Controls

display

save_dist

save_final_dist

profile

precision

track_balance

operator_timing

Contents of gf_out/