Skip to content

Ionization Operator (IonizationOperator)

This page documents the electron-impact ionization operator in GUERNICA and how it maps from the paper model to the actual implementation in IonizationOperator.hxx/.cxx, and how it is wired into the main driver (guernica.cxx).

1. Model being implemented

In the paper model, electron-impact ionization is treated as a pure loss process for neutrals. Under the standard edge-plasma approximation that electron thermal speeds dominate neutral speeds, the operator reduces to a linear relaxation/loss term

\[ C_{\mathrm{iz}}[f_n](v) = -\nu_{\mathrm{iz}}(x)\, f_n(v), \qquad \nu_{\mathrm{iz}}(x) = n_e(x)\,\langle \sigma_{\mathrm{iz}} v_e\rangle\!\left(T_e(x)\right). \]

In the discrete-velocity formulation used by GUERNICA, the ionization term is evaluated pointwise in velocity (i.e., it does not couple different velocity nodes):

\[ C_{\mathrm{iz}}[f_n](v_j) = -n_e\langle \sigma_{\mathrm{iz}} v_e\rangle(T_e)\, f_{n,j}. \]

So conceptually: multiply the distribution function by a (space-dependent) scalar rate and subtract it from the RHS.

2. What the code actually applies

2.1 Operator form

The code implements exactly the “multiply by \(-\nu\)” structure, but written in terms of the global kinetic state vector U:

\[ \frac{dU}{dt} = -\nu(x)\, U, \]

where:

  • \(\nu(x)\) is represented as an mfem::Coefficient (so it can be spatially varying),
  • the action is local in configuration space DOFs and local in velocity index (no face fluxes, no velocity integrals).

This is explicitly stated in the header comment:

dU/dt = -nu(x) * U (no face terms; diagonal in velocity & dofs)

2.2 Discretization level

This operator is applied at the configuration-space DOF level, not quadrature-point level:

  • nu_coeff_ is projected onto a GridFunction nu_gf_ on the same FiniteElementSpace as the DG unknown.
  • per element, the nodal DOF values of nu_gf_ are packed into a contiguous array nu_dof_dev_.
  • the kernel then multiplies U DOFs by the corresponding packed nu DOFs and accumulates into the RHS.

This makes ionization “cheap” and easy to keep in the element-major layout: it’s a local multiply-add.

3. Data layout and indexing (how it ties into GUERNICA’s global state)

GUERNICA stores the kinetic unknown in one big vector with an element-major layout (as described in the implementation section of the paper).

IonizationOperator matches that layout:

  • For each element e, let ld = fes.GetFE(e)->GetDof() be the number of local DG DOFs.
  • For each discrete velocity index iv ∈ [0, Nv), the unknown block is contiguous:
\[ U_{e,iv} \in \mathbb{R}^{ld}. \]

Internally the operator precomputes:

  • ldof_e_[e] = ld
  • elem_base_[e] = base offset into the global kinetic vector for the start of element e (for velocity iv=0)
  • the per-velocity offset is base + iv*ld

For the ionization rate:

  • nu_gf_ stores \(\nu\) at FES DOFs (collocated),
  • nu_dof_dev_ stores per-element packed DOF values (length sum_e ldof_e[e]),
  • nu_off_[e] gives the offset of element e’s packed DOFs.

So for each (e, iv) and local dof i:

\[ \mathrm{RHS}[e,iv,i] \mathrel{+}= -\nu[e,i]\; U[e,iv,i]. \]

4. Code walk-through (what to look at)

4.1 Public interface (IonizationOperator.hxx)

Key methods:

  • AddRHS(U, dUdt): accumulates the ionization contribution into an existing RHS.
  • Mult(U, dUdt): standalone “apply operator” form; it zeroes dUdt and then calls AddRHS.
  • ProjectNu(): (re)projects the Coefficient onto nu_gf_ and repacks to device.

Practical meaning:

  • If your plasma background changes in time (e.g., ne(x,t), Te(x,t)), you should update the coefficient and call ProjectNu() when it changes.

4.2 Implementation details (IonizationOperator.cxx)

There are three implementation ideas worth remembering:

  1. Projection step
    • nu_gf_.ProjectCoefficient(nu_coeff_)
    • PackNuToDevice() packs element DOFs into nu_dof_dev_
  2. Device-friendly indexing
    • BuildDeviceMirrors() creates device arrays for integer indexing (ldof_e_dev_, elem_base_dev_, nu_off_dev_)
  3. Parallel kernel
    • loops over n = 0..NE*Nv-1 (interpreted as (e,iv) pairs),
    • within each, loops over local DOFs i = 0..ld-1 and does the multiply-add.

This is consistent with the statement in the paper that ionization is a local operator applied without data rearrangement.

5. How it is wired into the main solver (guernica.cxx)

In guernica.cxx, ionization is currently enabled by a boolean config option and (for now) uses a constant coefficient:

  • ionization = config.Get<bool>("ionization", false);
  • nu0 = config.Get<double>("nu0", 0.0);

If enabled, GUERNICA does:

```cpp mfem::ConstantCoefficient nu_coeff(nu0); izOp = std::make_shared(fes, Nv, nu_coeff, t_final); rhs.Add(izOp);