CDR plasma model

In the CDR plasma model we are solving

(12)\[\begin{split}&\nabla\cdot\left(\epsilon_r\nabla\Phi\right) = -\frac{\rho}{\epsilon_0},\\ &\frac{\partial\sigma}{\partial t} = F_\sigma,\\ &\frac{\partial n}{\partial t} + \nabla\cdot\left(\mathbf{v} n - D\nabla n\right) = S,\end{split}\]

The above equations must be supported by additional boundary conditions on electrodes and insulating surfaces.

Radiative transport can be done either in the diffusive approximation or by means of Monte Carlo methods (see RtSolver). Diffusive RTE methods (see EddingtonSP1) involve solving

\[\partial_t\Psi + \kappa\Psi - \nabla\cdot\left(\frac{1}{3\kappa}\nabla\Psi\right) = \frac{\eta}{c},\]

where $\Psi$ is the isotropic photon density, $\kappa$ is an absorption length and $\eta$ is an isotropic source term. I.e., $\eta$ is the number of photons produced per unit time and volume. The time dependent term can be turned off and the equation can be solved stationary.

The module also supports discrete photons where photon transport and absorption is done by sampling discrete photons (see Monte Carlo solver). In general, discrete photon methods incorporate better physics (like shadows). They can also be more easily adapted to scattering media. They are, on the other hand, inherently stochastic which implies that some extra caution must be exercised when integrating the equations of motion.

The coupling that is (currently) available in chombo-discharge is

(13)\[\begin{split}\epsilon_r =& \epsilon_r(\mathbf{x}),\\ \mathbf{v} =& \mathbf{v}\left(t, \mathbf{x}, \mathbf{E}, n\right),\\ D =& \mathbf{v}\left(t, \mathbf{x}, \mathbf{E}, n\right),\\ S =& S\left(t, \mathbf{x}, \mathbf{E}, \nabla\mathbf{E}, n, \nabla n, \Psi\right),\\ \eta =& \eta\left(t, \mathbf{x}, \mathbf{E}, n\right),\\ F =& F(t, \mathbf{x}, \mathbf{E}, n),\end{split}\]

where $F$ is the boundary flux on insulators or electrodes (which must be separately implemented).

chombo-discharge works by embedding the equations above into an abstract C++ framework (see CdrPlasmaPhysics) that the user must implement or reuse existing pieces of, and then compile into an executable.

Tip

The CDR plasma model resides in /Physics/CdrPlasma and describes plasmas in the drift-diffusion approximation. This physics model also includes the following subfolders:

/Physics/CdrPlasma/PlasmaModel which contains various implementation of some plasma models that we have used.
/Physics/CdrPlasma/TimeSteppers contains various algorithms for advancing the equations of motion.
/Physics/CdrPlasma/CellTaggers contains various algorithms for flagging cells for refinement and coarsening.

See CdrPlasmaStepper for the overall C++ API.

This module uses the following solvers:

Advection-diffusion-reaction solver, CdrSolver.
Electrostatics solvers, FieldSolver.
Radiative transfer solver (either Monte-Carlo or continuum approximation), RtSolver.
Surface charge solver, see Surface ODE solver.

Time discretizations

Here, we discuss two discretizations of Eq. 12. Firstly, note that there are two layers to the time integrators:

A pure class CdrPlasmaStepper which inherits from TimeSteppers but does not implement an advance method. This class simply provides the base functionality for more easily developing time integrators. CdrPlasmaStepper contains methods that are necessary for coupling the solvers, e.g. calling the CdrPlasmaPhysics methods at the correct time.
Implementations of CdrPlasmaPhysics, which implement the advance method and can thus be used for advancing models.

The supported time integrators are located in $DISCHARGE_HOME/CdrPlasma/TimeSteppers. There are two integrators that are commonly used.

A Godunov operator splitting with either explicit or implicit diffusion. This integrator also supports semi-implicit formulations.
A spectral deferred correction (SDC) integrator with implicit diffusion. This integrator is an implicit-explicit.

Briefly put, the Godunov operator is our most stable integrator, while the SDC integrator is our most accurate integrator.

Godunov operator splitting

The CdrPlasmaGodunovStepper implements CdrPlasmaStepper and defines an operator splitting method between charge transport and plasma chemistry. It has a formal order of convergence of one. The source code is located in $DISCHARGE_HOME/Physics/CdrPlasma/TimeSteppers/CdrPlasmaGodunovStepper.

Warning

Splitting the terms yields splitting errors which can dominate for large time steps. Typically, the operator splitting discretization is not suitable for large time steps.

The basic advancement routine for CdrPlasmaGodunovStepper is as follows:

Advance the charge transport $\phi^k \rightarrow \phi^{k+1}$ with the source terms set to zero.
Compute the electric field.
Advance the plasma chemistry over the same time step using the field computed above I.e., advance $\partial_t\phi = S$ over a time step $\Delta t$.
Advance the radiative transport part. This can also involve discrete photons.

The transport/field steps can be done in various ways: The following transport algorithms are available:

Euler, where everything is advanced with the Euler rule.
Semi-implicit, where the Euler field/transport step is performed with an implicit coupling to the electric field.

In addition, diffusion can be treated

Explicitly, where all diffusion advances are performed with an explicit rule.
Implicitly, where all diffusion advances are performed with an implicit rule.
Automatically, where diffusion advances are performed with an implicit rule only if time steps dictate it, and explicitly otherwise.

Note

When setting up a new problem with the Godunov time integrator, the default setting is to use automatic diffusion and a semi-implicit coupling. These settings tend to work for most problems.

Specifying transport algorithm

To specify the transport algorithm, modify the flag CdrPlasmaGodunovStepper.transport, and set it to semi_implicit or euler. Everything else is an error.

Note that for the Godunov integrator, it is possible to center the advective discretization at the half time step. That is, the advancement algorithm is

\[n^{k+1} = n^{k} - \nabla\cdot\left(n^{k+1/2}\mathbf{v}\right) + \nabla\cdot\left(D\nabla\phi^k\right),\]

where $n^{k+1/2}$ is obtained by also including transverse slopes (i.e., extrapolation in time). See Trebotich and Graves [2015] for details. Note that the formal order of accuracy is still one, but the accuracy of the advective discretization is increased substantially.

Specifying diffusion

To specify how diffusion is treated, modify the flag CdrPlasmaGodunovStepper.diffusion, and set it to auto, explicit, or implicit. In addition, the flag CdrPlasmaGodunovStepper.diffusion_thresh must be set to a number.

When diffusion is set to auto, the integrator switches to implicit diffusion when

\[\frac{\Delta t_{\textrm{A}}}{\Delta t_{\textrm{AD}}} > \epsilon,\]

where $\Delta t_{\textrm{A}}$ is the advection-only limited time step and $\Delta t_{\textrm{AD}}$ is the advection-diffusion limited time step.

Note

When there are multiple species being advected and diffused, the integrator will perform extra checks in order to maximize the time steps for the other species.

Time step limitations

The basic time step limitations for the Godunov integrator are:

Manually set maximum and minimum time steps
Courant-Friedrichs-Lewy conditions, either on advection, diffusion, or both.
The dielectric relaxation time.

The user is responsible for setting these when running the simulation. Note when the semi-implicit scheme is used, it is not necessary to restrict the time step by the dielectric relaxation time.

Solver configuration

CdrPlasmaGodunovStepper contains several options that can be configured when running the solver, which are listed below:

# ====================================================================================================
# CdrPlasmaGodunovStepper class options
# ====================================================================================================
CdrPlasmaGodunovStepper.verbosity             = -1               ## Class verbosity
CdrPlasmaGodunovStepper.solver_verbosity      = -1               ## Individual solver verbosities
CdrPlasmaGodunovStepper.min_dt                = 0.               ## Minimum permitted time step
CdrPlasmaGodunovStepper.max_dt                = 1.E99            ## Maximum permitted time step
CdrPlasmaGodunovStepper.cfl                   = 0.8              ## CFL number
CdrPlasmaGodunovStepper.use_regrid_slopes     = false            ## Use slopes when regridding (or not)
CdrPlasmaGodunovStepper.filter_rho            = 0                ## Number of filterings of space charge
CdrPlasmaGodunovStepper.filter_compensate     = false            ## Use compensation step after filter or not
CdrPlasmaGodunovStepper.field_coupling        = semi_implicit    ## Field coupling. 'explicit' or 'semi_implicit'
CdrPlasmaGodunovStepper.advection             = muscl            ## Advection algorithm. 'euler', 'rk2', or 'muscl'
CdrPlasmaGodunovStepper.diffusion             = explicit         ## Diffusion. 'explicit', 'implicit', or 'auto'.
CdrPlasmaGodunovStepper.diffusion_thresh      = 1.2              ## Diffusion threshold. If dtD/dtA > this then we use implicit diffusion.
CdrPlasmaGodunovStepper.diffusion_order       = 2                ## Diffusion order.
CdrPlasmaGodunovStepper.relax_time            = 100.             ## Relaxation time. 100 <= is usually a "safe" choice.
CdrPlasmaGodunovStepper.fast_poisson          = 1                ## Solve Poisson every this time steps. Mostly for debugging.
CdrPlasmaGodunovStepper.fast_rte              = 1                ## Solve RTE every this time steps. Mostly for debugging.
CdrPlasmaGodunovStepper.fhd                   = false            ## Set to true if you want to add a stochastic diffusion flux
CdrPlasmaGodunovStepper.source_comp           = interp           ## Interpolated interp, or upwind X for species X
CdrPlasmaGodunovStepper.floor_cdr             = true             ## Floor CDR solvers to avoid negative densities
CdrPlasmaGodunovStepper.debug                 = false            ## Turn on debugging messages. Also monitors mass if it was injected into the system.
CdrPlasmaGodunovStepper.profile               = false            ## Turn on/off performance profiling.

Spectral deferred corrections

The CdrPlasmaImExSdcStepper uses implicit-explicit (ImEx) spectral deferred corrections (SDCs) to advance the equations. This integrator implements the advance method for CdrPlasmStepper, and is a high-order method with implicit diffusion.

SDC basics

First, we provide a quick introduction to the SDC procedure. Given an ordinary differential equation (ODE) as

\[\frac{\partial u}{\partial t} = F(u,t), \quad u(t_0) = u_0,\]

the exact solution is

\[u(t) = u_0 + \int_{t_0}^tF\left(u,\tau\right)d\tau.\]

Denote an approximation to this solution by $\widetilde{u}(t)$ and the correction by $\delta(t) = u(t) - \widetilde{u}(t)$. The measure of error in $\widetilde{u}(t)$ is then

\[R(\widetilde{u}, t) = u_0 + \int_{t_0}^tF(\widetilde{u}, \tau)d\tau - \widetilde{u}(t).\]

Equivalently, since $u = \widetilde{u} + \delta$, we can write

\[\widetilde{u} + \delta = u_0 + \int_{t_0}^t F\left(\widetilde{u}+\delta, \tau\right)d\tau.\]

This yields

\[\delta = \int_{t_0}^t\left[F\left(\widetilde{u}+\delta, \tau\right) - F\left(\widetilde{u}, \tau\right)\right]d\tau + R\left(\widetilde{u},t\right).\]

This is called the correction equation. The goal of SDC is to iteratively solve this equation in order to provide a high-order discretization.

The ImEx SDC method in chombo-discharge uses implicit diffusion in the SDC scheme. Coupling to the electric field is always explicit. The user is responsible for specifying the quadrature nodes, as well as setting the number of sub-intervals in the SDC integration and the number of corrections. In general, each correction raises the discretization order by one.

Time step limitations

The ImEx SDC integrator is limited by

The dielectric relaxation time.
An advective CFL conditions.

In addition to this, the user can specify maximum/minimum allowed time steps.

CdrPlasmaPhysics

Overview

CdrPlasmaPhysics is an abstract class which represents the plasma physics for the CDR plasma module, i.e. it provides the coupling functions in Eq. 13. The source code for the class resides in /Physics/CdrPlasma/CD_CdrPlasmaPhysics.H. Note that the entire class is an interface, whose implementations are used by the time integrators that advance the equations.

There are no default input parameters for CdrPlasmaPhysics, as users must generally implement their own kinetics. The class exists solely for providing the integrators with the necessary fundamentals for filling solvers with the correct quantities at the same time, for example filling source terms and drift velocities.

A successful implementation of CdrPlasmaPhysics has the following:

Instantiated a list of CdrSpecies. These become Convection-Diffusion-Reaction solvers and contain initial conditions and basic transport settings for the convection-diffusion-reaction solvers.
Instantiated a list RtSpecies. These become Radiative transfer solvers and contain metadata for the radiative transport solvers.
Implemented the core functionality that couple the solvers together.

chombo-discharge automatically allocates the specified number of convection-diffusion-reaction and radiative transport solvers from the list of species the is instantiated. For information on how to interface into the CDR solvers, see CdrSpecies. Likewise, see RtSpecies for how to interface into the RTE solvers.

Implementation of the core functionality is comparatively straightforward, but can lead to boilerplate code. For this reason we also provide an implementation layer JSON interface that provides a plug-and-play interface for specifying the plasma physics by using a JSON schema for description the physics.

Complete API

The full API for the CdrPlasmaPhysics class is given below:

Show/hide full 0D interface

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */

/**
   @file   CD_CdrPlasmaPhysics.H
   @brief  Declaration of the Physics::CdrPlasma::CdrPlasmaPhysics interface class
   @author Robert Marskar
*/

#ifndef CD_CDRPLASMAPHYSICS_H
#define CD_CDRPLASMAPHYSICS_H

// Chombo includes
#include <RealVect.H>
#include <RefCountedPtr.H>
#include <LoHiSide.H>

// Our includes
#include <CD_CdrSpecies.H>
#include <CD_RtSpecies.H>
#include <CD_NamespaceHeader.H>

namespace Physics {
  namespace CdrPlasma {

    /**
       @brief Abstract interface for specifying plasma kinetics in the CdrPlasma physics module.
       @details This class defines the coupling functions required by CdrPlasmaStepper to advance
       a system of convection-diffusion-reaction (CDR) equations for charged and neutral species
       together with radiative transfer equations (RTE) for photon densities, all in a
       self-consistent electric field obtained from Poisson's equation.

       Concrete subclasses must implement:
       - advanceReactionNetwork: compute CDR and RTE source terms from local state.
       - computeCdrDriftVelocities: drift velocities for each CDR species.
       - computeCdrDiffusionCoefficients: diffusion coefficients for each CDR species.
       - computeCdrElectrodeFluxes: EB boundary fluxes on electrode surfaces.
       - computeCdrDielectricFluxes: EB boundary fluxes on dielectric surfaces.
       - computeCdrDomainFluxes: fluxes through the computational domain boundary.
       - initialSigma: initial surface charge density on dielectrics.
       - computeAlpha: Townsend ionization coefficient (used by cell taggers).
       - computeEta: Townsend attachment coefficient (used by cell taggers).

       Species are registered by filling m_cdrSpecies and m_rtSpecies in the subclass constructor.
       The ordering of species in all input/output vectors throughout this interface matches the
       ordering in those two containers.
    */
    class CdrPlasmaPhysics
    {
    public:
      /**
         @brief Default constructor. Does nothing.
      */
      CdrPlasmaPhysics()
      {}

      /**
         @brief Destructor. Does nothing.
      */
      virtual ~CdrPlasmaPhysics()
      {}

      /**
         @brief Parse runtime options. Default implementation does nothing.
      */
      virtual void
      parseRuntimeOptions() {};

      /**
         @brief Return the number of plot variables provided by this physics class.
         @details CdrPlasmaStepper calls this to pre-allocate output storage. The returned value must
         match the length of the vectors returned by getPlotVariableNames() and getPlotVariables().
         @return Number of extra plot variables. Default is 0.
      */
      virtual int
      getNumberOfPlotVariables() const
      {
        return 0;
      }

      /**
         @brief Return the names of the extra plot variables.
         @details The length of the returned vector must match getNumberOfPlotVariables(). The
         ordering must be consistent with getPlotVariables().
         @return Vector of variable names.
      */
      virtual Vector<std::string>
      getPlotVariableNames() const
      {
        return Vector<std::string>(this->getNumberOfPlotVariables(), std::string("empty data"));
      }

      /**
         @brief Evaluate the extra plot variables at a single grid point.
         @details Called by CdrPlasmaStepper when writing plot files. The length of the returned
         vector must match getNumberOfPlotVariables().
         @param[in] a_cdrDensities  CDR species densities at this point.
         @param[in] a_cdrGradients  CDR species density gradients at this point.
         @param[in] a_rteDensities  RTE photon densities at this point.
         @param[in] a_E             Electric field vector at this point.
         @param[in] a_position      Physical position.
         @param[in] a_dx            Grid cell spacing.
         @param[in] a_dt            Time step used in the last advance.
         @param[in] a_time          Current simulation time.
         @param[in] a_kappa         Volume fraction of the grid cell.
         @return Vector of plot variable values.
      */
      virtual Vector<Real>
      getPlotVariables(const Vector<Real>&     a_cdrDensities,
                       const Vector<RealVect>& a_cdrGradients,
                       const Vector<Real>&     a_rteDensities,
                       const RealVect&         a_E,
                       const RealVect&         a_position,
                       const Real              a_dx,
                       const Real              a_dt,
                       const Real              a_time,
                       const Real              a_kappa) const
      {
        return Vector<Real>(this->getNumberOfPlotVariables(), 1.0);
      }

      /**
         @brief Compute the Townsend ionization coefficient alpha.
         @details Used primarily by cell tagging classes to identify high-ionization regions.
         @param[in] a_E        Electric field magnitude (V/m).
         @param[in] a_position Physical position.
         @return Townsend ionization coefficient (1/m).
      */
      virtual Real
      computeAlpha(const Real a_E, const RealVect& a_position) const = 0;

      /**
         @brief Compute the Townsend attachment coefficient eta.
         @details Used primarily by cell tagging classes to identify attachment-dominated regions.
         @param[in] a_E        Electric field magnitude (V/m).
         @param[in] a_position Physical position.
         @return Townsend attachment coefficient (1/m).
      */
      virtual Real
      computeEta(const Real a_E, const RealVect& a_position) const = 0;

      /**
         @brief Advance the reaction network over a time interval and return source terms.
         @details The advance is assumed to take the form phi^{k+1} = phi^k + S * dt. Implementations
         may either fill a_cdrSources and a_rteSources directly (explicit rule), or perform a
         fully implicit sub-step and back out the equivalent source S from the result.
         @param[out] a_cdrSources   CDR source terms (same ordering as m_cdrSpecies).
         @param[out] a_rteSources   RTE source terms (same ordering as m_rtSpecies).
         @param[in]  a_cdrDensities CDR species densities at this point.
         @param[in]  a_cdrGradients CDR species density gradients at this point.
         @param[in]  a_rteDensities RTE photon densities at this point.
         @param[in]  a_E            Electric field vector at this point.
         @param[in]  a_pos          Physical position.
         @param[in]  a_dx           Grid cell spacing.
         @param[in]  a_dt           Time step size.
         @param[in]  a_time         Current simulation time.
         @param[in]  a_kappa        Volume fraction of the grid cell.
      */
      virtual void
      advanceReactionNetwork(Vector<Real>&           a_cdrSources,
                             Vector<Real>&           a_rteSources,
                             const Vector<Real>&     a_cdrDensities,
                             const Vector<RealVect>& a_cdrGradients,
                             const Vector<Real>&     a_rteDensities,
                             const RealVect&         a_E,
                             const RealVect&         a_pos,
                             const Real              a_dx,
                             const Real              a_dt,
                             const Real              a_time,
                             const Real              a_kappa) const = 0;

      /**
         @brief Compute drift velocities for all CDR species.
         @param[in] a_time         Current simulation time.
         @param[in] a_pos          Physical position.
         @param[in] a_E            Electric field vector.
         @param[in] a_cdrDensities CDR species densities.
         @return Drift velocity vectors in the same order as m_cdrSpecies.
      */
      virtual Vector<RealVect>
      computeCdrDriftVelocities(const Real          a_time,
                                const RealVect&     a_pos,
                                const RealVect&     a_E,
                                const Vector<Real>& a_cdrDensities) const = 0;

      /**
         @brief Compute diffusion coefficients for all CDR species.
         @param[in] a_time         Current simulation time.
         @param[in] a_pos          Physical position.
         @param[in] a_E            Electric field vector.
         @param[in] a_cdrDensities CDR species densities.
         @return Scalar diffusion coefficients in the same order as m_cdrSpecies.
      */
      virtual Vector<Real>
      computeCdrDiffusionCoefficients(const Real          a_time,
                                      const RealVect&     a_pos,
                                      const RealVect&     a_E,
                                      const Vector<Real>& a_cdrDensities) const = 0;

      /**
         @brief Compute CDR species fluxes on electrode-gas EB interfaces.
         @details Used as a boundary condition on electrode EB cells. A positive flux is directed
         into the gas phase (outward from the electrode).
         @param[in] a_time            Current simulation time.
         @param[in] a_pos             EB centroid position.
         @param[in] a_normal          Outward unit normal pointing into the gas phase.
         @param[in] a_E               Electric field at the EB centroid.
         @param[in] a_cdrDensities    CDR densities at the EB centroid.
         @param[in] a_cdrVelocities   Normal component of CDR drift velocities at the EB.
         @param[in] a_cdrGradients    Normal gradient of CDR densities at the EB.
         @param[in] a_rteFluxes       Normal component of RTE fluxes at the EB.
         @param[in] a_extrapCdrFluxes CDR fluxes extrapolated from the gas-phase interior.
         @return Boundary fluxes in the same order as m_cdrSpecies.
      */
      virtual Vector<Real>
      computeCdrElectrodeFluxes(const Real          a_time,
                                const RealVect&     a_pos,
                                const RealVect&     a_normal,
                                const RealVect&     a_E,
                                const Vector<Real>& a_cdrDensities,
                                const Vector<Real>& a_cdrVelocities,
                                const Vector<Real>& a_cdrGradients,
                                const Vector<Real>& a_rteFluxes,
                                const Vector<Real>& a_extrapCdrFluxes) const = 0;

      /**
         @brief Compute CDR species fluxes on dielectric-gas EB interfaces.
         @details Used as a boundary condition on dielectric EB cells. A positive flux is directed
         into the gas phase (outward from the dielectric).
         @param[in] a_time            Current simulation time.
         @param[in] a_pos             EB centroid position.
         @param[in] a_normal          Outward unit normal pointing into the gas phase.
         @param[in] a_E               Electric field at the EB centroid.
         @param[in] a_cdrDensities    CDR densities at the EB centroid.
         @param[in] a_cdrVelocities   Normal component of CDR drift velocities at the EB.
         @param[in] a_cdrGradients    Normal gradient of CDR densities at the EB.
         @param[in] a_rteFluxes       Normal component of RTE fluxes at the EB.
         @param[in] a_extrapCdrFluxes CDR fluxes extrapolated from the gas-phase interior.
         @return Boundary fluxes in the same order as m_cdrSpecies.
      */
      virtual Vector<Real>
      computeCdrDielectricFluxes(const Real          a_time,
                                 const RealVect&     a_pos,
                                 const RealVect&     a_normal,
                                 const RealVect&     a_E,
                                 const Vector<Real>& a_cdrDensities,
                                 const Vector<Real>& a_cdrVelocities,
                                 const Vector<Real>& a_cdrGradients,
                                 const Vector<Real>& a_rteFluxes,
                                 const Vector<Real>& a_extrapCdrFluxes) const = 0;

      /**
         @brief Compute CDR species fluxes through computational domain boundaries.
         @details A positive flux is directed into the domain.
         @param[in] a_time            Current simulation time.
         @param[in] a_pos             Position on the domain face.
         @param[in] a_dir             Coordinate direction (0=x, 1=y, 2=z).
         @param[in] a_side            Low or high side of the domain.
         @param[in] a_E               Electric field at the domain face.
         @param[in] a_cdrDensities    CDR densities at the domain face.
         @param[in] a_cdrVelocities   Normal component of CDR drift velocities at the domain face.
         @param[in] a_cdrGradients    Normal component of CDR density gradients at the domain face.
         @param[in] a_rteFluxes       Normal component of RTE fluxes at the domain face.
         @param[in] a_extrapCdrFluxes CDR fluxes extrapolated from the interior.
         @return Boundary fluxes in the same order as m_cdrSpecies.
      */
      virtual Vector<Real>
      computeCdrDomainFluxes(const Real           a_time,
                             const RealVect&      a_pos,
                             const int            a_dir,
                             const Side::LoHiSide a_side,
                             const RealVect&      a_E,
                             const Vector<Real>&  a_cdrDensities,
                             const Vector<Real>&  a_cdrVelocities,
                             const Vector<Real>&  a_cdrGradients,
                             const Vector<Real>&  a_rteFluxes,
                             const Vector<Real>&  a_extrapCdrFluxes) const = 0;

      /**
         @brief Return the initial surface charge density on dielectric surfaces.
         @param[in] a_time Current simulation time (typically 0).
         @param[in] a_pos  Physical position on the dielectric surface.
         @return Initial surface charge density (C/m^2 in 3D, C/m in 2D).
      */
      virtual Real
      initialSigma(const Real a_time, const RealVect& a_pos) const = 0;

      /**
         @brief Return the list of CDR species.
         @return Const reference to the CDR species vector.
      */
      const Vector<RefCountedPtr<CdrSpecies>>&
      getCdrSpecies() const
      {
        return m_cdrSpecies;
      }

      /**
         @brief Return the list of RTE species.
         @return Const reference to the RTE species vector.
      */
      const Vector<RefCountedPtr<RtSpecies>>&
      getRtSpecies() const
      {
        return m_rtSpecies;
      }

      /**
         @brief Return the number of CDR species.
         @return Number of CDR species.
      */
      int
      getNumCdrSpecies() const
      {
        return m_cdrSpecies.size();
      }

      /**
         @brief Return the number of RTE equations.
         @return Number of RTE species.
      */
      int
      getNumRtSpecies() const
      {
        return m_rtSpecies.size();
      }

    protected:
      /**
         @brief CDR species definitions. Subclasses populate this in their constructor.
      */
      Vector<RefCountedPtr<CdrSpecies>> m_cdrSpecies;

      /**
         @brief RTE species definitions. Subclasses populate this in their constructor.
      */
      Vector<RefCountedPtr<RtSpecies>> m_rtSpecies;

      /**
         @brief Number of CDR species.
      */
      int m_numCdrSpecies;

      /**
         @brief Number of RTE species.
      */
      int m_numRtSpecies;
    };
  } // namespace CdrPlasma
} // namespace Physics

#include <CD_NamespaceFooter.H>

#endif

JSON interface

Since implementations of CdrPlasmaPhysics are usually boilerplate, we provide a class CdrPlasmaJSON which can initialize and parse various types of initial conditions and reactions from a JSON input file. This class is defined in $DISCHARGE_HOME/Physics/PlasmaModels/CdrPlasmaJSON.

CdrPlasmaJSON is a full implementation of CdrPlasmaPhysics which supports the definition of various species (neutral, plasma species, and photons) and methods of coupling them. We expect that CdrPlasmaJSON provides the simplest method of setting up a new plasma model. It is also comparatively straightforward to extend the class with further required functionality.

In the JSON interface, the radiative transfer solvers always solve for the number of photons that lead to photoionization events. This means that the interpretation of $\Psi$ is the number of photoionization events during the previous time step. This is true for both continuum and discrete radiative transfer models.

Usage

To use this plasma model, use -physics CdrPlasmaJSON when setting up a new plasma problem (see Setting up a new problem). When CdrPlasmaJSON is instantiated, the constructor will parse species, reactions, initial conditions, and boundary conditions from a JSON file that the user provides. In addition, users can parse transport data or reaction rates from tabulated ASCII files that they provide.

To specify the input plasma kinetics file, include

Specifying input file

CdrPlasmaJSON will read a JSON file specified by the input variable CdrPlasmaJSON.chemistry_file.

Discrete photons

There are two approaches when using discrete photons, and both rely on the user setting up the application with the Monte Carlo photon solver (rather than continuum solvers). For an introduction to the particle radiative transfer solver, see Monte Carlo solver.

The user must use one of the following:

Set the following class options:
CdrPlasmaJSON.discrete_photons = true McPhoto.photon_generation = deterministic McPhoto.source_type = number
When specifying CdrPlasmaJSON.discrete_photons = true, CdrPlasmaJSON will do a Poisson sampling of the number of photons that are generated in each cell and put this in the radiative transfer solvers’ source terms. This means that the radiative transfer solver source terms contain the physical number of photons generated in one time step. To turn off sampling inside the radiative transfer solver, we specify McPhoto.photon_generation = stochastic and set McPhoto.source_type = number to let the solver know that the source contains the number of physical photons.
Alternatively, set the following class options:
CdrPlasmaJSON.discrete_photons = false McPhoto.photon_generation = stochastic McPhoto.source_type = volume_rate
In this case the CdrPlasmaJSON class will fill the solver source terms with the volumetric rate, i.e. the number of photons produced per unit volume and time. When McPhoto generates the photons it will compute the number of photons generated in a cell through Poisson sampling $n = P\left(S_\gamma\Delta V\Delta t\right)$ where $P$ indicates a Poisson sampling operator.

Fundamentally, the two approaches differ only in where the Poisson sampling is performed. With the first approach, plotting the radiative transfer solver source terms will show the number of physical photons generated. In the second approach, the source terms will show the volume photo-generation rate.

Gas law and neutral background

General functionality

To include the gas law and neutral species, include a JSON object gas with the field law specified. Currently, law can be either ideal, troposphere, or table.

The purpose of the gas law is to set the temperature, pressure, and neutral density of the background gas. In addition, we specify the neutral species that are used through the simulation. These species are not stored on the mesh; we only store function pointers to their temperature, density, and pressure.

It is also possible to include a field plot which will then include the temperature, pressure, and density in plot files.

Ideal gas

To specify an ideal gas law, specify ideal gas law as follows:

{"gas":
  {
    "law": "ideal",
    "temperature": 300,
    "pressure": 1
  }
 }

In this case the gas pressure and temperatures will be as indicated, and the gas number density will be computed as

\[\rho = \frac{p_0^\prime N_{\textrm{A}}}{RT_0},\]

where $p^\prime$ is the pressure converted to Pascals.

Note that the input temperature should be specified in Kelvin, and the input pressure in atmospheres.

Troposphere

It is also possible to specify the pressure, temperature, and density to be functions of tropospheric altitude. In this case one must specify the extra fields

molar mass For specifying the molar mass (in $\textrm{g}\cdot\textrm{mol}^{-1}$) of the gas.
gravity Gravitational acceleration $g$.
lapse rate Temperature lapse rate $L$ in units of $\textrm{K}/\textrm{m}$.

In this case the gas temperature pressure, and number density are computed as

\[T(h) = T_0 - Lh\]

\[p(h) = p_0\left((1 - \frac{Lh}{T_0}\right)^{\frac{g M}{RL}}\]

\[\rho(h) = \frac{p^\prime(h) N_{\textrm{A}}}{RT(h)}\]

For example, specification of tropospheric conditions can be included by

{"gas":
  {
    "law": "troposphere",
    "temperature": 300,
    "pressure": 1,
    "molar_mass": 28.97,
    "gravity": 9.81,
    "lapse_rate": 0.0065,
    "plot": true
  }
}

Tabulated

To specify temperature, density, and pressure as function of altitude, set law to table and include the following fields:

file For specifying which file we read the data from.
height For specifying the column where the height is stored (in meters).
temperature For specifying the column where the temperature (in Kelvin) is stored.
pressure For specifying the column where the pressure (in Pascals) is stored.
density For specifying the column where the density (in $\textrm{kg}\cdot\textrm{m}^{-3}$) is stored.
molar mass For specifying the molar mass (in $\textrm{g}\cdot\textrm{mol}^{-1}$) of the gas.
min height For setting the minimum altitude in the chombo-discharge internal table.
max height For setting the minimum altitude in the chombo-discharge internal table.
res height For setting the height resolution in the chombo-discharge internal table.

For example, assume that our file MyAtmosphere.dat contains the following data:

# z [m]              rho [kg/m^3]    T [K]           p [Pa]
0000000E+00        1.2900000E+00   2.7210000E+02   1.0074046E+05
0000000E+03        1.1500000E+00   2.6890000E+02   8.8751220E+04
0000000E+03        1.0320000E+00   2.6360000E+02   7.8074784E+04
0000000E+03        9.2860000E-01   2.5690000E+02   6.8466555E+04
0000000E+03        8.3540000E-01   2.4960000E+02   5.9844569E+04

If we want to truncate this data to altitude $z \ in[1000\,\textrm{m}, 3000\,\textrm{m}]$ we specify:

{"gas":
  {
    "law": "table",
    "file": "ENMSIS_Atmosphere.dat",
    "molar mass": 28.97,
    "height": 0,
    "temperature": 2,
    "pressure": 3,
    "density": 1,
    "min height": 1000,
    "max height": 3000,
    "res height": 10
  }
}

Neutral species background

Neutral species are included by an array neutral species in the gas object. Each neutral species must have the fields

name Species name
molar fraction Molar fraction of the species.

If the molar fractions do not add up to one, they will be normalized.

Warning

Neutral species are not tracked on the mesh. They are simply stored as functions that allow us to obtain the (spatially varying) density, temperature, and pressure for each neutral species. If a neutral species needs to be tracked on the mesh (through e.g. a convection-diffusion-reaction solver) it must be defined as a plasma species. See Plasma species.

For example, a standard nitrogen-oxygen atmosphere will look like:

  {"gas":
    {
      "law": "ideal",
      "temperature": 300,
      "pressure": 1,
      "plot": true,
      "neutral species":
      [
        {
          "name": "O2",
          "molar_fraction": 0.2
        },
        {
          "name": "N2",
          "molar_fraction": 0.8
        }
    ]
}

Plasma species

The list of plasma species is included by an array plasma species. Each entry must have the entries

name (string) For identifying the species name.
Z (integer) Species charge number.
mobile (true/false) Mobile species or not.
diffusive (true/false) Diffusive species or not.

Optionally, the field initial data, can be included for providing initial data to the species Details are discussed further below.

For example, a minimum version would look like

{"plasma species":
  [
    {"name": "N2+", "Z":  1, "mobile": false, "diffusive": false},
    {"name": "O2+", "Z":  1, "mobile": false, "diffusive": false},
    {"name": "O2-", "Z": -1, "mobile": false, "diffusive": false}
  ]
}

Initial data

Initial data can be provided with

Function based densities.
Computational particles (deposited using a nearest-grid-point scheme).

Density functions

To provide initial data one include initial data for each species. Currently, the following fields are supported:

uniform For specifying a uniform background density. Simply the field uniform and a density (in units of $m^{-3}$)
gauss2 for specifying Gaussian seeds $n = n_0\exp\left(-\frac{\left(\mathbf{x}-\mathbf{x_0}\right)^2}{2R^2}\right)$. gauss2 is an array where each array entry must contain
- radius, for specifying the radius $R$:
- amplitude, for specifying the amplitude $n_0$.
- position, for specifying the seed position $\mathbf{x}$.
The position must be a 2D/3D array.
gauss2 for specifying Gaussian seeds $n = n_0\exp\left(-\frac{\left(\mathbf{x}-\mathbf{x_0}\right)^4}{2R^4}\right)$. gauss4 is an array where each array entry must contain
- radius, for specifying the radius $R$:
- amplitude, for specifying the amplitude $n_0$.
- position, for specifying the seed position $\mathbf{x}$.
The position must be a 2D/3D array.
height profile For specifying a height profile along $y$ in 2D, and $z$ in 3D. To include it, prepare an ASCII files with at least two columns. The height (in meters) must be specified in one column and the density (in units of $m^{-3}$) in another. Internally, this data is stored in a lookup table (see LookupTable1D). Required fields are
- file , for specifying the file.
- height, for specifying the column that stores the height.
- density, for specifying the column that stores the density.
- min height, for trimming data to a minimum height.
- max height, for trimming data to a maximum height.
- res height, for specifying the resolution height in the chombo-discharge lookup tables.
In addition, height and density columns can be scaled in the internal tables by including
- scale height for scaling the height data.
- scale density for scaling the density data.

Note

When multiple initial data fields are specified, chombo-discharge takes the superposition of all of them.

Initial particles

Initial particles can be included with the initial particles field. The current implementation supports

uniform For drawing initial particles randomly distributed inside a box. The user must specify the two corners lo corner and hi corner that indicate the spatial extents of the box, and the number of computational particles to draw. The weight is specified by a field weight. For example:

{"plasma species":
  [
    {
      "name": "e",
      "Z":  -1,
      "mobile": true,
      "diffusive": true,
      "initial particles": {
        "uniform": {
          "lo corner": [0,0,0],
          "hi corner": [1,1,1],
          "number": 100,
          "weight": 1.0
        }
      }
    }
  ]
}

sphere For drawing initial particles randomly distributed inside a sphere. Mandatory fields are

center for specifying the sphere center.
radius for specifying the sphere radius.
number for the number of computational particles.
weight for the initial particle weight.

{"plasma species":
  [
    {
      "name": "e",
      "Z":  -1,
      "mobile": true,
      "diffusive": true,
      "initial particles": {
        "sphere": {
          "center": [0,0,0],
          "radius": 1.0,
          "number": 100,
          "weight": 1.0
        }
      }
    }
  ]
}

copy For using an already initialized particle distribution. The only mandatory fields is copy, e.g.

{"plasma species":
[
  {
    "name": "e",
    "Z":  -1,
    "mobile": true,
    "diffusive": true,
    "initial particles": {
      "sphere": {
        "center": [0,0,0],
        "radius": 1.0,
        "number": 100,
        "weight": 1.0
      }
    }
  },
  {
    "name": "O2+",
    "Z": 1,
    "mobile": true,
    "diffusive": true,
    "initial particles": {
       "copy": "e"
    }
  }
]}

This will copy the particles from the species e to the species O2+.

Warning

The species one copies from must be defined before the species one copies to.

Complex example

For example, a species with complex initial data that combines density functions with initial particles can look like:

{"plasma species":
  [
    {
      "name": "N2+",
      "Z":  1,
      "mobile": false,
      "diffusive": false,
      "initial data": {
        "uniform": 1E10,
        "gauss2" :
          [
            {
               "radius": 100E-6,
               "amplitude": 1E18,
               "position": [0,0,0]
            },
            {
               "radius": 200E-6,
               "amplitude": 2E18,
               "position": [1,0,0]
            }
          ],
         "gauss4":
           [
             {
               "radius": 300E-6,
               "amplitude": 3E18,
               "position": [0,1,0]
             },
             {
               "radius": 400E-6,
               "amplitude": 4E18,
               "position": [0,0,1]
             }
           ],
         "height profile": {
           "file": "MyHeightProfile.dat",
           "height": 0,
           "density": 1,
           "min height": 0,
           "max height": 100000,
           "res height": 10,
           "scale height": 100,
           "scale density": 1E6
         }
      },
      "initial particles": {
        "sphere": {
          "center": [0,0,0],
          "radius": 1.0,
          "number": 100,
          "weight": 1.0
        }
      }
    }
  ]
}

Mobilities

If a species is specified as mobile, the mobility is set from a field mobility, and the field lookup is used to specify the method for computing it. Currently supported are:

Constant mobility.
Function-based mobility, i.e. $\mu = \mu(E,N)$.
Tabulated mobility, i.e. $\mu = \mu(E,N)$.

The cases are discussed below.

Constant mobility

Setting lookup to constant lets the user set a constant mobility. If setting a constant mobility, the field value is also required. For example:

{"plasma species":
  [
    {"name": "e", "Z":  -1, "mobile": true, "diffusive": false,
     "mobility": {
       "lookup" : "constant",
       "value": 0.05,
      }
    }
  ]
}

Function-based mobility

Setting lookup to function E/N lets the user set the mobility as a function of the reduced electric field. When setting a function-based mobility, the field function is also required.

Supported functions are:

ABC, in which case the mobility is computed as

\[\mu(E) = A \frac{E^B}{N^C}.\]

The fields A, B, and C must also be specified. For example:

{"plasma species":
  [
    {"name": "e", "Z":  -1, "mobile": true, "diffusive": false,
     "mobility": {
       "lookup" : "function E/N",
       "function": "ABC",
       "A": 1,
       "B": 1,
       "C": 1
      }
    }
  ]
}

Tabulated mobility

Specifying lookup to table E/N lets the user set the mobility from a tabulated value of the reduced electric field. BOLSIG-like files can be parsed by specifying the header which contains the tabulated data, and the columns that identify the reduced electric field and mobilities. This data is then stored in a lookup table, see LookupTable1D.

For example:

{"plasma species":
  [
    {"name": "e", "Z":  -1, "mobile": true, "diffusive": false,
     "mobility": {
       "lookup" : "table E/N",
       "file": "transport_file.txt",
       "header": "# Electron mobility (E/N, mu*N)",
       "E/N ": 0,
       "mu*N": 1,
       "min E/N": 10,
       "max E/N": 1000,
       "points": 100,
       "spacing": "exponential",
       "dump": "MyMobilityTable.dat"
      }
    }
  ]
}

In the above, the fields have the following meaning:

file The file where the data is found. The data must be stored in rows and columns.
header, the contents of the line preceding the table data.
E/N, the column that contains $E/N$.
mu*N, the column that contains $\mu\cdot E$.
min E/N, for trimming the data range.
max E/N, for trimming the data range.
points, for specifying the number of points in the lookup table.
spacing, for specifying how to regularize the table.
dump, an optional argument (useful for debugging) which will write the table to file.

Note that the input file does not need regularly spaced or sorted data. For performance reasons, the tables are always resampled, see LookupTable1D.

Diffusion coefficients

Setting the diffusion coefficient is done exactly in the same was as the mobility. If a species is diffusive, one must include the field diffusion as well as lookup. For example, the JSON input for specifying a tabulated diffusion coefficient is done by

{"plasma species":
  [
    {"name": "e", "Z":  -1, "mobile": false, "true": false,
     "diffusion": {
       "lookup" : "table E/N",
       "file": "transport_file.txt",
       "header": "# Electron diffusion coefficient (E/N, D*N)",
       "E/N ": 0,
       "D*N": 1,
       "min E/N": 10,
       "max E/N": 1000,
       "points": 1000,
       "spacing": "exponential"
      }
    }
  ]
}

Temperatures

Plasma species temperatures can set by including a field temperature for the plasma species.

Warning

If the temperature field is omitted, the species temperature will be set to the gas temperature.

Constant temperature

To set a constant temperature, include the field temperature and set lookup to constant and specify the temperature through the field value as follows:

{"plasma species":
  [
    {
      "name": "O2",
      "Z":  0,
      "mobile": false,
      "true": false,
      "temperature": {
        "lookup": "constant",
        "value": 300
      }
     }
  ]
}

Tabulated temperature

To include a tabulated temperature $T = T(E,N)$, set lookup to table E/N. The temperature is then computed as

\[T = \frac{2 \epsilon}{3k_{\textrm{B}}},\]

where $\epsilon$ is the energy and $k_{\textrm{B}}$ is the Boltzmann constant.

The following fields are required:

file for specifying which file the temperature is stored.
header for specifying where in the file the temperature is stored.
E/N for specifying in which column we find $E/N$.
eV for specifying in which column we find the species energy (in units of electron volts).
min E/N for trimming the data range.
max E/N for trimming the data range.
points for setting the number of points in the lookup table.
spacing for setting the grid point spacing type.
dump for writing the final table to file.

For a further explanation to these fields, see Mobilities.

A complete example is:

{"plasma species":
  [
    {
      "name": "e",
      "Z":  -1,
      "mobile": true,
      "true": true,
      "temperature": {
        "lookup": "table E/N",
        "file": "transport_data.txt",
        "header": "# Electron mean energy (E/N, eV)",
        "E/N": 0,
        "eV": 1,
        "min E/N": 10,
        "max E/N": 1000,
        "points": 1000,
        "spacing": "exponential",
        "dump": "MyTemperatureTable.dat"
      }
     }
  ]
}

Photon species

As for the plasma species, photon species (for including radiative transfer) are included by an array photon species. For each species, the required fields are

name For setting the species name.
kappa For specifying the absorption coefficient.

Currently, kappa can be either

constant Which lets the user set a constant absorption coefficient.
helmholtz Computes the absorption coefficient as

\[\kappa = \frac{p_X\lambda}{\sqrt{3}}\]

where $\lambda$ is a specified input parameter and $p_X$ is the partial pressure of some species $X$.
stochastic A which samples a random absorption coefficient as

\[\kappa = K_1 \left(\frac{K_2}{K_1}\right)^{\frac{f-f1}{f2-f1}}.\]

Here, $f_1$ and $f_2$ are frequency ranges, $K_1$ and $K_2$ are absorption coefficients, and $f$ is a stochastically sampled frequency. Note that this method is only sensible when using discrete photons.

Constant absorption coefficients

When specifying a constant absorption coefficient, one must include a field value as well. For example:

{"photon species":
   [
     {
       "name": "UVPhoton",
       "kappa": "constant",
       "value": 1E4
     }
   ]
}

Helmholtz absorption coefficients

The interface for the Helmholtz-based absorption coefficients are inspired by Bourdon et al. [2007] approach for computing photoionization. This method only makes sense if doing a Helmholtz-based reconstruction of the photoionization profile as a relation:

\[\left[\nabla^2 - \left(p_{\textrm{O}_2} \lambda\right)^2\right]S_\gamma = -\left(A p_{\textrm{O}_2}^2\frac{p_q}{p + p_q}\xi\nu\right)S_i,\]

where

$S_\gamma$ is the number of photoionization events per unit volume and time.
$A$ is a model coefficient.
$\frac{p_q}{p + p_q}$ is a quenching factor.
$\xi$ is a photoionization efficiency.
$\nu$ is a relative excitation efficiency.
$S_i$ is the electron impact ionization source term.

Since the radiative transfer solver is based on the Eddington approximation, the Helmholtz reconstruction can be written as

\[\kappa \Psi - \nabla\cdot\left(\frac{1}{3\kappa}\nabla \Psi\right) = \frac{\eta}{c}\]

where the absorption coefficient is set as

\[\kappa(\mathbf{x}) = \frac{p_{\textrm{O}_2}\lambda}{\sqrt{3}}.\]

The photogeneration source term is still

\[\eta = \frac{p_q}{p + p_q}\xi\nu S_i,\]

but the photoionization term is

\[S_\gamma = \frac{c A p_{\textrm{O}_2}}{\sqrt{3}\lambda}\Psi.\]

Note that the photoionization term is, in principle, not an Eddington approximation. Rather, the Eddington-like equations occur here through an approximation of the exact integral solution to the radiative transfer problem. In the pure Eddington approximation, on the other hand, $\Psi$ represents the total number of ionizing photons per unit volume, and we would have $S_\gamma = \frac{\Psi}{\Delta t}$ where $\Delta t$ is the time step.

When specifying the kappa field as helmholtz, the absorption coefficient is computed as

\[\kappa(\mathbf{x}) = \frac{p_X\left(\mathbf{x}\right)\lambda}{\sqrt{3}}\]

where $p_X$ is the partial pressure of a species $X$ and $\lambda$ is the same input parameter as in the Helmholtz reconstruction. These are specified through fields neutral and lambda as follows:

{"photon species":
   [
     {
       "name": "UVPhoton",
       "kappa": "helmholtz",
       "lambda": 0.0415,
       "neutral": "O2"
     }
   ]
}

This input will set $\kappa\left(\mathbf{x}\right) = \frac{p_{\textrm{O}_2}\left(\mathbf{x}\right)\lambda}{\sqrt{3}}$.

Note

The source term $\eta$ is specified when specifying the plasma reactions, see Plasma reactions.

Stochastic sampling

Setting the kappa field to stochastic A will stochastically sample the absorption length from

\[\kappa = K_1 \left(\frac{K_2}{K_1}\right)^{\frac{f-f1}{f2-f1}}.\]

where $K_1 = p_X\chi_{\textrm{min}}$, $K_1 = p_X\chi_{\textrm{max}}$, and $f_1$ and $f_2$ are frequency ranges. Like above, $p_X$ is the partial pressure of some species $X$. Note that all input parameters are given in SI units.

Stochastic sampling of the absorption length only makes sense when using discrete photons – this particular method is inspired by the method in Chanrion and Neubert [2008]. For example:

{"photon species":
   [
     {
       "name": "UVPhoton",
       "kappa": "stochastic A",
       "neutral": "O2",
       "f1":   2.925E15,
       "f2":   3.059E15,
       "chi min": 2.625E-2,
       "chi max": 1.5
     }
   ]
}

Plasma reactions

Plasma reactions are reactions between charged and neutral species and are written in the form

\[A + B + \ldots \rightarrow C + D + \ldots.\]

Importantly, the left hand side of the reaction can only consist of charged or neutral species. It is not permitted to put a photon species on the left hand side of these reactions; photo-ionization is handled separately by another set of reaction types (see Photo-reactions). However, photon species can appear on the left hand side of the equation.

When specifying reactions in this form, the reaction rate is computed as

\[R = k n_A n_B\ldots\]

When computing the source term for some species $X$, we subtract $R$ for each time $X$ appears on the left hand side of the reaction and add $R$ for each time $X$ appears on the right-hand side of the reaction.

Specifying reactions

Reactions of the above type are handled by a JSON array plasma reactions, with required fields:

reaction (string) containing the reaction process.
lookup (string) for determining how to compute the reaction rate.

{"plasma reactions":
  [
    {
      "reaction": "e + O2 -> e + e + O2+",
      "lookup": "constant",
      "rate": 1E-30
    }
  ]
}

This adds a reaction $\textrm{e} + \textrm{O}_2 \rightarrow \textrm{e} + \textrm{e} + \textrm{O}_2^+$ to the reaction set. We compute

\[R = kn_{\textrm{e}}n_{\textrm{O}_2^+}\]

and set

\[S_{\textrm{e}} = S_{\textrm{O}_2^+} = R.\]

Some caveats when setting the reaction string are:

Whitespace are separators. For example, O2+e will be interpreted as a species with string identifier O2+e, but O2 + e will interpreted as a reaction between O2 and e.
The reaction string must contain a left and right hand side separated by ->. An error will be thrown if this symbol can not be found.
The left-hand must consist only of neutral or plasma species. If the left-hand side consists of species that are not neutral or plasma species, an error will be thrown.
The right-hand side can consist of either neutral, plasma species, or photon species. Otherwise, an error will be thrown.
The reaction string will be checked for charge conservation.

Note that if a reaction involves a right-hand side that is not otherwise tracked, the user should omit the species from the right-hand side altogether. For example, if we have a model which tracks the species $e$ and $\textrm{O}_2^+$ but we want to include the dissociative recombination reaction $e + \textrm{O}_2^+ \rightarrow O + O$, this reaction should be added to the reaction with an empty right-hand side:

{"plasma reactions":
  [
    {
      "reaction": "e + O2 -> e + e + O2+",
      "lookup": "constant",
      "rate": 1E-30
    },
    {
      "reaction": "e + O2+ -> ",
      "lookup": "constant",
      "rate": 1E-30
    }
  ]
}

Wildcards

Reaction specifiers may include the wildcard @ which is a placeholder for another species. The wildcards must be specified by including a JSON array @ of the species that the wildcard is replaced by. For example:

{"plasma reactions":
  [
    {
      "reaction": "N2+ + N2 + @ -> N4+ + @",
      "@": ["N2", "O2"],
      "lookup": "constant",
      "rate": 1E-30
    }
  ]
}

The above code will add two reactions to the reaction set: $N_2 + N_2 + N_2 \rightarrow N_4^+ + N_2$ and $N_2 + N_2 + \textrm{O}_2 \rightarrow N_4^+ + \textrm{O}_2$. It is not possible to set different reaction rates for the two reactions.

Specifying reaction rates

Constant reaction rates

To set a constant reaction rate for a reaction, set the field lookup to "constant" and specify the rate. For example:

{"plasma reactions":
  [
    {
      "reaction": "e + O2 -> e + e + O2+",
      "lookup": "constant",
      "rate": 1E-30
    }
  ]
}

Single-temperature rates

functionT A To set a rate dependent on a single species temperature in the form $k(T) = c_1T^{c_2}$, set lookup to functionT A. The user must specify the species from which we compute the temperature $T$ by including a field T. The constants $c_1$ and $c_2$ must also be included.

For example, in order to add a reaction $e + \textrm{O}_2 \rightarrow \varnothing$ with rate $k = 1.138\times 10^{-11}T_{\textrm{e}}^{-0.7}$ we can add the following:
{"plasma reactions": [ { "reaction": "e + M+ ->", "lookup": "functionT A", "T": "e", "c1": 1.138 "c2": -0.7 } ] }

Two-temperature rates

functionT1T2 A To set a rate dependent on two species temperature in the form $k(T_1, T_2) = c_1\left(T_1/T_2\right)^{c_2}$, set lookup to functionT1T2 A. The user must specify which temperatures are involved by specifying the fields T1, T2, as well as the constants through fields c1 and c2. For example, to include the reaction $e + \textrm{O}_2 + \textrm{O}_2 \rightarrow \textrm{O}_2^- + O2$ in the set, with this reaction having a rate

\[k = 2.4\times 10^{-41}\left(\frac{T_{\textrm{O}_2}}{T_e}\right),\]

we add the following:
```
{"plasma reactions":
  [
    {
      "reaction": "e + O2 + O2 -> O2- + O2",
      "lookup": "functionT1T2 A",
      "T1": "O2",
      "T2": "e",
      "c1": 2.41E-41,
      "c2": 1
    }
  ]
}
```

Townsend ionization and attachment

To set standard Townsend ionization and attachment reactions, set lookup to alpha*v and eta*v, respectively. This will compute the rate constant $k = \alpha \left|\mathbf{v}\right|$ where $\mathbf{v}$ is the drift velocity of some species. To specify the species one includes the field species.

For example, to include the reactions $e \rightarrow e + e + M^+$ and $e \rightarrow M^-$ one can specify the reactions as

{"plasma reactions":
  [
    {
      "reaction": "e -> e + e + M+",
      "lookup": "alpha*v",
      "species": "e"
    },
    {
      "reaction": "e -> M-",
      "lookup": "eta*v",
      "species": "e"
    }
  ]
}

Tabulated rates

To set a tabulated rate with $k = k(E,N)$, set the field lookup to table E/N and specify the file, header, and data format to be used. For example:

{"plasma reactions":
  [
    {
      "reaction": "e + O2 -> e + e + O2+",
      "lookup": "table E/N",
      "file": "transport_file.txt",
      "header": "# O2 ionization (E/N, rate/N)",
      "E/N ": 0,
      "rate/N": 1,
      "min E/N": 10,
      "max E/N": 1000,
      "spacing": "exponential",
      "points": 1000,
      "plot": true,
      "dump": "O2_ionization.dat"
    }
  ]
}

The file field specifies which field to read the reaction rate from, while header indicates where in the file the reaction rate is found. The file parser will read the files below the header line until it reaches an empty line. The fields E/N and rate/N indicate the columns where the reduced electric field and reaction rates are stored.

The final fields min E/N, max E/N, and points are formatting fields that trim the range of the data input and organizes the data along a table with points entries. As with the mobilities (see Mobilities), the spacing argument determines whether or not the internal interpolation table uses uniform or exponential grid point spacing. Finally, the dump argument will tell chombo-discharge to dump the table to file, which is useful for debugging or quality assurance of the tabulated data.

Modifying reactions

Collisional quenching

To quench a reaction, include a field qenching_pressure and specify the quenching pressure (in atmospheres). When computing reaction rates, the rate for the reaction will be modified as

\[k \rightarrow k\frac{p_q}{p_q + p}\]

where $p^q$ is the quenching pressure and $p = p(\mathbf{x})$ is the gas pressure.

Important

The quenching pressure should be specified in Pascal.

For example:

{"plasma reactions":
  [
    {
      "reaction": "e + N2 -> e + N2 + Y",
      "lookup": "table E/N",
      "file": "transport_file.txt",
      "header": "# N2 ionization (E/N, rate/N)",
      "E/N ": 0,
      "rate/N": 1,
      "min E/N": 10,
      "max E/N": 1000,
      "points": 1000,
      "spacing": "exponential",
      "quenching pressure": 4000
    }
  ]
}

Reaction efficiencies

To modify a reaction efficiency, include a field efficiency and specify it. This will modify the reaction rate as

\[k \rightarrow \nu k\]

where $\nu$ is the reaction efficiency. For example:

{"plasma reactions":
  [
    {
      "reaction": "e + N2 -> e + N2 + Y",
      "lookup": "table E/N",
      "file": "transport_file.txt",
      "header": "# N2 ionization (E/N, rate/N)",
      "E/N ": 0,
      "rate/N": 1,
      "min E/N": 10,
      "max E/N": 1000,
      "points": 1000,
      "spacing": "exponential",
      "efficiency": 0.6
    }
  ]
}

Scaling reactions

Reactions can be scaled by including a scale argument to the reaction. This works exactly like the efficiency field outlined above.

Energy correction

Occasionally, it can be necessary to incorporate an energy correction to models, accounting e.g. for electron energy loss near strong gradients. The JSON interface supports the correction in Soloviev and Krivtsov [2009]. To use it, include an (optional) field soloviev and specify correction and species. For example:

{"plasma reactions":
  [
    {
      "reaction": "e + N2 -> e + N2 + Y",
      "lookup": "table E/N",
      "file": "transport_file.txt",
      "header": "# N2 ionization (E/N, rate/N)",
      "E/N ": 0,
      "rate/N": 1,
      "min E/N": 10,
      "max E/N": 1000,
      "points": 1000,
      "spacing": "exponential",
      "efficiency": 0.6,
      "soloviev": {
        "correction": true,
        "species": "e"
      }
    }
  ]
}

When this energy correction is enabled, the rate coefficient is modified as

\[k \rightarrow k\left(1 + \frac{\mathbf{E}\cdot D_s\nabla n_s}{\mu_s n_s E^2}\right),\]

where $s$ is the species specified in the soloviev field, $n_s$ is the density and $D_s$ and $\mu_s$ are diffusion and mobility coefficients. We point out that the correction factor is restricted such that the reaction rate is always non-negative. Note that this correction makes sense when rates are dependent only on the electric field, see Soloviev and Krivtsov [2009].

Note

When using the energy correction, the specifies species must be both mobile and diffusive.

Plotting reactions

It is possible to have CdrPlasmaJSON include the reaction rates in the HDF5 output files by including a field plot as follows:

{"plasma reactions":
  [
    {
      "reaction": "e + O2 -> e + e + O2+",
      "plot": true,
      "lookup": "constant",
      "rate": 1E-30,
    }
  ]
}

Plotting the reaction rate can be useful for debugging or analysis. Note that it is, by extension, also possible to add useful data to the I/O files from reactions that otherwise do not contribute to the discharge evolution. For example, if we know the rate $k$ for excitation of nitrogen to a specific excited state, but do not otherwise care about tracking the excited state, we can add the reaction as follows:

{"plasma reactions":
  [
    {
      "reaction": "e + N2 -> e + N2",
      "plot": true,
      "lookup": "constant",
      "rate": 1E-30,
    }
  ]
}

This reaction is a dud in terms of the discharge evolution (the left and right hand sides are the same), but it can be useful for plotting the excitation rate.

Note

This functionality should be used with care because each reaction increases the I/O load.

Warnings and caveats

Note that the JSON interface always computes reactions as if they were specified by the deterministic reaction rate equation

\[\partial_t n_i = \sum_r k_r n_j n_k n_l\ldots,\]

where the fluid source term for any reaction $r$ is $S_r = k_r n_j n_k n_l\ldots$. Caution should always be exercised when defining a reaction set.

Higher-order reactions

Usually, many rate coefficients depend on the output of other software (e.g., BOLSIG+) and the scaling of rate coefficients is not immediately obvious. This is particularly the case for three-body reactions with BOLSIG+ that may require scaling before running the Boltzmann solver (by scaling the input cross sections), or after running the Boltzmann solver, in which case the rate coefficients themselves might require scaling. In any case the user should investigate the cross-section file that BOLSIG+ uses, and figure out the required scaling.

Important

For two-body reactions, e.g. $A + B\rightarrow \varnothing$ the rate coefficient must be specified in units of $\textrm{m}^3\textrm{s}^{-1}$, while for three-body reactions $A + B + C\rightarrow\varnothing$ the rate coefficient must have units of $\textrm{m}^6\textrm{s}^{-1}$.

For three-body reactions the units given by BOLSIG+ in the output file may or may not be incorrect (depending on whether or not the user scaled the cross sections).

Townsend coefficients

Townsend coefficients are not fundamentally required for specifying the reactions, but as with the higher-order reactions some of the output rates for three-body reactions might be inconsistently represented in the BOLSIG+ output files. For example, some care might be required when using the Townsend attachment coefficient for air when the reaction $\textrm{e} + \textrm{O}_2 + \textrm{O}_2\rightarrow \textrm{O}_2^- + \textrm{O}_2$ is included because the rate constant might require proper scaling after running the Boltzmann solver, but this scaling is invisible to the BOLSIG+’s calculation of the attachment coefficient $\eta/N$.

Warning

The JSON interface does not guard against inconsistencies in the user-provided chemistry, and provision of inconsistent $\eta/N$ and attachment reaction rates are quite possible.

Photo-reactions

Photo-reactions are reactions between charged/neutral and photons in the form

\[A + B + \gamma + \ldots \rightarrow C + D + \ldots.\]

where species $A, B, \ldots$ are charged and neutral species and $\gamma$ is a photon. The left hand side can contain only one photon species, and the right-hand side can not contain a photon species. In other words, two-photon absorption is not supported, and photons that are absorbed on the mesh cannot become new photons. This is not a fundamental limitation, but a restriction imposed by the JSON interface.

Specifying reactions

Reactions of the above type are handled by a JSON array photo reactions, with required fields:

reaction (string) containing the reaction process.
lookup (string) for determining how to compute the reaction rate.

For example:

{"photo reactions":
  [
    {"reaction": "Y + O2 -> e + O2+"},
  ]
}

The rules for specifying reaction strings are the same as for the plasma reactions, see Plasma reactions. Wildcards also apply, see Wildcards.

Default behavior

Since the radiative transfer solvers solve for the number of ionizing photons, the CDR solver source terms are incremented by

\[S \rightarrow S + \frac{\Psi}{\Delta t}.\]

where $\Psi$ is the number of ionizing photons per unit volume (i.e., the solution $\Psi$).

Helmholtz reconstruction

When performing a Helmholtz reconstruction the photoionization source term is

\[S = \frac{c A p_{\textrm{O}_2}}{\sqrt{3}\lambda}\Psi.\]

To modify the source term for consistency with Helmholtz reconstruction specify the field helmholtz with variables

A. the $A$ coefficient.
lambda. the $\lambda$ coefficient. This value will also be specified in the photon species, but it is not retrieved automatically.
neutral. The neutral species for which we obtain the partial pressure.

For example:

{"photo reactions":
  [
    {
      "reaction": "Y + O2 -> e + O2+",
      "helmholtz": {
        "A": 1.1E-4,
        "lambda": 0.0415,
        "neutral": "O2"
      }
    }
  ]
}

Scaling reactions

Photo-reactions can be scaled by including a scale argument. For example, to completely turn off the photoreaction above:

{"photo reactions":
  [
    {
      "reaction": "Y + O2 -> e + O2+",
      "helmholtz": {
        "A": 1.1E-4,
        "lambda": 0.0415,
        "neutral": "O2"
      },
      "scale": 0.0
    }
  ]
}

EB boundary conditions

Boundary conditions on the embedded boundary are included by the fields

electrode reactions, for specifying secondary emission on electrodes.
dielectric reactions, for specifying secondary emission on dielectrics.

To include secondary emission, the user must specify a reaction string in the form $A \rightarrow B$, and also include an emission rate. Currently, we only support constant emission rates (i.e., secondary emission coefficients). This is likely to change in the future.

The following points furthermore apply:

By default, standard outflow boundary conditions. When electrode reactions or dielectric reactions are specified, the user only controls the inflow back into the domain.
Wildcards can appear on the left hand side of the reaction.
If one specifies $A + B \rightarrow C$ for a surface reaction, this is the same as specifying two reactions $A \rightarrow C$ and $B\rightarrow C$. The same emission coefficient will be used for both reactions.
Both photon species and plasma species can appear on the left hand side of the reaction.
Photon species can not appear on the right-hand side of the reaction; we do not include surface sources for photoionization.
To scale reactions, include a modifier scale.

For example, the following specification will set secondary emission efficiencies to $10^{-3}$:

{"electrode reactions":
  [
    { "reaction": "@ -> e",
      "@": ["N2+", "O2+", "N4+", "O4+", "O2+N2"],
      "lookup": "constant",
      "value": 1E-4
    }
  ],
 "dielectric reactions":
  [
    { "reaction": "@ -> e",
      "@": ["N2+", "O2+", "N4+", "O4+", "O2+N2"],
      "lookup": "constant",
      "value": 1E-3
    }
  ]

Setting up a new problem

New problems that use the CdrPlasma physics model are best set up by using the Python tools provided with the module. A full description is available in the README.md file contained in the folder:

# Physics/CdrPlasma
This physics module solves the minimal plasma model

The folders contains

* **PlasmaModels** Implementations of various plasma models. 
* **TimeSteppers** Implementations of the various time integration algorithms.
* **CellTaggers** Implementations of various cell taggers.
* **Data** Various data (e.g. transport data). 

## Setting up a new application
To set up a new problem, use the Python script. For example:

```shell
python setup.py -base_dir=/home/foo/MyApplications -app_name=MyPlasmaModel -geometry=Vessel
```

To install within chombo-discharge:

```shell
python setup.py -base_dir=$DISCHARGE_HOME/MyApplications -app_name=MyPlasmaModel -geometry=Vessel
```

The application will then be installed to $DISCHARGE_HOME/MyApplications/MyPlasmaApplication
The user will need to modify the geometry and set the initial conditions through the inputs file. 

## Modifying the application
Users are free to modify this application, e.g. adding new initial conditions and flow fields.

To see the list of available options type

cd $DISCHARGE_HOME/Physics/CdrPlasma
python setup.py --help

CDR plasma model

Time discretizations

Godunov operator splitting

Specifying transport algorithm

Specifying diffusion

Time step limitations

Solver configuration

Spectral deferred corrections

SDC basics

Time step limitations

CdrPlasmaPhysics

Overview

Complete API

JSON interface

Usage

Specifying input file

Discrete photons

Gas law and neutral background

General functionality

Ideal gas

Troposphere

Tabulated

Neutral species background

Plasma species

Initial data

Density functions

Initial particles

Complex example

Mobilities

Diffusion coefficients

Temperatures

Photon species

Plasma reactions

Specifying reactions

Wildcards

Specifying reaction rates

Constant reaction rates

Single-temperature rates

Two-temperature rates

Townsend ionization and attachment

Tabulated rates

Modifying reactions

Collisional quenching

Reaction efficiencies

Scaling reactions

Energy correction

Plotting reactions

Warnings and caveats

Higher-order reactions

Townsend coefficients

Photo-reactions

Specifying reactions

Default behavior

Helmholtz reconstruction

Scaling reactions

EB boundary conditions

Adaptive mesh refinement

Setting up a new problem