Methods Overview
GUERNICA solves the kinetic equation for the neutral distribution function \(f_n(\mathbf x, \mathbf v)\) in two spatial dimensions and two or three velocity dimensions,
where the right-hand side collects the active operators (ionization, recombination, charge exchange, neutral–neutral collisions, sources) and boundary conditions enter through the advection face terms.
Configuration space
In configuration space, GUERNICA uses a nodal discontinuous Galerkin spectral element method (DGSEM) on quadrilateral meshes. Each element uses a tensor-product basis on Gauss–Lobatto–Legendre (GLL) nodes, and the same nodes are used as quadrature points, so the mass matrix is diagonal. The polynomial order is set with order under [grid] (order >= 1), giving order + 1 nodes per element in each direction. Faces use an upwind flux. See Advection for details, and Hesthaven and Warburton, Nodal Discontinuous Galerkin Methods for background.
Velocity space
In velocity space, GUERNICA uses the discrete velocity method (DVM). The velocity domain is a uniform Cartesian grid, with the number of points and the bounds in each direction set in [vx], [vy], and [vz]. Setting num = 1 on an axis collapses it to a single node at zero, which is the usual way to run a 2D2V problem ([vz] num = 1).
Velocity moments are computed with a product of 1D quadrature rules:
- odd
num: composite Simpson's 1/3 rule, - even
num(≥ 4): the average of the two mirrored combinations of Simpson's 3/8 rule on three intervals and Simpson's 1/3 rule on the rest, which keeps the weights symmetric, num = 2: trapezoid rule,num = 1: weight 1.
Every operator uses these same weights, so discrete moments, sources, and collision terms are consistent with each other.
Units
GUERNICA works in normalized units. Temperatures are in eV, lengths in metres (\(L_0 = 1\) m), densities in units of \(n_0 = 10^{18}\ \mathrm{m^{-3}}\), and velocities in units of
where \(m\) is the neutral mass set by [units] mass_amu (default 1, giving \(v_0 \approx 9823\) m/s). The time unit is \(t_0 = L_0 / v_0\). With these choices a Maxwellian of temperature \(T\) in eV is \(\exp(-|\mathbf v|^2 / 2T)\). The normalization of each run is recorded in diag/normalization.txt.
Time integration
GUERNICA has two run modes, selected by mode under [time].
Transient (mode = transient)
The solution is advanced in physical time with an explicit third-order strong-stability-preserving Runge–Kutta method (SSP-RK3, solver = rk3). The run ends at t_final, with step dt, and outputs every vis_steps steps. All operators, including collisions, are treated explicitly.
The step is limited by a CFL condition. At startup GUERNICA prints an estimate of the stable step,
where \(p\) is the polynomial order and \(r_e\) is the outflow rate of element \(K_e\). The log line CFL estimate: dt_max = ... (dt = ..., ratio = ...) is followed by a note when dt exceeds the bound and a warning when it exceeds twice the bound. Runs above about twice the bound have been observed to diverge. This estimate covers advection only, so stiff collision rates can require a smaller step.
Steady state (mode = steady)
The steady solve marches in pseudo-time with a forward-Euler update (solver = pseudo_euler),
where \(F(U)\) is the full right-hand side (advection, boundaries, sources, and collisions). The pseudo-time step is local to each element and in-plane velocity:
Here dt is a reference step and \(r_e\) is the element outflow rate defined above. The minimum, mean, and maximum \(\Delta\tau\) are printed at startup.
Every check_every iterations the unscaled residual \(F(U)\) is evaluated in three norms: the element-local relative \(L^2\) norms of the density and energy residuals (maximum over elements), and the global relative \(L^2\) norm of the distribution residual. The run stops when the largest of the three falls below steady_tol, or after max_iters iterations. The final state is always written. If the state becomes non-finite, or its norm grows by more than \(10^8\) relative to the first check, the run stops with a divergence error asking for a smaller dt.
See Post Processing for plotting residual histories.