Reactions and their components#
What reactions (the abstraction) are#
As noted in the Introduction, we use reactions abstraction to represent the various physical collisional and reactive processes. Here we expand on those ideas and show the components of reactions as well as some examples.
We refer to reactions as any process that involves one or more ingoing particles (physical or otherwise) interacting with other particles in the simulation or fields (stored as ParticleDats), as well as one or more of the following:
The production of new particles in the simulation (with their own velocities, weights, and various internal states)
The modification of ingoing particle properties (weights, etc.)
Feedback on fields (such as particle/energy sources - assumed stored on the ingoing particle and projected onto the mesh separately)
In the abstract, a reaction is fully defined by:
Ingoing and outgoing particle IDs (these are notionally unique integer labels associated with particle species)
Any (per particle) data and associated calculation methods needed to apply the reaction - most notably the reaction rate
How the properties (NESO-Particles ParticleDats) of the parents and children are modified/generated
We make a distinction between linear and non-linear reactions in the particle sense. A linear reaction is any reaction where only one of the reactants is represented as a particle. There are no constraints on the number of reaction products in linear reactions. In contrast, a non-linear reaction is a reaction where two or more reactants are represented as particles in the simulation. An example of a non-linear reaction would be an elastic collision between two neutrals of the same species. There exist linearisation techniques for some of these reactions, so the initial focus of the library is linear reactions.
Linear Reaction structure#
The main components of reactions are the reaction data and reaction kernel objects. Their overall responsibilities are as follows:
Reaction data - calculate the per particle data required for the application of the reaction. This could be reaction rates, values randomly sampled from some distributions, etc.
Reaction kernels - define the properties of the products of the reaction (velocities, weights, internal states), as well as the feedback on fields and the parent particle
The key idea behind this separation of concern is the ability to separate the data and the physics, and allow the combination of different data calculation methods and different reaction physics. For example, the physics of an ionisation reaction is the same regardless of the reaction rate or the energy cost of the reaction, and the goal of flexibility in Reactions has lead to the data+kernel design.
The implementation of reactions, as well as reaction kernels will be covered in the developer guide, as it involves considerations of SYCL host and device types, as well as NESO-Particles ParticleLoop constructs.
Both data and kernels, in executing their responsibilities, access particle data, and use the property map system.
Reaction data and kernels are invoked in the two main ParticleLoop-containing methods in the linear reaction class, with the idea that data is calculated first, storing anything needed for the application of the kernels or for the global management of reaction application. Both loops are assumed to be invoked cell-wise, which allows for the reuse of various buffers.
Reaction data and the LinearReaction data loop#
Reaction data objects calculate a fixed number of components per particle. For example, data objects used to calculate reaction rates have a single component, while an object that is sampling velocities from a distribution might have two or more components. Each reaction needs at least one reaction data object - responsible for the calculation of the rate for the reaction. Further data calculation can be bundled using the DataCalculator container of reaction data.
Invoking the rate loop on a reaction object does the following:
Calculates the reaction rate and stores it in a local buffer used to apply the reaction using the kernels
Adds the calculated reaction rate to a total reaction rate
ParticleDat- used in the global management of reaction applicationCalculates and stores any
DataCalculatorresults for use by the kernels
Reaction kernels and the product loop#
With any data required to apply the reaction calculated and stored for some particles, the next step is to apply the reaction, which might involve feedback on fields and the ingoing particle, as well as some specification of product properties following the reaction. This is (semi-)independently specified by choosing a reaction kernel. The only requirement on the data that a kernel might have is that any required data exists, i.e. that the total dimensionality of data conforms with whatever the kernel requires. For example, if a kernel requires 2 sampled velocities, the DataCalculator must produce a total of 2 data values per perticle. Other than this requirement, data and kernels are independent.
Each kernel object consists of four kernel functions, in order to allow for extensibility. These are:
The scattering kernel - nominally specifies the velocities of the products
The weight kernel - nominally specifies the weights of the products
The transformation kernel - nominally specifies any complex internal state changes of the product
The feedback kernel - nominally specifies any feedback on the parent particle an on any fields, such as fluid sources
The above are applied in that order, and the product loop stores any products/children into a separate particle group in order to allow for any transformations before they are added into the group with the parents (see ReactionController documentation).
As noted above, both kernels and data specify their required particle properties using the property enum+map system, so that on-the-fly remaping of required variables is possible.
Putting a linear reaction together#
As noted above, to construct a linear reaction, we need to know the state IDs of the ingoing (parent/reactant) and outgoing (children/products) particles, as well as the data and kernels.
An example with the built-in charge-exchange kernels using fixed values for all of the data is given below. It demonstrates the pipeline needed to build a linear reaction object, as well as some of the method calls on the object relating to the two loops described above. More details on the individual data and kernel objects and their required properties will be presented below.
inline void linear_reaction_CX_example(ParticleGroupSharedPtr particle_group) {
auto particle_spec = particle_group->get_particle_spec();
// We take two species, one with internal state ID = 0 and one with ID = 1
// We shall treat the ID = 0 species as the projectile in a CX event
auto projectile_species = Species("ION", 1.2, 0.0, 0);
auto target_species = Species("ION2", 2.0, 0.0, 1);
// All reactions are applied to a subgroup. Here we intend to apply the
// reaction only to those particles with ID = 0, since those are the
// projectiles
//
// We use the following marking strategy, but avoid hardcoding the internal
// state ID by using the default map
auto prop_map = get_default_map();
auto mark_id_zero = make_marking_strategy<
ComparisonMarkerSingle<INT, EqualsComp>>(
Sym<INT>(
prop_map[default_properties.internal_state]), // this will result in
// "INTERNAL_STATE" when
// using the default map
projectile_species.get_id()); // projectile species internal state id = 0
// The resulting subgroup will have only particles with ID=0
auto input_subgroup = std::make_shared<ParticleSubGroup>(particle_group);
auto particle_sub_group = mark_id_zero->make_marker_subgroup(input_subgroup);
// For this example, we will use the FixedRateData reaction data class
// it simply sets the rate to a fixed number
auto rate_data =
FixedRateData(1.0); // the values used here are somewhat arbitrary,
// normally they would depend on the normalisation
// Since the CX kernels require ndim values per particle in additional
// reaction data we build a DataCalculator that will produce ndim values
//
// Here we choose ndim = 2, and the two data values per particle expected by
// the 2D CX kernel are the x and y values of the ion velocities sampled from
// the ion distribution
//
// For this example, we mimic a beam in 2D, by using two FixedRateData objects
// in the DataCalculator
//
// In this case, the beam would flow to the bottom right of the domain
auto vx_beam_data = FixedRateData(1.0);
auto vy_beam_data = FixedRateData(-1.0);
// DataCalculators are templated against their contents, so in this case
auto data_calculator =
DataCalculator<FixedRateData, FixedRateData>(vx_beam_data, vy_beam_data);
// Finally, the CXReactionKernels class only requires the dimensionality of
// the velocity space, as well as the two species
auto cx_kernel = CXReactionKernels<2>(
target_species, projectile_species,
prop_map // The property map used to remap any of the required properties
// in the kernel The CX kernel requires
// default_properties.weight, default_properties.velocity as well
// as some of the source properies (see the CXReactionKernels
// docs)
);
// We can now assemble the linear reaction using the base class constructor
auto cx_reaction = LinearReactionBase<
1, // The number of outgoing particles - CX will produce one neutral of
// the target species
FixedRateData, // The reaction data class to be used for the rate
// calculation
CXReactionKernels<2>, // The kernel class used
DataCalculator<FixedRateData, FixedRateData> // DataCalculator class used
>(particle_group
->sycl_target, // The reaction class needs access to the sycl_target
// of the group whose subgroups it's to be applied to
projectile_species.get_id(), // Ingoing partice state id
std::array<int, 1>{static_cast<int>(
target_species
.get_id())}, // State IDs of all the products - here just one
rate_data, // Reaction data used for the rate calculation
cx_kernel, // Reaction kernel object to be used
data_calculator,
prop_map // Map used for getting weight and total rate syms
);
// The following is normally handled by the ReactionController
//
// We will loop over all cells and generate products
int cell_count = particle_group->domain->mesh->get_cell_count();
auto product_group = std::make_shared<ParticleGroup>(
particle_group->domain, particle_spec, particle_group->sycl_target);
for (int i = 0; i < cell_count; i++) {
// This is the rate loop, here the reaction rates are calculated,
// they are added to a total reaction rate, and the DataCalculator
// performs any calculations in needs to
cx_reaction.calculate_rates(particle_sub_group, i, i + 1);
// For the product loop, the reaction needs to know the timestep (here
// arbitrarily set to 0.1) and the product group
//
// The timestep is used to calculate the total particle weight participating
// in the reaction as rate * timestep
cx_reaction.apply(particle_sub_group, i, i + 1, 0.1,
product_group);
}
return;
}
Reaction data types#
Reactions offers a number of built-in data types. These will be covered here in the following format:
Dimensionality - the number of data values produced by the data object per particle
Required properties - all reaction data objects provided by Reactions use the default properties enum and their required properties will be listed (both simple and species properties, where applicable)
Details - any explantion of the calculations done by the data object, e.g. formulae, restrictions, etc.
Example - where the constructor of the object is non-trivial an example of how to construct it is given
Fixed rate data#
Dimensionality: 1
Required properties: none
Details: The rate is simply set to a fixed value
Example: See the example in the previous section
Fixed rate coefficient data#
Dimensionality: 1
Required properties: Simple props: weight; Species props: none
Details: Given a coefficient \(K\), and a particle weight \(w\), the rate is given as \(Kw\), with \(K\) being fixed.
Example:
inline void fixed_rate_coeff_example() {
// In case we would like to remap the used weight Sym
auto used_map = get_default_map();
auto fixed_coeff_rate =
FixedCoefficientData(5.0, // K - fixed rate coefficient
used_map);
return;
}
AMJUEL 1D rate fit#
Dimensionality: 1
Required properties: Simple props: fluid_density, fluid_temperature, weight; Species props: none
Details: Uses the following fit for the rate coefficient from AMJUEL
\[K=\ln\langle\sigma v\rangle = \sum_{n=0}^N b_n (\ln T)^n\]
where the number of coefficients \(N\) and the coefficients \(b_n\) are set on construction. The final output rate is given as \(nKw\), where \(n\) here is the fluid density and \(w\) is the particle weight. All normalisation is set in the constructor (see the example). The rate is assumed to evolve some quantity \(q\), and requires the knowledge of the normalisation of that quantity. For example, if evolving the weight it should be left at 1.0 while if evolving a background energy field (e.g. providing an energy source) it would require the normalisation of the energy density (see below for the assumed normalisation in case of built-in kernels).
Example:
inline void amjuel_1d_example() {
// In case we wish to remap the default weight, fluid_temperature,
// fluid_density
auto used_map = get_default_map();
auto coeffs = std::array<REAL, 3>{1.0, 1.0, 1.0}; // b_n coefficients
// We pass the number of fit coefficients to the constructor as a template
// parameter
auto amjuel_data =
AMJUEL1DData<3>(1.0, // The normalisation of the evolved quantity
// (density, energy, particle weight, etc.)
1e19, // Normalisation of density in m^{-3}
1.0, // Temperature normalistion in eV
1e-8, // Time normalisation in seconds
coeffs, // fit coefficients
used_map); // Optional property map
return;
}
AMJUEL 2D rate fit (n,T)#
Dimensionality: 1
Required properties: Simple props: fluid_density, fluid_temperature, weight; Species props: none
Details: Uses the following fit for the rate coefficient from AMJUEL
\[K=\ln\langle\sigma v\rangle = \sum_{n=0}^N \sum_{m=0}^M \alpha_{n,m}(\ln \tilde{n})^m (\ln T)^n\]
where the numbers of coefficients \(N\) and \(M\), and the coefficients \(\alpha_{n,m}\) are set on construction. \(\tilde{n}\) is density rescaled to \(10^{14} m^{-3}\). Density dependence is dropped below \(\tilde{n}=1\), and only the \(m=0\) coefficients are used (this is the Coronal approximation). The LTE limit is not implemented yet (for densities above \(10^{22} m^{-3}\)). The final output rate is given as \(nKw\), where \(n\) here is the fluid density and \(w\) is the particle weight. Normalisation as in the 1D fit case.
Example:
inline void amjuel_2d_example() {
// In case we wish to remap the default weight, fluid_temperature,
// fluid_density
auto used_map = get_default_map();
// alpha coefficients (the inner array size is the number of density
// coefficients, the outer is temperature)
auto coeffs = std::array<std::array<REAL, 2>, 2>{
std::array<REAL, 2>{1.0, 0.02}, std::array<REAL, 2>{0.01, 0.02}};
// We pass the number of fit coefficients to the constructor as a template
// parameter
auto amjuel_data =
AMJUEL2DData<2, 2>(1.0, // The normalisation of the evolved quantity
// (density, energy, particle weight, etc.)
1e19, // Normalisation of density in m^{-3}
1.0, // Temperature normalistion in eV
1e-8, // Time normalisation in seconds
coeffs, // fit coefficients
used_map); // Optional property map
return;
}
AMJUEL 2D rate fit (E,T)#
Dimensionality: 1
Required properties: Simple props: fluid_density, fluid_temperature, fluid_flow_speed, weight, velocity; Species props: none
Details: Uses the following fit for the rate coefficient from AMJUEL section H.3
\[K=\ln\langle\sigma v\rangle = \sum_{n=0}^N \sum_{m=0}^M \alpha_{n,m}(\ln E)^m (\ln T)^n\]
where the numbers of coefficients \(N\) and \(M\), and the coefficients \(\alpha_{n,m}\) are set on construction. The neutral energy \(E\) is relative to the fluid flow speed. The final output rate is given as \(nKw\), where \(n\) here is the fluid density and \(w\) is the particle weight. Normalisation as in the 1D fit case, with the added normalisation of the velocity and the requirement for the neutral energy to be specified in amus.
Example:
inline void amjuel_2d_H3_example() {
// In case we wish to remap the default weight, fluid_temperature,
// fluid_density, fluid_flow_speed, or particle velocity
auto used_map = get_default_map();
// alpha coefficients (the inner array size is the number of neutral energy
// coefficients, the outer is temperature)
auto coeffs = std::array<std::array<REAL, 2>, 2>{
std::array<REAL, 2>{1.0, 0.02}, std::array<REAL, 2>{0.01, 0.02}};
// We pass the number of fit coefficients to the constructor as a template
// parameter
auto amjuel_data =
AMJUEL2DDataH3<2, 2>(1.0, // The normalisation of the evolved quantity
// (density, energy, particle weight, etc.)
1e19, // Normalisation of density in m^{-3}
1.0, // Temperature normalistion in eV
1e-8, // Time normalisation in seconds
1e6, // Velocity normalisation in m/s
1.0, // Neutral mass in amus
coeffs, // fit coefficients
used_map); // Optional property map
return;
}
Filtered Maxwellian sampler#
Dimensionality: variable - corresponds to fluid flow field dimensionality
Required properties: Simple props: fluid_temperature, fluid_flow_speed, velocity; Species props: none
Details: Produces velocity components sampled from a drifting Maxwellian at the local fluid temperature and with local mean fluid flow. Optionally filters the sampled velocities based on an interaction cross section using a rejection method, effectively sampling from
\[f(\vec{v}) \propto \sigma(v_{rel}) f_M(\vec{v},\vec{u},T)\]
where \(\sigma(v_{rel})\) is the interaction cross-section evaluated at the relative velocity \(|\vec{v}-\vec{u}|\), and \(\vec{u}\) and \(T\) are the fluid flow speed and temperature, respectively. By default, the cross-section is assumed constant, which just leads to sampling from a drifting Maxwellian. See below for cross-section objects.
Example:
inline void maxwellian_sampler_example() {
// In case we wish to remap the fluid_temperature, fluid_flow_speed, or
// velocity
auto used_map = get_default_map();
// Default cross-section object - results in sampling from an unfiltered
// drifting Maxwellian
auto default_cs = ConstantRateCrossSection(1.0);
// The sampler needs a NESO-Particles rng_kernel
// In general, it will need the atomic block kernel (see NESO-Particles
// documentation) This is because rejection sampling has an a priori unknown
// number of samples
//
// Here we use an arbitrary lambda, but this should in general be a standard
// uniform distribution
auto rng_lambda = [&]() -> REAL { return 0.5; };
auto rng_kernel = host_atomic_block_kernel_rng<REAL>(rng_lambda, 1000);
// The sampler is templated against velocity space dimensionality - here 2D
auto sampler_data = FilteredMaxwellianSampler<2>(
1.0, // This is the ratio between kT_0 and mv_0^2,
// where T_0 is the temperature normalisation,
// m is the mass of the species whose distribution is sampled,
// and v_0 is the velocity normalisation
default_cs, // Optional cross-section object - here explicitly defaulted
rng_kernel, // RNG kernel used to perform Box-Muller sampling and
// rejection sampling
used_map); // Optional property map
return;
}
Cross-section objects#
Currently, cross-section objects are restricted to being used by the above sampler. For that purpose, they can be evaluated at a given relative velocity, and have an associated maximum \(\sigma v_{rel}\), used for rejection sampling.
Constant rate cross-section#
This is the simplest cross section, with \(\sigma v_{rel} = c\). It always accepts all samples. See the above sampler example.
AMJUEL H.1 cross-section fit#
These are single parameter fits in lab energy from AMJUEL of the form
where we convert to the centre-of-mass energy. The coefficients \(a_n\) can be given for asymptotic values of the energy, as well, both high or low.
inline void amjuel_h1_cs_example() {
// These are example values
auto coeffs = std::array<REAL, 3>{1.0, 1.0, 1.0}; // a_n coefficients
auto l_coeffs =
std::array<REAL, 3>{1.0, 1.0, 1.0}; // left asymptote a_n coefficients
auto r_coeffs =
std::array<REAL, 3>{1.0, 1.0, 1.0}; // right asymptote a_n coefficients
REAL E_lab_max = 1e3; // energy value after which the r_coeffs are used
REAL E_lab_min = 1e-1; // energy value after which the l_coeffs are used
REAL reduced_mass_amu = 1.0;
auto cs = AMJUELFitCrossSection<
3, // Dimensionality of bulk energy fit
3, // Dimensionality of left asymptote fit - no asymptote if 0
3 // Dimensionality of right asymptote fit - no asymptote if 0
>(1e6, // velocity normalisation
1e-4, // cross-section normalisation in m^2
reduced_mass_amu,
coeffs, // Bulk fit coefficients
l_coeffs, // Left asymptote coefficients - set to std::array<REAL,0>{}
// if no left asymptote
r_coeffs, // Right asymptote coefficents - set to std::array<REAL,0>{}
// if no right asymptote
E_lab_min, // Left asymptote energy threshold - ignored if l_coeffs of
// size 0
E_lab_max, // Right asymptote energy threshold - ignored if r_coeffs of
// size 0
1e5); // this->Maximum expected energy - used to evaluate the maximum
// value of the cross-section for rejection sampling - after this
// value, the cross-section decays as 1/v_rel
return;
}
Reaction kernel types#
VANTAGE-Reactions offers several built-in reaction kernels. These are presented in the following format:
Overview - general description, number of products, required
DataCalculatortotal dimensionality, etc.Required properties - the required properties from the default properties enum (as for reaction data) - here split into parent and descendant
Scattering kernel - if there are any products, how their velocities are calculated
Weight kernel - if there are any products, how weight is distributed amongst them
Transformation kernel - if there are any products, how aspects of their internal states are set
Feedback kernel - determines how the parent weight is affected, as well as how the various source
ParticleDatvalues on the parent are setExample - example of constructing the kernel
Kernels that produce products have a set of specified descendant particle required properties. These are usually the particle velocities and weights, and are modified by the kernel. All other properties are copied from the parent, so care should be taken if some of these need zeroing (sources, etc.).
NOTE: Reactions assumes all sources are ParticleDat objects on particles. All pre-built kernels also assume that the sources are not rates, i.e. that the user will divide them by the timestep
lenght to get the rate after applying reactions. This is so that different length timesteps could be used for different reactions, or so that operator splitting can be done without worrying about the individual steps.
Base ionisation kernels#
Overview: These are general ionisation kernels with the fewest possible assumptions. Since ionisation is an absorption process, there are no descendant particles. This implementation allows for different electron, projectile, and target species, i.e. it represents projectile-impact target ionisation. It expects at least one
DataCalculatorvalue, representing the energy loss rate of the projectile species in the process. Optionally, a momentum loss rate can be included, with the momentum being transferred to the target species. Electron momentum is assumed negligible. NOTE: The units of the energy and momentum sources are tied to the velocity normalisation via the weight and amu - e.g. the energy source normalisation is assumed to be \(w_0m_0 v_0^2\), where \(m_0\) is the amu, \(v_0\) is the velocity normalisation, and \(w_0\) represents the weight normalisation (for example a number of particles associated with unit weight)Required properties:
Parent: Simple props: weight, velocity; Species props: source_density, source_energy, source_momentum
Descendant: N/A
Scattering kernel: N/A
Weight kernel: N/A
Transformation kernel: N/A
Feedback kernel: The total weight participating in the reaction is removed from the parent particle. The first
DataCalculatorvalue is used as the energy rate, and if a momentum rate is marked as set, the second value is interpreted as that. Density sources for the target and electron species are set to that same weight value. The target momentum source always includes the momentum of the ionised neutral, regardless of the presence of a momentum kernel.Example:
inline void ionisation_kernels_example() {
// In case we would like to remap the used Syms
auto used_map = get_default_map();
auto electron_species = Species("ELECTRON", // name
5.5e-4, // electron mass in amu
-1.0 // charge
);
auto target_species = Species(
"ION", 2.0, 0.0, 1); // This is the target species, i.e. the ID
// corresponds to the neutral being ionised and the
// species name corresponds to the ion fluid
auto ion_kernels = IoniseReactionKernels<
2, // velocity dat dimensionality
2, // momentum source dat dimensionality (defaults to the velocity dat
// dimensionality)
false // set to true if there is an expected momentum source rate data in
// the data calculator //
>(target_species, // target species - neutral species and resulting ion
// species
electron_species, // electron species - in general used only to store a
// density source
electron_species, // projectile species - energy and momentum sources
// (so here electrons will have an energy and a
// particle source contribution)
used_map // Optional map for remapping property names
);
return;
}
Base charge-exchange kernels#
Overview: These kernels perform direct charge-exchange with a pre-sampled ion. The ion velocities are assumed to be set in the accompanying
DataCalculatorobject. As such, this kernel is not in charge of the sampling process (use, for example, theFilteredMaxwellianSampler). These kernels assume one reaction product, which is the resulting charge-exchanged neutral particle. NOTE: Energy and momentum source normalisation are the same here as in the ionisation kernels.Required properties:
Parent: Simple props: weight, velocity; Species props: source_density, source_energy, source_momentum
Descendant: Simple props: weight, velocity, internal_state; Species props: N/A
Scattering kernel: Sets the product velocities to the pre-calculated velocities from the
DataCalculator.Weight kernel: The product gets the full weight that participated in the reaction
Transformation kernel: The products internal_state is set to the correct species ID
Feedback kernel: The total weight participating in the reaction is removed from the parent particle. The energy and momentum sources are computed from the participating particles’ velocities (the parent and the sample ion)
Example: See the above example on putting together a linear reaction for an example of a CX kernel being constructed and used
Base recombination kernels#
Overview: These kernels allow for implementing recombination using pseudo-particles (also referred to as markers). The self-consistent calculation of the rates used by the kernels is left to the users, as it depends on the mesh properties (i.e. the mapping of ion densities to the marker weights). Like the ionisation kernels, assumes that the first value calculated by the
DataCalculatoris the electron energy loss rate (not including the potential energy). Similarly to the CX kernel above, this kernel assumes pre-sampled ion velocities set by aDataCalculator(after the energy loss rate). Recombination produces a single product, and does not modify the weights of the parents/marker particles. NOTE: Energy and momentum source normalisation are the same here as in the previous two kernels.Required properties:
Parent: Simple props: weight; Species props: source_density, source_energy, source_momentum
Descendant: Simple props: weight, velocitym internal_state; Species props: N/A
Scattering kernel: Sets the product velocities to the pre-calculated velocities from the
DataCalculator(excluding the first entry)Weight kernel: The product gets the full weight that participated in the reaction
Transformation kernel: The products internal_state is set to the correct species ID
Feedback kernel: Weight is not removed from the parent, but the particle sources is updated as if it were. The energy and momentum source of the target species (the ions) are computed from the sampled velocities. The projectile species (electron) momentum source is assumed to be negligible, while the energy cost is calculated using the energy loss rate \(K_E\) and the normalised (to \(m_0 v_0^2\)) ionisation potential energy \(\epsilon_i\) as \(- K_E \Delta t - \epsilon_i \Delta w\), where the timestep and weight participating in reaction are set self-consistently.
Example:
inline void recombination_kernels_example() {
// In case we would like to remap the used Syms
auto used_map = get_default_map();
auto electron_species = Species("ELECTRON", // name
5.5e-4, // electron mass in amu
-1.0 // charge
);
auto target_species = Species("ION", 2.0, 0.0, -1); // This is the target species, i.e. the ID corresponds to the
// marker species and the species name corresponds to the ion fluid
auto reaction_energy_rate = FixedRateData(1.0); // Energy rate for
// projectile energy loss
// per recombination event
// Assuming the ions are a beam with given x and y speeds
auto ion_vel_x = FixedRateData(1.0);
auto ion_vel_y = FixedRateData(2.0);
auto normalised_potential_energy = 1.0; // Set for convenience, otherwise should
// be normalised to m_0 * v_0^2
// where m_0 is the mass normalisation (usually amu)
// and v_0 is the velocity normalisation
// Used data calculator
auto data_calculator =
DataCalculator<decltype(reaction_energy_rate),
decltype(ion_vel_x),
decltype(ion_vel_y)>(reaction_energy_rate,
ion_vel_x,
ion_vel_y);
auto recombination_kernels = RecombReactionKernels<2, // velocity dat dimensionality
2 // momentum source dat dimensionality (defaults to the velocity dat dimensionality)
>(
target_species, // target species - marker species corresponding to the ions
electron_species, // projectile species
normalised_potential_energy, // normalised ionisation potential to be included in projectile energy source
used_map // Optional map for remapping property names
);
return;
}
Pre-built reactions#
VANTAGE-Reactions offers pre-built reaction classes that bundle commonly used options together. It should be noted that these can be completely reproduced by users from the base reaction class and data and kernels.
Electron-impact ionisation#
Given that the most commonly treated class of ionisation reactions is electron-impact ionisation, the library offers a streamlined way of constructing electron-impact ionisation reactions. See below for an example of such a reaction.
inline void electron_impact_ion_example(ParticleGroupSharedPtr particle_group) {
auto used_map = get_default_map();
auto electron_species = Species("ELECTRON", // name
5.5e-4, // electron mass in amu
-1.0 // charge
);
auto target_species = Species(
"ION", 2.0, 0.0, 1); // This is the target species, i.e. the ID
// corresponds to the neutral being ionised and the
// species name corresponds to the ion fluid
auto test_data = FixedRateData(1.0); // Example fixed rate data
auto ionise_reaction = ElectronImpactIonisation<
FixedRateData, // Reaction data class used for the reaction rate
FixedRateData // Reaction data class used for the energy rate
>(particle_group->sycl_target, // Reactions need access to the SYCL target
test_data, // Reaction rate data
test_data, // Energy rate data
target_species, // Ionisation target species
electron_species, // Electron species object (projectile)
used_map // Weight and total reaction rate remapping - here the default
// map
);
return;
}
Recombination#
Like electron-impact ionisation, a recombination reaction can be constructed directly without using LinearReactionBase:
inline void recombination_reaction_example(ParticleGroupSharedPtr particle_group) {
// In case we would like to remap the used Syms
auto used_map = get_default_map();
auto electron_species = Species("ELECTRON",
5.5e-4,
-1.0
);
auto marker_species = Species("ION", 2.0, 0.0, -1); // This is the target species, i.e. the ID corresponds to the
// marker species and the species name corresponds to the ion fluid
auto neutral_species = Species("ION", 2.0, 0.0, 0); // This is the recombined species, i.e. the ID corresponds to the
// neutral species and the species name corresponds to the ion fluid
auto reaction_rate = FixedRateData(1.0); // Reaction rate
// NOTE: this would in reality
// have to account for the background
// electrons and ions
// Following the recombination_kernels_example
auto reaction_energy_rate = FixedRateData(1.0);
auto ion_vel_x = FixedRateData(1.0);
auto ion_vel_y = FixedRateData(2.0);
auto normalised_potential_energy = 1.0;
auto data_calculator =
DataCalculator<decltype(reaction_energy_rate),
decltype(ion_vel_x),
decltype(ion_vel_y)>(reaction_energy_rate,
ion_vel_x,
ion_vel_y);
auto recombination_reaction =
Recombination<decltype(reaction_rate), decltype(data_calculator)>(
particle_group->sycl_target, // Reactions need to know the used SYCL target
reaction_rate, // Reaction rate data object
data_calculator,// Data calculator containing the energy rate and the velocity sampling
marker_species, // Marker pseudo-particle species - ingoing particle representing ions
electron_species, // Electron species - will have energy loss set by given rate and potential energy
neutral_species, // Species into which the ions recombine - product particle
normalised_potential_energy // Normalised ionisation potential energy (to m_0 v_0^2)
);
return;
}
NOTE: To properly use recombination, especially with built-in AMJUEAL reaction data, care must be taken that the marker weights are updated in a way consistent with the background ion densities.