Îto-KMC plasma model
Underlying model
Plasma transport
The Îto-KMC model uses an Îto solver for some of the species, i.e. the particles drift and diffuse according to
where \(\mathbf{X}\) is the particle position and \(\mathbf{V}\) and \(D\) is the particle drift velocity and diffusion coefficients, respectively. These are obtained by interpolating the macroscopic drift velocity and diffusion coefficients to the particle position. Further details regarding the Îto method are given in ItoSolver.
Not all species in the Îto-KMC model need to be defined by particle solvers, as some species can be tracked by more conventional advection-diffusion-reaction solvers, see CdrSolver for discretization details.
Photoionization
Photoionization in the Îto-KMC model is done using discrete photons. These are generated and advanced in every time step, and interact with the plasma through user-defined photo-reactions. The underlying solver is the discrete Monte Carlo photon transport solver, see Monte Carlo sampling.
Field interaction
The Îto-KMC model uses an electrostatic approximation where the field is obtained by solving the Poisson equation for the potential. See Electrostatic solver for further details.
Chemistry
Kinetic Monte Carlo (KMC) is used within grid cells for resolving the plasma chemistry. The algorithmic concepts are given in Kinetic Monte Carlo.
Algorithms
Time stepping
Reaction network
Particle management
0D chemistry
The user input interface to the Îto-KMC model consists of a zero-dimensional plasma kinetics interface called ItoKMCPhysics.
This interface consists of the following main functionalities:
User-defined species, and which solver types (Îto or CDR) to use when tracking them.
User-defined initial particles, reaction kinetics, and photoionization processes.
ItoKMCPhysics
The complete C++ interface specification is given below.
Because the interface is fairly extensive, chombo-discharge also supplies a JSON-based implementation called ItoKMCJSON (see JSON 0D chemistry interface) for defining these things through file input.
The difference between ItoKMCPhysics and its implementation ItoKMCJSON is that the JSON-based implementation class only implements a subset of potential features supported by ItoKMCPhysics.
Show/hide full 0D interface
/* chombo-discharge
* Copyright © 2021 SINTEF Energy Research.
* Please refer to Copyright.txt and LICENSE in the chombo-discharge root directory.
*/
/*!
@file CD_ItoKMCPhysics.H
@brief Main file for describing Ito-based plasma physics
@author Robert Marskar
*/
#ifndef CD_ItoKMCPhysics_H
#define CD_ItoKMCPhysics_H
// Std includes
#include <memory>
#include <vector>
// Chombo includes
#include <RealVect.H>
#include <RefCountedPtr.H>
#include <List.H>
// Our includes
#include <CD_ItoSpecies.H>
#include <CD_CdrSpecies.H>
#include <CD_RtSpecies.H>
#include <CD_Photon.H>
#include <CD_ItoParticle.H>
#include <CD_PointParticle.H>
#include <CD_ItoKMCPhotoReaction.H>
#include <CD_ItoKMCSurfaceReactionSet.H>
#include <CD_KMCDualState.H>
#include <CD_KMCDualStateReaction.H>
#include <CD_KMCSolver.H>
#include <CD_NamespaceHeader.H>
namespace Physics {
namespace ItoKMC {
/*!
@brief Numerical representation of the KMC state. Can be floating-point or integer type
*/
using FPR = Real;
/*!
@brief KMC state used in the Kinetic Monte Carlo advancement
*/
using KMCState = KMCDualState<FPR>;
/*!
@brief KMC reaction used in the Kinetic Monte Carlo advancement
*/
using KMCReaction = KMCDualStateReaction<KMCState, FPR>;
/*!
@brief KMC solverl type used in the Kinetic Monte carlo advancement
*/
using KMCSolverType = KMCSolver<KMCReaction, KMCState, FPR>;
/*!
@brief Map to species type
@details This is just for distinguishing between species that are treated with an Ito or CDR formalism
*/
enum SpeciesType
{
Ito,
CDR
};
/*!
@brief Base class for interaction between Kinetic Monte Carlo and Ito-based plasma solvers.
@details When using this class, the user should populate m_kmcReactions with the appropriate reactions AND implement routines
for computing the mobility/diffusion coefficeints.
*/
class ItoKMCPhysics
{
public:
/*!
@brief Constructor. Does nothing.
*/
inline ItoKMCPhysics() noexcept;
/*!
@brief Destructor. Does nothing.
*/
inline virtual ~ItoKMCPhysics() noexcept;
/*!
@brief Define the KMC solver and state
*/
inline void
defineKMC() const noexcept;
/*!
@brief Kill the KMC solver
*/
inline void
killKMC() const noexcept;
/*!
@brief Compute a time step to use.
@param[in] a_E Electric field.
@param[in] a_pos Position
@param[in] a_numParticles Number of reactive particles of each species
@return Default returns infinity but user can implement a physics based time step.
*/
inline virtual Real
computeDt(const RealVect a_E, const RealVect a_pos, const Vector<FPR> a_numParticles) const noexcept;
/*!
@brief Compute Townsend ionization coefficient
@param[in] a_E Electric field magnitude
@param[in] a_x Physical coordinate
@return Should return the Townsend ionization coefficient.
*/
virtual Real
computeAlpha(const Real a_E, const RealVect a_x) const = 0;
/*!
@brief Compute Townsend attachment coefficient
@param[in] a_E Electric field magnitude
@param[in] a_x Physical coordinate
@return Should return the Townsend attachment coefficient.
*/
virtual Real
computeEta(const Real a_E, const RealVect a_x) const = 0;
/*!
@brief Get all particle drift-diffusion species
@return m_itoSpecies
*/
const Vector<RefCountedPtr<ItoSpecies>>&
getItoSpecies() const;
/*!
@brief Get all fluid drift-diffusion species
@return m_cdrSpecies
*/
const Vector<RefCountedPtr<CdrSpecies>>&
getCdrSpecies() const;
/*!
@brief Get all photon species
@return m_rtSpecies
*/
const Vector<RefCountedPtr<RtSpecies>>&
getRtSpecies() const;
/*!
@brief Get number of plot variables
*/
virtual Vector<std::string>
getPlotVariableNames() const noexcept;
/*!
@brief Get plot variables
@param[in] a_E Electric field
@param[in] a_pos Physical position
@param[in] a_phi Plasma species densities
@param[in] a_gradPhi Density gradients for plasma species.
@param[in] a_dx Grid resolution
@param[in] a_kappa Cut-cell volume fraction
*/
virtual Vector<Real>
getPlotVariables(const RealVect a_E,
const RealVect a_pos,
const Vector<Real>& a_phi,
const Vector<RealVect>& a_gradPhi,
const Real a_dx,
const Real a_kappa) const noexcept;
/*!
@brief Get number of plot variables
*/
virtual int
getNumberOfPlotVariables() const noexcept;
/*!
@brief Return number of Ito solvers.
*/
inline int
getNumItoSpecies() const;
/*!
@brief Return number of CDR solvers.
*/
inline int
getNumCdrSpecies() const;
/*!
@brief Return total number of plasma species.
*/
inline int
getNumPlasmaSpecies() const;
/*!
@brief Return number of RTE solvers.
*/
inline int
getNumPhotonSpecies() const;
/*!
@brief Return true/false if physics model needs species gradients.
*/
virtual bool
needGradients() const noexcept;
/*!
@brief Get the internal mapping between plasma species and Ito solvers
*/
inline const std::map<int, std::pair<SpeciesType, int>>&
getSpeciesMap() const noexcept;
/*!
@brief Parse run-time options
*/
inline virtual void
parseRuntimeOptions() noexcept;
/*!
@brief Set initial surface charge. Default is 0, override if you want.
@param[in] a_time Simulation time
@param[in] a_pos Physical coordinate
*/
inline virtual Real
initialSigma(const Real a_time, const RealVect a_pos) const;
/*!
@brief Compute the Ito solver mobilities.
@param[in] a_time Time
@param[in] a_pos Position
@param[in] a_E Electric field
@return Must return a vector of non-negative mobility coefficients for the plasma species
*/
virtual Vector<Real>
computeMobilities(const Real a_time, const RealVect a_pos, const RealVect a_E) const noexcept = 0;
/*!
@brief Compute the Ito solver diffusion coefficients
@param[in] a_time Time
@param[in] a_pos Position
@param[in] a_E Electric field
@return Must return a vector of non-negative diffusion coefficients for the plasma species
*/
virtual Vector<Real>
computeDiffusionCoefficients(const Real a_time, const RealVect a_pos, const RealVect a_E) const noexcept = 0;
/*!
@brief Resolve secondary emission at the EB.
@details Routine is here to handle charge injection, secondary emission etc.
@param[out] a_secondaryParticles Outgoing plasma species particles.
@param[out] a_cdrFluxes CDR fluxes for CDR species.
@param[out] a_secondaryPhotons Photons injected through the EB.
@param[in] a_primaryParticles Particles that left the computational domain through the EB.
@param[in] a_cdrFluxesExtrap Extrapolated CDR fluxes.
@param[in] a_primaryPhotons Photons that left the computational domain through the EB.
@param[in] a_newNumParticles Total number of particles in the cut-cell AFTER the transport step.
@param[in] a_oldNumParticles Total number of particles in the cut-cell BEFORE the transport step.
@param[in] a_electricField Electric field.
@param[in] a_physicalCellCenter Physical position of the cell center.
@param[in] a_cellCentroid Cell centroid relative to the cell center (not multiplied by dx).
@param[in] a_bndryCentroid EB face centroid relative to the cell center (not multiplied by dx).
@param[in] a_bndryNormal Cut-cell normal vector.
@param[in] a_bndryArea Cut-cell boundary area - not multiplied by dx (2D) or dx^2 (3D).
@param[in] a_dx Grid resolution on this level.
@param[in] a_dt Time step.
@param[in] a_isDielectric Dielectric or electrode.
@param[in] a_matIndex Material index (taken from computationalGeometry).
*/
virtual void
secondaryEmissionEB(Vector<List<ItoParticle>>& a_secondaryParticles,
Vector<Real>& a_cdrFluxes,
Vector<List<Photon>>& a_secondaryPhotons,
const Vector<List<ItoParticle>>& a_primaryParticles,
const Vector<Real>& a_cdrFluxesExtrap,
const Vector<List<Photon>>& a_primaryPhotons,
const RealVect& a_E,
const RealVect& a_physicalCellCenter,
const RealVect& a_cellCentroid,
const RealVect& a_bndryCentroid,
const RealVect& a_bndryNormal,
const Real a_bndryArea,
const Real a_dx,
const Real a_dt,
const bool a_isDielectric,
const int a_matIndex) const noexcept = 0;
/*!
@brief Advance particles.
@param[inout] a_numParticles Number of physical particles
@param[out] a_numNewPhotons Number of new physical photons to generate (of each type)
@param[in] a_phi Plasma species densities.
@param[in] a_gradPhi Plasma species density gradients.
@param[in] a_dt Time step
@param[in] a_E Electric field
@param[in] a_pos Physical position
@param[in] a_dx Grid resolution
@param[in] a_kappa Cut-cell volume fraction.
*/
inline void
advanceKMC(Vector<FPR>& a_numParticles,
Vector<FPR>& a_numNewPhotons,
const Vector<Real>& a_phi,
const Vector<RealVect>& a_gradPhi,
const Real a_dt,
const RealVect a_E,
const RealVect a_pos,
const Real a_dx,
const Real a_kappa) const;
/*!
@brief Reconcile the number of particles.
@details This will add/remove particles and potentially also adjust the particle weights.
@param[inout] a_particles Computational particles
@param[in] a_newNumParticles New number of particles (i.e., after the KMC advance)
@param[in] a_oldNumParticles Previous number of particles (i.e., before the KMC advance)
@param[in] a_cellPos Cell center position
@param[in] a_centroidPos Cell centroid position
@param[in] a_loCorner Low corner of minimum box enclosing the cut-cell
@param[in] a_hiCorner High corner of minimum box enclosing the cut-cell
@param[in] a_bndryCentroid Cut-cell boundary centroid
@param[in] a_bndryNormal Cut-cell normal (pointing into the domain)
@param[in] a_dx Grid resolution
@param[in] a_kappa Cut-cell volume fraction.
@note Public because this is called by ItoKMCStepper
*/
inline void
reconcileParticles(Vector<List<ItoParticle>*>& a_particles,
const Vector<FPR>& a_newNumParticles,
const Vector<FPR>& a_oldNumParticles,
const RealVect a_cellPos,
const RealVect a_centroidPos,
const RealVect a_lo,
const RealVect a_hi,
const RealVect a_bndryCentroid,
const RealVect a_bndryNormal,
const Real a_dx,
const Real a_kappa) const noexcept;
/*!
@brief Generate new photons.
@details This will add photons
@param[in] a_newPhotons New photons
@param[in] a_numNewPhotons Number of physical photons to be generated.
@param[in] a_cellPos Cell center position
@param[in] a_centroidPos Cell centroid position
@param[in] a_loCorner Low corner of minimum box enclosing the cut-cell
@param[in] a_hiCorner High corner of minimum box enclosing the cut-cell
@param[in] a_bndryCentroid Cut-cell boundary centroid
@param[in] a_bndryNormal Cut-cell normal (pointing into the domain)
@param[in] a_dx Grid resolution
@param[in] a_kappa Cut-cell volume fraction.
@note Public because this is called by ItoKMCStepper
*/
inline void
reconcilePhotons(Vector<List<Photon>*>& a_newPhotons,
const Vector<FPR>& a_numNewPhotons,
const RealVect a_cellPos,
const RealVect a_centroidPos,
const RealVect a_lo,
const RealVect a_hi,
const RealVect a_bndryCentroid,
const RealVect a_bndryNormal,
const Real a_dx,
const Real a_kappa) const noexcept;
/*!
@brief Reconcile photoionization reactions.
@param[inout] a_itoPrticles Particle products placed in Ito solvers.
@param[inout] a_cdrParticles Particle products placed in CDR solvers.
@param[in] a_absorbedPhotons Photons absorbed on the mesh.
@details This runs through the photo-reactions and associates photo-ionization products.
*/
inline void
reconcilePhotoionization(Vector<List<ItoParticle>*>& a_itoParticles,
Vector<List<PointParticle>*>& a_cdrParticles,
const Vector<List<Photon>*>& a_absorbedPhotons) const noexcept;
protected:
/*!
@brief Enum for switching between KMC algorithms
@details 'SSA' is the Gillespie algorithm, 'Tau' is tau-leaping and 'Hybrid' is the Cao et. al. algorithm.
*/
enum class Algorithm
{
SSA,
TauPlain,
TauMidpoint,
HybridPlain,
HybridMidpoint,
};
/*!
@brief Enum for switching between various particle placement algorithms
*/
enum class ParticlePlacement
{
Random,
Centroid
};
/*!
@brief Algorithm to use for KMC advance
*/
Algorithm m_algorithm;
/*!
@brief Particle placement algorithm
*/
ParticlePlacement m_particlePlacement;
/*!
@brief Map for associating a plasma species with an Ito solver or CDR solver.
*/
std::map<int, std::pair<SpeciesType, int>> m_speciesMap;
/*!
@brief Class name. Used for options parsing
*/
std::string m_className;
/*!
@brief Turn on/off debugging
*/
bool m_debug;
/*!
@brief Is defined or not
*/
bool m_isDefined;
/*!
@brief Is the KMC solver defined or not.
*/
static thread_local bool m_hasKMCSolver;
/*!
@brief Kinetic Monte Carlo solver used in advanceReactionNetwork.
*/
static thread_local KMCSolverType m_kmcSolver;
/*!
@brief KMC state used in advanceReactionNetwork.
*/
static thread_local KMCState m_kmcState;
/*!
@brief KMC reactions used in advanceReactionNetowkr
@note This is set up via setupKMC in order toi ensure OpenMP thread safety when calling advanceReactionNetwork. The vector
is later depopulated in killKMC().
*/
static thread_local std::vector<std::shared_ptr<const KMCReaction>> m_kmcReactionsThreadLocal;
/*!
@brief List of reactions for the KMC solver
*/
std::vector<KMCReaction> m_kmcReactions;
/*!
@brief List of photoionization reactions
*/
std::vector<ItoKMCPhotoReaction> m_photoReactions;
/*!
@brief Random number generators for photoionization pathways.
@details The first index is the photon index, i.e. entry in m_photonSpecies. The second index in the
map is an RNG generator for selecting photo-reactions, and a map which associates the returned
reaction from the RNG generator with an index in m_photoReactions.
*/
std::map<int, std::pair<std::discrete_distribution<int>, std::map<int, int>>> m_photoPathways;
/*!
@brief Surface reactions.
*/
ItoKMCSurfaceReactionSet m_surfaceReactions;
/*!
@brief List of solver-tracked particle drift-diffusion species.
*/
Vector<RefCountedPtr<ItoSpecies>> m_itoSpecies;
/*!
@brief List of solver-tracked fluid drift-diffusion species.
*/
Vector<RefCountedPtr<CdrSpecies>> m_cdrSpecies;
/*!
@brief List of solver-tracked photon species.
*/
Vector<RefCountedPtr<RtSpecies>> m_rtSpecies;
/*!
@brief Maximum new number of particles generated by the chemistry advance
*/
int m_maxNewParticles;
/*!
@brief Maximum new number of photons generated by the chemistry advance
*/
int m_maxNewPhotons;
/*!
@brief Solver setting for the Cao et. al algorithm.
@details Determines critical reactions. A reaction is critical if it is m_Ncrit firings away from depleting a reactant.
*/
int m_Ncrit;
/*!
@brief Solver setting for the Cao et. al algorithm.
@details Maximum number of SSA steps to run when switching into SSA-based advancement for non-critical reactions.
*/
int m_NSSA;
/*!
@brief Solver setting for the Cao et. al. algorithm.
@details Equal to the maximum permitted change in the relative propensity for non-critical reactions
*/
Real m_SSAlim;
/*!
@brief Solver setting for the Cao et. al. algorithm.
@details Equal to the maximum permitted change in the relative propensity for non-critical reactions
*/
Real m_eps;
/*!
@brief Define method -- defines all the internal machinery
*/
inline void
define() noexcept;
/*!
@brief Build internal representation of how we distinguish the Ito and CDR solvers
@note This should ALWAYS be called after initializing the species since ItoKMCStepper will rely on it.
*/
inline void
defineSpeciesMap() noexcept;
/*!
@brief Define pathways for photo-reactions
*/
inline void
definePhotoPathways() noexcept;
/*!
@brief Parse the maximum number of particles generated per cell
*/
inline void
parsePPC() noexcept;
/*!
@brief Parse the maximum number of particles generated per cell
*/
inline void
parseDebug() noexcept;
/*!
@brief Parse reaction algorithm
*/
inline void
parseAlgorithm() noexcept;
/*!
@brief Update reaction rates
@param[out] a_kmcReactions Reaction rates to be set.
@param[in] a_E Electric field
@param[in] a_pos Physical position
@param[in] a_phi Plasma species densities
@param[in] a_gradPhi Density gradients for plasma species.
@param[in] a_dx Grid resolution
@param[in] a_kappa Cut-cell volume fraction
@note Must be implemented by the user.
*/
virtual void
updateReactionRates(std::vector<std::shared_ptr<const KMCReaction>>& a_kmcReactions,
const RealVect a_E,
const RealVect a_pos,
const Vector<Real>& a_phi,
const Vector<RealVect>& a_gradPhi,
const Real a_dx,
const Real a_kappa) const noexcept = 0;
/*!
@brief Remove particles from the input list.
@details This will remove weight from the input particles if we can. Otherwise we remove full particles.
@param[inout] a_particles List of (super-)particles to remove from.
@param[in] a_numToRemove Number of physical particles to remove from the input list
*/
inline void
removeParticles(List<ItoParticle>& a_particles, const long long a_numToRemove) const;
};
} // namespace ItoKMC
} // namespace Physics
#include <CD_NamespaceFooter.H>
#include <CD_ItoKMCPhysicsImplem.H>
#endif
Species definitions
Species are defined either as input species for CDR solvers (see CdrSolver) or Îto solvers (see ItoSolver).
It is sufficient to populate the ItoKMCPhysics species vectors m_cdrSpecies and m_itoSpecies if one only wants to initialize the solvers.
See ItoKMCPhysics for their C++ definition.
In addition to actually populating the vectors, users will typically also include initial conditions for the species.
The interfaces permit initial particles for the Îto solvers, while the CDR solvers permit particles and initial density functions.
Additionally, one must define all species associated with radiative transfer solvers.
Plasma reactions
Plasma reactions in the Îto-KMC model are represented stoichiometrically as
where \(c\) is the KMC rate coefficient, i.e. not the rate coefficient from the reaction rate equation. An arbitrary number of reactions is supported, but currently the reaction mechanisms are limited to:
The left-hand side must consist only of species that are tracked by a CDR or Îto solver.
The right-hand side can consist of species tracked by a CDR solver, an Îto solver, or a radiative transfer solver.
These rules apply only to the internal C++ interface; implementations of this interface will typically also allow species defined as ‘’background species’’, i.e. species that are not tracked by a solver. In that case, the interface implementation must absorb the background species contribution into the rate constant.
In the KMC algorithm, one operates with chemical propensities rather than rate coefficients. For example, the chemical propensity \(a_1\) for a reaction \(A + B \xrightarrow{c_1}\varnothing\) is
and from the macroscopic limit \(X_A \gg 1\), \(X_B \gg 1\) one may in fact also derive that \(c_1 = k_1/\Delta V\) where \(k_1\) is the rate coefficient from the reaction rate equation, i.e. where \(\partial_t n_A = -k_1n_An_B\). Note that if the reaction was \(A + A\xrightarrow{c_2}\varnothing\) then the propensity is
since there are \(\frac{1}{2}X_A\left(X_A-1\right)\) distinct pairs of particles of type \(A\). Likewise, the fluid rate coefficient would be \(k_2 = c_2\Delta V/2\).
The distinction between KMC and fluid rates is an important one; the reaction representation used in the Îto-KMC model only operates with the KMC rates \(c_j\), and it is us to the user to ensure that these are consistent with the fluid limit.
Internally, these reactions are implemented through the dual state KMC implementation, see KMCDualState.
During the reaction advance the user only needs to update the \(c_j\) coefficients (typically done via an interface implementation); the calculation of the propensity is automatic and follows the standard KMC rules (e.g., the KMC solver accounts for the number of distinct pairs of particles).
This must be done in the routine updateReactionRates(...), see ItoKMCPhysics for the complete specification.
Photoionization
Photo-reactions are also represented stoichiometrically as
where \(\xi\) is the photo-reaction efficiency, i.e. the probability that the photon \(\gamma\) causes a reaction when it is absorbed on the mesh. Currently, photons can only be generated through plasma reactions, see Plasma reactions, and we do not permit reactions between the photon and plasma species.
Important
We currently only support constant photoionization probabilities \(\xi\).
ItoKMCPhysics adds some flexibility when dealing with photo-reactions as it permits pre-evaluation of photoionization probabilities.
That is, when generating photons through reactions of the type \(A + B \rightarrow \gamma\), one may scale the reaction by the photoionization probability \(\xi\) so that only ionizing photons are generated and then write the photo-reaction as \(\gamma \xrightarrow{\xi=1} A + B + \ldots\).
However, this process becomes more complicated when dealing with multiple photo-reaction pathways, e.g.,
where \(\xi_1\) and \(\xi_2\) are the probabilities that the photon \(\gamma\) triggers the reaction. If only a single physical photon is absorbed, then one must stochastically determine which pathway is triggered. There are two ways of doing this:
Pre-scale the photon generation by \(\xi_1 + \xi_2\), and then determine the pathway through the relative probabilities \(\xi_1/\left(\xi_1+\xi_2\right)\) and \(\xi_2/\left(\xi_1+\xi_2\right)\).
Do not scale the photon generation reaction and determine the pathway through the absolute probabilities \(\xi_1\) and \(\xi_2\). This can imply a large computational cost since one will have to track all photons that are generated.
The former method is normally the preferred way as it can lead to reductions in computational complexity and particle noise, but for flexibility ItoKMCPhysics supports both of these.
The latter method can be of relevance if users wants precise descriptions of photons that trigger both photoionizing reactions and surface reactions (e.g., secondary electron emission).
An alternative is to split the photon type \(\gamma\) into volumetrically absorbed and surface absorbed photon species \(\gamma_v\) and \(\gamma_s\). In this case one may pre-scale the photon-generating reactions by the photo-reaction probabilities \(\xi\) (for volume absorption) and \(\zeta\) (for surface absorption) as follows:
Volumetric and surface absorption is then treated independently
This type of pre-evaluation of the photo-reaction pathways is sensible in a statistical sense, but loses meaning if only a single photon is involved.
Surface reactions
Warning
Surface reactions are supported by ItoKMCPhysics but not implemented in the JSON interface (yet).
Transport coefficients
Species mobilities and diffusion coefficients should be computed just as they are done in the fluid approximation.
The ItoKMCPhysics interface requires implementations of two functions that define the coefficients as functions of \(\mu = \mu\left(t,\mathbf{x}, \mathbf{E}\right)\),
and likewise for the diffusion coefficients.
Note that these functions should return the fluid coefficients.
Important
There is currently no support for computing \(\mu\) as a function of the species densities (e.g., the electron density), but this only requires modest extensions of the Îto-KMC module.
JSON 0D chemistry interface
The JSON-based chemistry interface (called ItoKMCJSON) simplifies the definition of the plasma kinetics and initial conditions by permitting specifications through a JSON file.
In addition, the JSON specification will permit the definition of background species that are not tracked as solvers (but that simply exist as a density function).
The JSON interface thus provides a simple entry point for users that have no interest in defining the chemistry from scratch.
Several mandatory fields are required when specifying the plasma kinetics, such as:
Gas pressure, temperature, and number density.
Background species (if any).
Townsend coefficients.
Plasma species and their initial conditions, which solver to use when tracking them.
Photon species and their absorption lengths.
Plasma reactions.
Photoionization reactions.
These specifications follow a pre-defined JSON format which is discussed in detail below.
Tip
While JSON files do not formally support comments, the JSON submodule used by chombo-discharge allows them.
Gas law
The JSON gas law represents a loose definition of the gas pressure, temperature, and number density.
Multiple gas laws can be defined in the JSON file, and the user can then select which one to use.
This is defined within a JSON entry gas and a subentry law.
In the gas/law field one must specify an id field which indicates which gas law to use.
One may also set an optional field plot to true/false if one wants to include the pressure, temperature, and density in the HDF5 output files.
Each user-defined gas law exists as a separate entry in the gas/law field, and the minimum requirements for these are that they contain a type specifier.
Currently, the supported types are
ideal(for constant ideal gas)table vs heightfor tabulated density, temperature and pressure vs some height/axis.
The specification rules for these are different, and are discussed below.
An example JSON specification is given below.
Here, we have defined two gas law my_ideal_gas and tabulated_atmosphere and specified through the id parameter that we will use the my_ideal_gas specification.
Example gas law specification
{
"gas" : {
// Specify which gas law we use. The user can define multiple gas laws and then later specify which one to use.
"law" : {
"id" : "my_ideal_gas", // Specify which gas law we use.
"plot" : true, // Turn on/off plotting.
"my_ideal_gas" : {
// Definition for an ideal gas law. Temperature is always in Kevlin the pressure is in bar. The neutral density
// is derived using an ideal gas law.
"type" : "ideal",
"temperature" : 300,
"pressure" : 1.0
},
"tabulated_atmosphere" : {
// Tabulated gas law. The user must supply a file containing the pressure, temperature and number density and
// specify the table resolution that will be used internally.
"type" : "table vs height", // Specify that our gas law contains table-vs-height data
"file" : "ENMSIS_Atmosphere.dat", // File name
"axis" : "y", // Associated Cartesian axis for the "height"
"temperature column" : 0, // Column containing the temperature
"pressure column" : 1, // Column containing the pressure
"density column" : 2, // Column containing the number density
"min height" : 0.0, // Minium height kept when resampling the table
"max height" : 250000, // Maximum height kept when resampling the table
"res height" : 500, // Table resolution
"dump tables" : true, // Optional argument for dumping tables to disk (useful for debugging)
"T scale" : 1.0, // Optional argument for user-specified temperature data scaling.
"P scale" : 1.0, // Optional argument for user-specified pressure data scaling.
"Rho scale" : 1.0 // Optional argument for user-specified density data scaling.
}
}
}
}
Ideal gas
When specifying that the gas law type is ideal, the user must further specify the temperature and pressure of the gas.
Example ideal gas specification
{
"gas" : {
"law" : {
"id" : "my_ideal_gas", // Specify which gas law we use.
"my_ideal_gas" : {
"type" : "ideal", // Ideal gas law type
"temperature" : 300, // Temperature must be in Kelvin
"pressure" : 1.0 // Pressure must be in bar
},
}
}
}
Table vs height
Background species
Background species \(i\) are defined as continuous background densities \(N_i\) given by
where \(\chi_i\left(\mathbf{x}\right)\) is the molar fraction of species \(i\) (typically, but not necessarily, a constant).
When specifying background species, the user must include an array of background species in the gas JSON field.
For example,
{
"gas" : {
"background species" : [
{
// First species goes here
},
{
// Second species goes here
}
]
}
}
The order of appearance of the background species does not matter.
Name
The species name is specified by including an id field, e.g.
{
"gas" : {
"background species" : [
{
"id": "O2"
}
]
}
}
Molar fraction
The molar fraction can be specified in various forms by including a field molar fraction.
This field also requires the user to specify the type of molar fraction.
In principle, the molar fraction can have a positional dependency.
Warning
There’s no internal renormalization of the molar fractions input by the user. Internal inconsitencies will occur if the user supplies inputs molar fractions that sum to a number different than one.
Various checks will be enabled if the user compiles in debug-mode (see Installation), but these checks are not guaranteed to catch all cases.
Constant
To specify a constant molar fraction, set the type specifier to constant and specify the value.
An example JSON specification is
{
"gas" : {
"background species" : [
{
"id": "O2" // Species name
"molar fraction": { // Molar fraction specification
"type": "constant", // Constant molar fraction
"value": 0.2 // Molar fraction value
}
}
]
}
}
Tabulated versus height
The molar fraction can be set as a tabulated value versus one of the Cartesian coordinate axis by setting the type specifier to table vs height.
The input data should be tabulated in column form, e.g.
# height molar fraction
0 0.1
1 0.1
2 0.1
The file parser (see LookupTable) will ignore the header file if it starts with a hashtag (#). Various other inputs are then also required:
fileFile name containing the height vs. molar fraction data (required).axisCartesian axis to associate with the height coordinate (required).height columnColumn in data file containing the height (optional, defaults to 0).molar fraction columnColumn in data file containing the molar fraction (optional, defaults to 1).height scaleScaling of height column (optional).fraction scaleScaling (optional).min heightTruncate data to minimum height (optional, applied after scaling).max heightTruncate data to maximum height (optional, applied after scaling).num pointsNumber of points to keep in table representation (optional, defaults to 1000).spacingSpacing table representation (optional, defaults to “linear”.dumpOption for dumping tabulated data to file.
An example JSON specification is
{
"gas" : {
"background species" : [
{
"id": "O2" // Species name
"molar fraction": { // Molar fraction specification
"type": "table vs height", // Constant molar fraction
"file" : "O2_molar_fraction.dat", // File name containing molar fraction
"dump" : "debug_O2_fraction.dat", // Optional debugging hook for dumping table to file
"axis" : "y", // Axis which represents the "height"
"height column" : 0, // Optional specification of column containing the height data (defaults to 0)
"molar fraction column" : 1, // Optional specification of column containing the height data (defaults to 1)
"height scale" : 1.0, // Optional scaling of height column
"fraction scale" : 1.0, // Optional scaling of molar fraction column
"min height" : 0.0, // Optional truncation of minimum height kept in internal table (occurs after scaling)
"max height" : 2.0, // Optional truncation of maximum height kept in internal table (occurs after scaling)
"num points" : 100, // Optional specification of number of data points kept in internal table (defaults to 1000)
"spacing" : "linear" // Optional specification of table representation. Can be 'linear' or 'exponential' (defaults to linear)
}
}
]
}
}
Plotting
Background species densities can be included in HDF5 plots by including an optional field plot.
For example
{
"gas" : {
"background species" : [
{
"id": "O2" // Species name
"molar fraction": { // Molar fraction specification
"type": "constant", // Constant molar fraction
"value": 0.2 // Molar fraction value
},
"plot": true // Plot number density in HDF5 or not
}
]
}
}
Townsend coefficients
Townsend ionization and attachment coefficients \(\alpha\) and \(\eta\) must be specified, and have two usages:
Flagging cells for refinement (users can override this).
Potential usage in plasma reactions (which requires particular care, see Warnings and caveats.
These are specified by including JSON entries alpha and eta, e.g.
{
"alpha": {
// alpha-coefficient specification goes here
},
"eta" : {
// eta-coefficient specification goes here
}
}
There are various way of specifying these, as discussed below:
Automatic
In auto-mode the Townsend coefficients for the electrons is automatically derived from the user-specified list of reactions. E.g., the Townsend ionization coefficient is computed as
where \(\sum k\) is the sum of all ionizing reactions, also incorporating the neutral density. The advantage of this approach is that one may choose to discard some of the reactions from e.g. BOLSIG+ output but without recomputing the Townsend coefficients.
To automatically set Townsend coefficients, set type to auto and specify the species that is involved, e.g.
{
"alpha": {
"type": "auto",
"species": "e"
}
}
The advantage of using auto-mode for the coefficients is that one automatically ensures consistency between the user-specified list of reactions and the
Constant
To set a constant Townsend coefficient, set type to constant and then specify the value, e.g.
{
"alpha": {
"type": "constant",
"value": 1E5
}
}
Tabulated vs \(E/N\)
To set the coefficient as functions \(f = f\left(E/N\right)\), set the type specifier to table vs E/N and specify the following fields:
fileFor specifying the file containing the input data, which must be organized as column data, e.g.# E/N alpha/N 0 1E5 100 1E5 200 1E5
headerFor specifying where in the file one starts reading data. This is an optional argument intended for use with BOLSIG+ output data where the user specifies that the data is contained in the lines below the specified header.E/N columnFor setting which column in the data file contains the values of \(E/N\) (optional, defaults to 0).alpha/N columnFor setting which column in the data file contains the values of \(\alpha/N\). Note that this should be replaced byeta/N columnwhen parsing the attachment coefficient. This is an optional argument that defaults to 1.scale E/NOptional scaling of the column containing the \(E/N\) data.scale alpha/NOptional scaling of the column containing the \(\alpha/N\) data.min E/NOptional truncation of the table representation of the data (applied after scaling).max E/NOptional truncation of the table representation of the data (applied after scaling).num pointsNumber of data points in internal table representation (optional, defaults to 1000).spacingSpacing in internal table representation. Can belinearorexponential. This is an optional argument that defaults toexponential.dumpA string specifier for writing the table representation of \(E/N, \alpha/N\) to file.
An example JSON specification that uses a BOLSIG+ output file for parsing the data is
{
"alpha": {
"type" : "table vs E/N", // Specify that we read in alpha/N vs E/N
"file" : "bolsig_air.dat", // File containing the data
"dump" : "debug_alpha.dat", // Optional dump of internalized table to file. Useful for debugging.
"header" : "E/N (Td)\tTownsend ioniz. coef. alpha/N (m2)", // Optional argument. Contains line immediately preceding the data to be read.
"E/N column" : 0, // Optional specification of column containing E/N (defaults to 0)
"alpha/N column" : 1, // Optional specification of column containing alpha/N (defaults to 1)
"min E/N" : 1.0, // Optional truncation of minium E/N kept when resampling the table (occurs after scaling)
"max E/N" : 1000.0, // Optional truncation of maximum E/N kept when resampling the table (happens after scaling)
"num points" : 1000, // Optional number of points kept when resamplnig the table (defaults to 1000)
"spacing" : "exponential", // Optional spcification of table representation. Defaults to 'exponential' but can also be 'linear'
"scale E/N" : 1.0, // Optional scaling of the column containing E/N
"scale alpha/N" : 1.0 // Optional scaling
}
}
Plotting
To include the Townsend coefficients as mesh variables in HDF5 files, include the plot specifier, e.g.
{
"alpha": {
"type": "auto",
"species": "e",
"plot": true
}
}
Plasma species
Plasma species are species that are tracked using either a CDR or Îto solver. The user must specify the following information:
An ID/name for the species.
The species’ charge number.
The solver type, i.e. whether or not it is tracked by a particle or fluid solver.
Whether or not the species is mobile/diffusive. If the species is mobile/diffusive, the mobility/diffusion coefficient must also be specified.
Optionally, the species temperature. If not specified, the temperature will be set to the background gas temperature.
Initial particles for the solvers.
Basic definition
Species are defined by an array of entries in a plasma species JSON field.
Each species must specify the following fields
id(string) For setting the species name.Z(integer) For setting the species’ charge number.solver(string) For specifying whether or not the species is tracked by a CDR or Îto solver. Acceptable values are cdr and ito.mobile(boolean) For specifying whether or not the species is mobile.mobile(boolean) For specifying whether or not the species is diffusive.
An example specification is
"plasma species" :
[
// List of plasma species that are tracked. This is an array of species
// with various identifiers, some of which are always required (id, Z, type, mobile, diffusive) and
// others which are secondary requirements.
{
"id": "e", // Species ID
"Z" : -1, // Charge number
"solver" : "ito", // Solver type. Either 'ito' or 'cdr'
"mobile" : true, // Mobile or not
"diffusive" : true // Diffusive or not
},
{
// Definition of O2+ plasma species.
"id": "O2+", // Species ID
"Z" : 1, // Charge number
"solver" : "cdr", // CDR solver.
"mobile" : false, // Not mobile.
"diffusive" : false // Not diffusive
}
]
Note that the order of appearance of the various species is irrelevant. However, if a species is specified as mobile/diffusive, the user must also specify the corresponding transport coefficients.
Mobility and diffusion coefficients
Mobility and diffusion coefficients are specified by including fields mobility and diffusion that further specifies how the transport coefficients are calculated.
Note that the input to these fields should be fluid transport coefficients.
These fields can be specified in various forms, as shown below.
Constant
To set a constant coefficient, set the type specifier to constant and then assign the value.
For example,
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : -1, // Charge number
"solver" : "ito", // Solver type. Either 'ito' or 'cdr'
"mobile" : true, // Mobile or not
"diffusive" : true // Diffusive or not,
"mobility": {
"type": "constant", // Set constant mobility
"value": 0.02
},
"diffusion": {
"type": "constant", // Set constant diffusion coefficient
"value": 2E-4
}
}
]
Constant \(*N\)
To set a coefficient that is constant vs \(N\), set the type specifier to constant mu*N or constant D*N and then assign the value.
For example,
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : -1, // Charge number
"solver" : "ito", // Solver type. Either 'ito' or 'cdr'
"mobile" : true, // Mobile or not
"diffusive" : true // Diffusive or not,
"mobility": {
"type": "constant mu*N", // Set mu*N to a constant
"value": 1E24
},
"diffusion": {
"type": "constant D*N", // Set D*N to a constant
"value": 5E24
}
}
]
Table vs \(E/N\)
To set the transport coefficients as functions \(f = f\left(E/N\right)\), set the type specifier to table vs E/N and specify the following fields:
fileFor specifying the file containing the input data, which must be organized as column data, e.g.# E/N mu/N 0 1E5 100 1E5 200 1E5
headerFor specifying where in the file one starts reading data. This is an optional argument intended for use with BOLSIG+ output data where the user specifies that the data is contained in the lines below the specified header.E/N columnFor setting which column in the data file contains the values of \(E/N\) (optional, defaults to 0).mu*N columnFor setting which column in the data file contains the values of \(\mu N\) (or alternatively \(DN\) for the diffusion coefficient). This is an optional argument that defaults to 1.E/N scaleOptional scaling of the column containing the \(E/N\) data.mu*N scaleOptional scaling of the column containing the \(\mu N\) data (or alternatively \(DN\) for the diffusion coefficient).min E/NOptional truncation of the table representation of the data (applied after scaling).max E/NOptional truncation of the table representation of the data (applied after scaling).num pointsNumber of data points in internal table representation (optional, defaults to 1000).spacingSpacing in internal table representation. Can belinearorexponential. This is an optional argument that defaults toexponential.dumpA string specifier for writing the table representation of \(E/N, \mu N\) to file.
An example JSON specification that uses a BOLSIG+ output file for parsing the data for the mobility is
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : -1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : true, // Mobile
"diffusive" : false, // Not diffusive
"mobility" : {
"type" : "table vs E/N", // Specification of tabulated mobility lookup method
"file" : "bolsig_air.dat", // File containg the mobility data
"dump" : "debug_mobility.dat", // Optional argument for dumping table to disk (useful for debugging)
"header" : "E/N (Td)\tMobility *N (1/m/V/s)", // Line immediately preceding the colum data
"E/N column" : 0, // Column containing E/N
"mu*N column" : 1, // Column containing mu*N
"min E/N" : 1, // Minimum E/N kept when resampling table
"max E/N" : 2E6, // Maximum E/N kept when resampling table
"points" : 1000, // Number of grid points kept when resampling the table
"spacing" : "exponential", // Grid point spacing when resampling the table
"E/N scale" : 1.0, // Optional argument for scaling mobility coefficient
"mu*N scale" : 1.0 // Optional argument for scaling the E/N column
}
}
]
Tip
The parser for the diffusion coefficient is analogous; simply replace mu*N by D*N.
Temperature
The species temperature is only parametrically attached to the species (i.e., not solved for), and can be specified by the user.
Note that the temperature is mostly relevant for transport coefficients (e.g., reaction rates).
The temperature can be specified through the temperature field for each species, and various forms of specifying this is available.
Important
If the specifies temperature is not specified by the user, it will be set equal to the background gas temperature.
Background gas
To set the temperature equal to the background gas temperature, set the type specifier field to gas.
For example:
"plasma species" :
[
{
"id": "O2+", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"temperature": { // Specify temperature
"type": "gas" // Temperature same as gas
}
}
]
Constant
To set a constant temperature, set the type specifier to constant and then specify the temperature.
For example:
"plasma species" :
[
{
"id": "O2+", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"temperature": { // Specify temperature
"type": "constant", // Constant temperature
"value": 300 // Temperature in Kelvin
}
}
]
Table vs \(E/N\)
To set the species temperature as a function \(T = T\left(E/N\right)\), set the type specifier to table vs E/N and specify the following fields:
fileFor specifying the file containing the input data, which must be organized as column data versus the mean energy, e.g.# E/N energy (eV) 0 1 100 2 200 3
headerFor specifying where in the file one starts reading data. This is an optional argument intended for use with BOLSIG+ output data where the user specifies that the data is contained in the lines below the specified header.E/N columnFor setting which column in the data file contains the values of \(E/N\) (optional, defaults to 0).eV columnFor setting which column in the data file contains the energy vs \(E/N\). This is an optional argument that defaults to 1.E/N scaleOptional scaling of the column containing the \(E/N\) data.eV scaleOptional scaling of the column containing the energy (in eV).min E/NOptional truncation of the table representation of the data (applied after scaling).max E/NOptional truncation of the table representation of the data (applied after scaling).num pointsNumber of data points in internal table representation (optional, defaults to 1000).spacingSpacing in internal table representation. Can belinearorexponential. This is an optional argument that defaults toexponential.dumpA string specifier for writing the table representation of \(E/N, eV\) to file.
An example JSON specification that uses a BOLSIG+ output file for parsing the data for the mobility is
"plasma species" :
[
{
"id": "O2+", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : true, // Not mobile
"diffusive" : false, // Not diffusive
"temperature": { // Specification of temperature
"type": "table vs E/N", // Tabulated
"file": "bolsig_air.dat", // File name
"dump": "debug_temperature.dat", // Dump to file
"header" : "E/N (Td)\tMean energy (eV)", // Header preceding data
"E/N column" : 0, // Column containing E/N
"eV column" : 1, // Column containing the energy
"min E/N" : 10, // Truncation of table
"max E/N" : 2E6, // Truncation of table
"E/N scale": 1.0, // Scaling of input data
"eV scale": 1.0, // Scaling of input data
"spacing" : "exponential", // Table spacing
"points" : 1000 // Number of points in table
}
}
]
Initial particles
Initial particles for species are added by including an initial particles field.
The field consists of an array of JSON entries, where each array entry represents various ways of adding particles.
Important
The initial particles field is incrementing, each new entry will add additional particles.
For example, to add two particles with particle weights of 1, one may specify
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : false, // Not mobile
"diffusive" : false, // Not diffusive
"initial particles": [ // Initial particles
{
"single particle": { // Single physical particle at (0,0,0)
"position": [0, 0, 0],
"weight": 1.0
}
},
{
"single particle": { // Single physical particle at (0,1,0)
"position": [0, 1, 0],
"weight": 1.0
}
}
]
}
]
The various supported ways of additional initial particles are discussed below.
Important
If a solver is specified as a CDR solver, the initial particles are deposited as densities on the mesh using an NGP scheme.
Single particle
Single particles are placed by including a single particle JSON entry.
The user must specify
positionThe particle position.weightThe particle weight.
For example:
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : false, // Not mobile
"diffusive" : false, // Not diffusive
"initial particles": [ // Initial particles
{
"single particle": { // Single physical particle at (0,0,0)
"position": [0, 0, 0],
"weight": 1.0
}
}
]
}
]
Uniform distribution
Particles can be randomly drawn from a uniform distribution (a rectangular box in 2D/3D) by including a uniform distribution.
The user must specify
low cornerThe lower left corner of the box.high cornerThe lower left corner of the box.num particlesNumber of computational particles that are drawn.weightParticle weights.
For example:
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : false, // Not mobile
"diffusive" : false, // Not diffusive
"initial particles": [ // Initial particles
{
"uniform distribution" : { // Particles inside a box
"low corner" : [ -0.04, 0, 0 ], // Lower-left physical corner of sampling region
"high corner" : [ 0.04, 0.04, 0 ], // Upper-right physical corner of sampling region
"num particles": 1000, // Number of computational particles
"weight" : 1 // Particle weights
}
}
]
}
]
Sphere distribution
Particles can be randomly drawn inside a circle (sphere in 3D) by including a sphere distribution.
The user must specify
centerThe sphere center.radiusThe sphere radius.num particlesNumber of computational particles that are drawn.weightParticle weights.
For example:
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : false, // Not mobile
"diffusive" : false, // Not diffusive
"initial particles": [ // Initial particles
{
"sphere distribution" : { // Particles inside sphere
"center" : [ 0, 7.5E-3, 0 ], // Sphere center
"radius" : 1E-3, // Sphere radius
"num particles": 1000, // Number of computational particles
"weight" : 1 // Particle weights
}
}
]
}
]
Gaussian distribution
Particles can be randomly drawn from a Gaussian distribution by including a gaussian distribution.
The user must specify
centerThe sphere center.radiusThe sphere radius.num particlesNumber of computational particles that are drawn.weightParticle weights.
For example:
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : false, // Not mobile
"diffusive" : false, // Not diffusive
"initial particles": [ // Initial particles
{
"gaussian distribution" : { // Particles drawn from a gaussian
"center" : [ 0, 7.5E-3, 0 ], // Center
"radius" : 1E-3, // Radius
"num particles": 1000, // Number of computational particles
"weight" : 1 // Particle weights
}
}
]
}
]
List/file
Particles can be read from a file by including a list specifier.
The particles must be organized as rows containing the coordinate positions and weights, e.g.
# x y z w
0 0 0 1
The user must specify
fileFile name.xColumn containing the \(x\) coordinates (optional, defaults to 0).yColumn containing the \(y\) coordinates (optional, defaults to 1).zColumn containing the \(z\) coordinates (optional, defaults to 2).wColumn containing the particle weight (optional, defaults to 3).
For example:
"plasma species" :
[
{
"id": "e", // Species ID
"Z" : 1, // Charge number
"solver" : "ito", // Solver type.
"mobile" : false, // Not mobile
"diffusive" : false, // Not diffusive
"initial particles": [ // Initial particles
{
"list": {
"file": "initial_particles.dat", // File containing the particles
"x column": 0,
"y column": 1,
"z column": 2,
"w column": 3
}
}
]
}
]
Photon species
Photon species are species that are tracked using a radiative transfer solver.
These are defined by an array of entries in a photon species JSON entry, and must define the following information:
An ID/name for the species.
The absorption coefficient for the photon type.
Basic definition
A basic definition of a single photon species is
"photon species":
[
{
"id": "Y", // Photon species id. Must be unique
"kappa": { // Specification of absorption coefficient.
"type": "constant", // Specify constant absorption coefficient
"value": 1E4 // Value of the absorption coefficeint
}
}
]
Multiple photon species are added by appending with more entries, e.g.
"photon species":
[
{
"id": "Y1", // Photon species id. Must be unique
"kappa": { // Specification of absorption coefficient.
"type": "constant", // Specify constant absorption coefficient
"value": 1E4 // Value of the absorption coefficeint
}
},
{
"id": "Y2", // Photon species id. Must be unique
"kappa": { // Specification of absorption coefficient.
"type": "constant", // Specify constant absorption coefficient
"value": 1E4 // Value of the absorption coefficeint
}
}
]
Absorption coefficient
Constant
In order to set a constant absorption coefficient, set type to constant and then specify the value field.
For example:
"photon species":
[
{
"id": "Y",
"kappa": {
"type": "constant",
"value": 1E4
}
}
]
stochastic A
The stochastic A specifier computes a random absorption length from the expression
where \(p_i\) is the partial pressure of some species \(i\), and \(p_i\chi_{\textrm{min}}\) and \(p_i\chi_{\textrm{max}}\) are minimum and maximum absorption lengths on the frequency interval \(f\in[f_1,f_2]\).
The user must specify \(\chi_{\textrm{min}}\), \(\chi_{\textrm{max}}\), \(f_1\), \(f_2\), and the background species used when computing \(p_i\).
This is done through fields f1, f_2, chi min, chi max, and neutral.
For example:
"photon species":
[
{
"id": "Y",
"kappa": {
"type": "stochastic A",
"f1": 2.925E15,
"f2": 3.059E15,
"chi min": 2.625E-2,
"chi max": 1.5,
"neutral": "O2"
}
}
]
stochastic B
The stochastic B specifier computes a random absorption length from the expression
where \(p_i\) is the partial pressure of some species \(i\), and \(p_i\chi_{\textrm{min}}\) and \(p_i\chi_{\textrm{max}}\) are minimum and maximum absorption lengths, and \(u\) is a random number between 0 and 1.
Note that that this is just a simpler way of using the stochastic A specifier above.
The user must specify \(\chi_{\textrm{min}}\), \(\chi_{\textrm{max}}\), and the background species used when computing \(p_i\).
This is done through fields chi min, chi max, and neutral.
For example:
"photon species":
[
{
"id": "Y",
"kappa": {
"type": "stochastic B",
"chi min": 2.625E-2,
"chi max": 1.5,
"neutral": "O2"
}
}
]
Plasma reactions
Plasma reactions are always specified stoichiometrically in the form \(S_A + S_B + \ldots \xrightarrow{k} S_C + S_D + \ldots\). The user must supply the following information:
Species involved on the left-hand and right-hand side of the reaction.
How to compute the reaction rate.
Optionally specify gradient corrections to the reaction.
Optionally specify whether or not the reaction rate will be written to file.
Important
The JSON interface expects that the reaction rate \(k\) is the same as in the reaction rate equation. Internally, the rate is converted to a format that is consistent with the KMC algorithm.
Plasma reactions should be entered in the JSON file using an array (the order of reactions does not matter).
Basic reaction example specifiers
"plasma reactions": [
{
// First reaction goes here,
},
{
// Next reaction goes here,
}
]
Reaction specifiers
Each reaction must include an entry reaction in the list of reactions.
This entry must consist of a single string with a left-hand and right-hand sides separated by ->.
For example:
"plasma reactions": [
{
"reaction": "e + N2 -> e + e + N2+"
}
]
which is equivalent to the reaction \(\text{e} + \text{N}_2 \rightarrow \text{e} + \text{e} + \text{N}_2^+\).
Each species must be separated by whitespace and a addition operator.
For example, the specification A + B + C will be interpreted as a reaction \(\text{A} + \text{B} + \text{C}\) but the specification A+B + C will be interpreted as a reaction between one species A+B and another species C .
Internally, both the left- and right-hand sides are checked against background and plasma species. The program will perform a run-time exit if the species are not defined in these lists. Note that the right-hand side can additonal contain photon species.
Tip
Products not not defined as a species can be enclosed by parenthesis ( and ).
The reaction string generally expects all species on each side of the reaction to be defined.
It is, however, occasionally useful to include products which are not defined.
For example, one may with to include the reaction \(\text{e} + \text{O}_2^+ \rightarrow \text{O} +\text{O}\) but not actually track the species \(\text{O}\).
The reaction string can then be defined as the string e + O2+ -> but as it might be difficult for users to actually remember the full reaction, we may also enter it as e + O2+ -> (O) + (O).
Rate calculation
Reaction rates can be specified using multiple formats (constant, tabulated, function-based, etc).
To specify how the rate is computed, one must specify the type keyword in the JSON entry.
Constant
To use a constant reaction rate, the type specifier must be set to constant and the reaction rate must be specified through the value keyword.
A JSON specification is e.g.
"plasma reactions": [
{
"reaction": "e + N2 -> e + e + N2+" // Specify reaction
"type": "constant", // Reaction rate is constant
"value": 1.E-30 // Reaction rate
}
]
Function vs E/N
Some rates given as a function \(k = k\left(E/N\right)\) are supported, which are outlined below.
Important
In the below function-based rates, \(E/N\) indicates the electric field in units of Townsend.
function E/N exp A
This specification is equivalent to a fluid rate
The user must specify the constants \(c_i\) in the JSON file. An example specification is
"plasma reactions": [
{
"reaction": "A + B -> " // Example reaction string.
"type": "function E/N exp A", // Function based rate.
"c1": 1.0, // c1-coefficient
"c2": 1.0, // c2-coefficient
"c3": 1.0, // c3-coefficient
"c4": 1.0, // c4-coefficient
"c5": 1.0 // c5-coefficient
}
]
Temperature-dependent
Some rates can be given as functions \(k = k(T)\) where \(T\) is some temperature.
function T A
This specification is equivalent to a fluid rate
Mandatory input variables are \(c_1, c_2\), and the specification of the species corresponding to \(T_i\). This can correspond to one of the background species. An example specification is
"plasma reactions": [
{
"reaction": "A + B -> " // Example reaction string.
"type": "function T A", // Function based rate.
"c1": 1.0, // c1-coefficient
"c2": 1.0, // c2-coefficient
"T": "A" // Which species temperature
}
]
function TT A
This specification is equivalent to a fluid rate
Mandatory input variables are \(c_1, c_2\), and the specification of the species corresponding to \(T_1\) and \(T_2\). This can correspond to one of the background species. An example specification is
"plasma reactions": [
{
"reaction": "A + B -> " // Example reaction string.
"type": "function TT A", // Function based rate.
"c1": 1.0, // c1-coefficient
"c2": 1.0, // c2-coefficient
"T1": "A", // Which species temperature for T1
"T2": "B" // Which species temperature for T2
}
]
Townsend rates
Reaction rates can be specified to be proportional to \(\alpha\left|\mathbf{v}_i\right|\) where \(\alpha\) is the Townsend ionization coefficient and \(\left|\mathbf{v}_i\right|\) is the drift velocity for some species \(i\). This type of reaction is normally encountered when using simplified chemistry, e.g.
which is representative of the reaction \(S_i \rightarrow S_i + S_i\). To specify a Townsend rate constant, one can use the following:
alpha*vfor setting the rate constant proportional to the Townsend ionization rate.eta*vfor setting the rate constant proportional to the Townsend attachment rate.
One must also specify which species is associated with \(\left|\mathbf{v}\right|\) by specifying a species flag. A complete JSON specification is
{
"plasma reactions":
[
{
"reaction": "e -> e + e + M+", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e" // Species for v
}
]
}
Warning
When using the Townsend coefficients for computing the rates, one should normally not include any neutrals on the left hand side of the reaction.
The reason for this is that the Townsend coefficients \(\alpha\) and \(\eta\) already incorporate the neutral density.
By specifying e.g. a reaction string \(\text{e} + \text{N}_2 \rightarrow \text{e} + \text{e} + \text{N}_2^+\) together with the alpha*v or eta*v specifiers, one will end up multiplying in the neutral density twice, which will lead to an incorrect rate.
Tabulated vs E/N
Scaling
Reactions can be scaled by including a scale field in the JSON entry.
This will scale the reaction coefficient by the input factor, e.g. modify the rate constant as
where \(\nu\) is the scaling factor. This is useful when scaling reactions from different units, or for completely turning off some input reactions. An example JSON specification is
"plasma reactions":
[
{
"reaction": "e -> e + e + M+", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"scale": 0.0 // Scaling factor
}
]
Tip
If one turns off a reaction by setting scale to zero, the KMC algorithm will still use the reaction but no reactants/products are consumed/produced.
Efficiency
Reactions efficiencies can be modified in the same way as one do with the scale field, e.g. modify the rate constant as
where \(\nu\) is the reaction efficiency. Specifications of the efficiency can be achieved in the forms discussed below.
Constant efficiency
An example JSON specification for a constant efficiency is:
"plasma reactions":
[
{
"reaction": "e -> e + e + M+",
"efficiency": 0.5
}
]
Efficiency vs E/N
An example JSON specification for an efficiency computed versus \(E/N\) is
"plasma reactions":
[
{
"reaction": "e -> e + e + M+",
"efficiency vs E/N": "efficiencies.dat"
}
]
where the file efficiencies.dat must contain two-column data containing values of \(E/N\) along the first column and efficiencies along the second column.
An example file is e.g.
0 0
100 0.1
200 0.3
500 0.8
1000 1.0
This data is then internally convered to a uniformly spaced lookup table (see LookupTable).
Efficiency vs E
An example JSON specification for an efficiency computed versus \(E\) is
"plasma reactions":
[
{
"reaction": "e -> e + e + M+",
"efficiency vs E": "efficiencies.dat"
}
]
where the file efficiencies.dat must contain two-column data containing values of \(E\) along the first column and efficiencies along the second column.
This method follows the same as efficiency vs E/N where the data in the input file is put in a lookup table.
Quenching
Reaction quenching can be achieved in the following forms:
Pressure-based
The reaction rate can modified by a factor \(p_q/\left[p_q + p\left(\mathbf{x}\right)\right]\) where \(p_q\) is a quenching pressure and \(p\left(\mathbf{x}\right)\) is the gas pressure. This will modify the rate as
An example JSON specification is
"plasma reactions":
[
{
"reaction": "e + N2 -> e + N2 + Y", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"quenching pressure": 4000.0 // Quenching pressure
}
]
Important
The quenching pressure must be specified in units of Pascal.
Rate-based
The reaction rate can be modified by a factor \(k_r/\left[k_r + k_p + k_q\right]\). The intention behind this scaling is that reaction \(r\) occurs only if it is not predissociated (by rate \(k_p\)) or quenched (by rate \(k_q\)). Such processes can occur, for example, in excited molecules. This will modifiy the rate constant as
where the user will specifiy \(k_r\), \(k_p\), and \(k_q/N\).
The JSON specification must contain quenching rates, for example:
"plasma reactions":
[
{
"reaction": "e + N2 -> e + N2 + Y", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"quenching rates": { // Specify relevant rates
"kr": 1E9,
"kp": 1E9,
"kq/N": 0.0
}
}
]
Gradient correction
In LFA-based models it is frequently convenient to include energy-corrections to the ionization rate. In this case one modifies the rate as
To include this correction one may include a specifier gradient correction in the JSON entry, in which case one must also enter the species \(n_i\).
A JSON specification that includes this
"plasma reactions":
[
{
"reaction": "e -> e + e + M+", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"gradient correction": "e" // Specify gradient correction using species "e"
}
]
Understanding reaction rates
The JSON specification takes fluid rates as input to the KMC algorithm. Note that subsequent application of multiple scaling factors are multiplicative. For example, the specification
"plasma reactions":
[
{
"reaction": "e -> e + e + M+", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"gradient correction": "e" // Specify gradient correction using species "e"
"quenching pressure": 4000, // Quenching pressure
"efficiency": 0.1 // Reaction efficiency
}
]
is equivalent to the rate constant
Internally, conversions between the fluid rate \(k\) and the KMC rate \(c\) are made. The conversion factors depend on the reaction order, and ensure consistency between the KMC and fluid formulations.
Plotting rates
Reaction rates can be plotted by including an optional plot specifier.
If the specifier included, the reaction rates will be included in the HDF5 output.
The user can also include a description string which will be used when plotting the reaction.
The following two specifiers can be included:
plot(true/false) for plotting the reaction.description(string) for naming the reaction in the HDF5 output.
If the plot description is left out, the reaction string will be used as a variable name in the plot file. A JSON specification that includes these
"plasma reactions":
[
{
"reaction": "e -> e + e + M+", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"plot": true, // Plot this reaction
"description": "Townsend ionization" // Variable name in HDF5
}
]
Printing rates
ItoKMCJSON can print the plasma reaction rates (including photon-generating ones) to file.
This is done through an optional input argument print_rates in the input file (not the JSON specification).
The following options are supported:
ItoKMCJSON.print_rates = true
ItoKMCJSON.print_rates_minEN = 1
ItoKMCJSON.print_rates_maxEN = 1000
ItoKMCJSON.print_rates_num_points = 200
ItoKMCJSON.print_rates_spacing = exponential
ItoKMCJSON.print_rates_filename = fluid_rates.dat
ItoKMCJSON.print_rates_pos = 0 0 0
Setting ItoKMCJSON.print_rates to true in the input file will write all reaction rates as column data \(E/N, k(E/N)\).
Here, \(k\) indicates the fluid rate, so for a reaction \(A + B + C \xrightarrow{k}\ldots\) it will include the rate \(k\).
Reactions are ordered identical to the order of the reactions in the JSON specification.
This feature is mostly used for debugging or development efforts.
Photoionization
Photoionization reactions are specified by including a photoionization array of JSON entries that specify each reaction.
The left-hand side of the reaction must contain a photon species (plasma and background species can be present but will be ignored), and the right-hand side can consist of any species except other photon species.
Each entry must contain a reaction string specifier and optionally also an efficiency (which defaults to 1).
Assume that photons are generated through the reaction
"plasma reactions":
[
{
"reaction": "e + N2 -> e + N2 + Y", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"quenching pressure": 4000.0 // Quenching
"efficiency": 0.6 // Excitation events per ionization event
"scale": 0.1 // Photoionization events per absorbed photon
}
]
where we assume an excitation efficiency of \(0.6\) and photoionization efficiency of \(0.1\). An example photoionization specification is then
"photoionization":
[
{
"reaction": "Y + O2 -> e + O2+"
}
]
Note that the efficiency of this reaction is 1, i.e. the photoionization probabilities was pre-evaluated. As discussed in ItoKMCPhysics, we can perform late evaluation of the photoionization probability by specifically including a efficiency. In this case we modify the above into
"plasma reactions":
[
{
"reaction": "e + N2 -> e + N2 + Y", // Reaction string
"type": "alpha*v", // Rate is alpha*v
"species": "e", // Species for v,
"quenching pressure": 4000.0 // Quenching
"efficiency": 0.6 // Excitation events per ionization event
}
],
"photoionization":
[
{
"reaction": "Y + O2 -> e + O2+",
"efficiency": 0.1
},
{
"reaction": "Y + O2 -> (null)",
"efficiency": 0.9
}
]
where a null-absorption model has been added for the photon absorption. When multiple pathways are specified this way, and they have probabilities \(\xi_1, \xi_2, \ldots\), the reaction is stochastically determined from a discrete distribution with relative probabilities \(p_i = \xi_i/(\xi_1 + \xi_2+\ldots)\).
Important
The null-reaction model is not automatically added when using late evaluation of the photoionization probability.
Secondary emission
Secondary emission is supported for particles, including photons, but not for fluid species.
Specification of secondary emission is done separately on dielectrics and electrodes by specifying JSON entries electrode emission and dielectric emission.
The internal specification for these are identical, so all examples below use dielectric emission.
Secondary emission reactions are specified similarly to photoionization reactions.
However, the left-hand side must consist of a single species (photon or particle), while the right-hand side can consist of multiple species.
Wildcards @ are supported, as is the (null) species which enables null-reactions, as further discussed below.
An example JSON specification that enables secondary electron emission due to photons and ion impact is
"dielectric emission":
[
{
"reaction": "Y -> e",
"efficiency": 0.1
},
{
"reaction": "Y -> (null)",
"efficiency": 0.9
}
]
The reaction field specifies which reaction is triggered: The left hand side is the primary (outgoing) species and the left hand side must contain the secondary emissions.
The left-hand side can consist of either photons (e.g., Y) or particle (e.g., O2+) species.
The efficiency field specifies the efficiency of the reaction.
When multiple reactions are specified, we randomly sample the reaction according to a discrete distribution with probabilities
where \(\nu_i\) are the efficiencies.
Note that the efficiencies do not need to sum to one, and if only a single reaction is specified the efficiency specifier has not real effect.
The above reactions include the null reaction in order to ensure that the correct secondary emission probability is reached, where the (null) specifier implies that no secondary emission takes place.
In the above example the probability of secondary electron emission is 0.1, while the probability of a null-reaction (outgoing particle is absorbed without any associated emission) is 0.9.
The above example can be compressed by using a wildcard and an efficiencies array as follows:
"dielectric emission":
[
{
"reaction": "Y -> @",
"@": ["e", "(null)"],
"efficiencies": [1,9]
}
]
where for the sake of demonstration the efficiencies are set to 1 and 9 (rather than 0.1 and 0.9). This has no effect on the probabilities \(p(i)\) given above.
Warnings and caveats
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 inconsitencies in the user-provided chemistry, and provision of inconsistent \(\eta/N\) and attachment reaction rates are quite possible.
Example programs
Example programs that use the Îto-KMC module are given in
$DISCHARGE_HOME/Exec/Examples/ItoKMC/AirBasicfor a basic streamer discharge in atmospheric air.$DISCHARGE_HOME/Exec/Examples/ItoKMC/AirDBDfor a streamer discharge over a dielectric.