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 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 intantiated. 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:

/*!
  @brief Abstract base class for specifying plasma kinetics. This is the base class used by CdrPlasmaStepper when advancing the minimal plasma model. 
*/
class CdrPlasmaPhysics
{
public:
  /*!
    @brief Default constructor. Does nothing. 
  */
  CdrPlasmaPhysics()
  {}

  /*!
    @brief Base destructor. Does nothing.
  */
  virtual ~CdrPlasmaPhysics()
  {}

  /*!
    @brief Parse run-time options
  */
  virtual void
  parseRuntimeOptions() {

  };

  /*!
    @brief Get number of plot variables for this physics class. 
    @details This is used by CdrPlasmaStepper for pre-allocating data that will be put in a plot file. When overriding
    this method then this routine should return the number of plot variables that will be plotted. The returned number
    should thus be the same as the vectors that are returned from getPlotVariableNames and getPlotVariables.
  */
  virtual int
  getNumberOfPlotVariables() const
  {
    return 0;
  }

  /*!
    @brief Get plot variable names. The
    @details This function will return the plot variable names that CdrPlasmaPhysics will plot. The length of the returned
    vector must be the same as the returned value from getNumberOfPlotVariables. Note that the names/positions between this 
    routine and getPlotVariables should be consistent.
  */
  virtual Vector<std::string>
  getPlotVariableNames() const
  {
    return Vector<std::string>(this->getNumberOfPlotVariables(), std::string("empty data"));
  }

  /*!
    @brief Provide plot variables. This is used by CdrPlasmaStepper when writing plot files. 
    @details The length of the returned vector should be the same as the returned value from getNumberOfPlotVariables. The names
    for the returned variables are given by the returned vector in getPlotVariableNames().
    @param[in]  a_cdrDensities  Grid-based density for particle species.
    @param[in]  a_cdrGradients  Grid-based gradients for particle species.
    @param[in]  a_rteDensities  Grid-based densities for photons.
    @param[in]  a_E             Electric field.
    @param[in]  a_pos           Position in space.
    @param[in]  a_dx            Grid resolution. 
    @param[in]  a_dt            Advanced time.
    @param[in]  a_time          Current time.
    @param[in]  a_kappa         Grid cell unit volume. 
  */
  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 alpha. Should return Townsend ionization coefficient. 
    @details This function is mostly used for the cell tagging classes
    @param[in] a_E        Electric field.
    @param[in] a_position Physical coordinates
  */
  virtual Real
  computeAlpha(const Real a_E, const RealVect a_position) const = 0;

  /*!
    @brief Compute eta. Should return Townsend attachment coefficient. 
    @details This function is mostly used for the cell tagging classes
    @param[in] a_E        Electric field.
    @param[in] a_position Physical coordinates
  */
  virtual Real
  computeEta(const Real a_E, const RealVect a_position) const = 0;

  /*!
    @brief Routine intended for advancing a reaction network over a time a_dt. 
    @details This routine assumes that the subsequent advance is in the form phi^(k+1) = phi^k + S*a_dt. Thus, this routine exists such that users can EITHER
    fill a_cdrSources and a_rteSources directly with an explicit rule, OR they can perform a fully implicit advance within this routine and set S from that. 
    @param[out] a_cdrSources    Source terms for CDR equations.
    @param[out] a_rteSources    Source terms for RTE equations.
    @param[in]  a_cdrDensities  Grid-based density for particle species.
    @param[in]  a_cdrGradients  Grid-based gradients for particle species.
    @param[in]  a_rteDensities  Grid-based densities for photons.
    @param[in]  a_E             Electric field.
    @param[in]  a_pos           Position in space.
    @param[in]  a_dx            Grid resolution. 
    @param[in]  a_dt            Advanced time.
    @param[in]  a_time          Current time.
    @param[in]  a_kappa         Grid cell unit volume. 
  */
  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 velocities for the CDR equations
    @param[in] a_time         Time
    @param[in] a_pos          Position
    @param[in] a_E            Electric field
    @param[in] a_cdrDensities CDR densities
    @return Returns the drift velocities for each CDR species. The vector ordering is the same 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 the CDR equations. 
    @param[in] a_time         Time
    @param[in] a_pos          Position
    @param[in] a_E            Electric field
    @param[in] a_cdrDensities CDR densities
    @return Returns the diffusion coefficients for each CDR species. The vector ordering is the same 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 fluxes on electrode-gas interfaces. This is used as a boundary condition in the CDR equations. 
    @param[in] a_time            Time
    @param[in] a_pos             Position
    @param[in] a_normal          Boundary normal vector. This points into the gas phase. 
    @param[in] a_E               Electric field
    @param[in] a_cdrVelocities   CDR velocities. Normal component only. 
    @param[in] a_cdrDensities    CDR densities. 
    @param[in] a_cdrGradients    Normal gradients of cdr densities
    @param[in] a_rteFluxes       RTE fluxes (normal component only)
    @param[in] a_extrapCdrFluxes Extrapolated fluxes from the gas side. 
    @return Returns the flux on an electrode interface cell. The vector ordering must be the same as m_cdrSpecies. 
    @note A positive flux is a flux INTO the domain. 
  */
  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 fluxes on dielectric-gas interfaces. This is used as a boundary condition in the CDR equations. 
    @param[in] a_time            Time
    @param[in] a_pos             Position
    @param[in] a_normal          Normal vector. This points into the gas phase. 
    @param[in] a_E Electric      field
    @param[in] a_cdrDensities    CDR densities (on the EB)
    @param[in] a_cdrVelocities   Normal component of CDR velocities (on the EB).
    @param[in] a_cdrGradients    Normal gradients of cdr densities
    @param[in] a_rteFluxes       RTE fluxes (normal component only)
    @param[in] a_extrapCdrFluxes Extrapolated fluxes from the gas side. 
    @return Returns the flux on a dielectric interface cell. The vector ordering must be the same as m_cdrSpecies. 
    @note A positive flux is a flux INTO the domain. 
  */
  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 fluxes through domain sides. This is used as a boundary condition in the CDR equations.
    @param[in] a_time            Time
    @param[in] a_pos             Position
    @param[in] a_dir             Direction (0 = x, 1=y etc)
    @param[in] a_side            Side (low or high side)
    @param[in] a_E               Electric field
    @param[in] a_cdrDensities    CDR densities. 
    @param[in] a_cdrVelocities   CDR velocities (normal component only). 
    @param[in] a_cdrGradients    CDR gradients (normal component only)
    @param[in] a_rteFluxes       RTE fluxes (normal component only)
    @param[in] a_extrapCdrFluxes Extrapolated fluxes from the gas side. 
    @note A positive flux is a flux INTO the domain. 
  */
  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 Set the initial surface charge
    @param[in] a_time Time
    @param[in] a_pos Position
  */
  virtual Real
  initialSigma(const Real a_time, const RealVect a_pos) const = 0;

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 incldue 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