Reactions and their components#
The reaction abstraction#
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 application
Reaction kernels and the product loop#
With any rate 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. Note that data calculated by the DataCalculator is calculated at the time of application of the reaction. This is to avoid unnecessary computations (for example when randomly selecting which particles have reacted we do not want to calculate all of the data for particles that did not react).
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.
void linear_reaction_CX_example(NP::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
auto prop_map = get_default_map();
auto spec_id =
projectile_species.get_id(); // projectile species internal state id = 0
// The resulting subgroup will have only particles with ID=0
// using the default properties name for the internal state
auto input_subgroup = std::make_shared<NP::ParticleSubGroup>(particle_group);
auto particle_subgroup = particle_sub_group(
particle_group, [=](auto id) { return id[0] == spec_id; },
NP::Access::read(
NP::Sym<INT>(prop_map[default_properties.internal_state])));
// 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<NP::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_subgroup, 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_subgroup, i, i + 1, 0.1, product_group);
}
return;
}
Reaction data types#
As noted above, reaction data objects are used both for reaction rate calculations as well as any other data that might be required for the desired application of reactions.
Broadly, the data can be split up into the following groups:
Rate-specific data - data originally intended to be used for various reaction rate calculations, including other size 1 data
Multi-dimensional data - data used for such things as sampling velocities from a distribution
Composite data - data objects that represent compositions of other data objects. These allow for such things as pipeline construction
Surface reaction data - objects designed for calculating various surface interaction data, such as post-reflection velocities
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 (if the data is a composite and requires a specific dimensionality of input data this will be noted here).
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 explanation 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.
Rate-specific and size-1 data#
The following data objects all return size 1 data, and the majority of them are meant to be used as reaction rate data.
Fixed rate data#
Dimensionality: 1
Required properties: none
Details: The rate is simply set to a fixed value \(K\), so that the weight evolution equation (assuming deterministic evolution) is:
\[\frac{dw}{dt} = -K\]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. Leads to the following deterministic weight evolution equation
\[\frac{dw}{dt} = -kw\]Example:
void fixed_rate_coeff_example() {
// In case we would like to remap the used weight NP::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). When used as the deterministic reaction rate (evolving weight), leads to
\[\frac{dw}{dt} = -nkw\]Example:
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 and effective deterministic evolution equation as in the 1D fit case.
Example:
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. Deterministic evolution equation as in the 1D fit case.
Example:
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;
}
Arrhenius data#
Dimensionality: 1
Required properties: Simple props: weight, fluid_temperature; Species props: none
Details: Given two coefficients \(a\) and \(b\), returns an Arrhenius form rate \(a T^b w\), with \(T\) being a temperature, and \(w\) being the particle weight. NOTE: for many reactions this will need to be multiplied by one or more densities - see below entries on composite data for how one can do that. Leads to the following deterministic weight evolution equation (assuming no additional density multiplication)
\[\frac{dw}{dt} = -aT^bw\]Example:
void arrhenius_example() {
// In case we wish to remap the default weight, fluid_temperature
auto used_map = get_default_map();
// The following will calculate the rate as 2 * T^3 * w
auto arrhenius_data =
ArrheniusData(2.0, // The a coefficient in the Arrhenius formula
3.0, // The b coefficient in the Arrhenius formula
used_map); // Optional property map
return;
}
Sampler data#
Dimensionality: 1
Required properties: Simple props: none; Species props: none
Details: Returns a sample from the contained random number generator kernel. Not meant to be used as a rate data object. Instead, it was implemented as a way of allowing random sample input to composite data objects.
Example:
void sampler_example() {
// In case we wish to remap the default panic flag - used in case the sampling
// fails
auto used_map = get_default_map();
// The sampler needs a NESO-Particles rng_kernel
//
// Here we use an arbitrary lambda, but this can be anything
auto rng_lambda = [&]() -> REAL { return 0.5; };
auto rng_kernel = NP::host_atomic_block_kernel_rng<REAL>(rng_lambda, 1000);
// The following will just sample from the above kernel
auto sampler_data = SamplerData(rng_kernel,
used_map); // Optional property map
return;
}
Multi-dimensional data#
The following data objects allow for generating multidimensional data, with the most common use case being velocity generation (e.g. post-scattering values), or calculating inputs into composite data objects.
Fixed array data#
Dimensionality: any
Required properties: Simple props: none; Species props: none
Details: Returns a fixed array of values. Useful for providing fixed inputs to composite data objects.
Example:
void fixed_array_data_example() {
auto data_array = std::array<REAL, 3>{1.0, 2.0, 3.0};
// The following will just return the above array
// NOTE: templated on array size so can return any length array
auto fixed_array = FixedArrayData(data_array);
return;
}
Array lookup data#
Dimensionality: any
Required properties: Simple props: custom key (specified by
Sym); Species props: noneDetails: Uses an integer property on the particle as a key for a map, returning arrays based on the value of the property, with the option to specify a default return value if the key isn’t in the map.
Example:
void array_lookup_example(NP::ParticleGroupSharedPtr particle_group) {
// The arrays can be any size
// The default array is what is returned if the key value isn't
// found
auto default_array = std::array<REAL, 1>{1.0};
std::map<int, std::array<REAL, 1>> lookup_table_map;
lookup_table_map[0] = std::array<REAL, 1>{2.0};
lookup_table_map[1] = std::array<REAL, 1>{3.0};
// The following data will check the first component of the particle's
// INTERNAL_STATE value, use it as the lookup key in the above map, and
// return the corresponding map value if found, otherwise returning
// the default array
auto array_lookup_data =
ArrayLookupData<1>(NP::Sym<INT>("INTERNAL_STATE"), // The key NP::Sym to
// be used for lookup
0, // The component of the
// NP::ParticleDat referred
// to by the key NP::Sym to
// be used for the key
// value
lookup_table_map, default_array,
particle_group->sycl_target); // A sycl target is
// needed to store the
// on-device lookup table
return;
}
Extractor data#
Dimensionality: variable (up to the dimensionality of the extracted property)
Required properties: Simple props: custom key (specified by
Sym); Species props: noneDetails: Returns the first N components of a given REAL particle property. Useful for providing inputs into composite data objects.
Example:
void extractor_example() {
// Extract the first 2 components (template arg) of the particle POSITION
auto extracted_data = ExtractorData<2>(NP::Sym<REAL>("POSITION"));
// Alternatively
auto extracted_data_quick = extract<2>("POSITION");
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:
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 = NP::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;
}
One-way Maxwellian flux sampler#
Dimensionality: 3
Required properties: Simple props: fluid_temperature, fluid_flow_speed, surface_basis_e1, surface_basis_e2, surface_basis_pi; Species props: none
Details: Generate particle velocities sampled from a one-way Maxwellian distribution defined with respect to a surface. The two components tangential to the surface are sampled from a drifting Maxwellian, while the component along
surface_basis_piis sampled from the one way Maxwellian, referred to as truncated drifting Maxwellian in section “1.5 Recycling surface sources” of the EIRENE documentation. Thus, with \(v_\pi \geq 0\), the normal component is sampled from\[f(v_\pi) \propto v_\pi \exp\left[-\frac{(v_\pi-v_{d,\pi})^2}{2\sigma^2}\right], v_\pi \gt 0\]where \(v_{d,\pi}\) is the drift/fluid flow velocity component along
surface_basis_piand \(\sigma^2 = T\,\mathtt{norm\_ratio}\). The three surface-basis vectors are assumed to be orthonormal. Sampling of the normal component uses rejection sampling following Makkonen, Airila, and Kurki-Suonio (2015), doi:10.1088/0031-8949/90/1/015204.
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.
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;
}
Composite data#
The following data objects all have in common that they contain other data objects, and that they perform all operations on device, i.e. without saving intermediate states in a buffer such as the one used by the DataCalculator.
Since the purely composite data objects in this section require particle properties only through their contained objects, that section of the description will be omitted for brevity.
Concatenator data#
Dimensionality: sum of the dimensions of the contained objects
Details: Takes as aruments any number of data objects, evaluates their outputs, and concatenates the result, similarly to the effects of the
DataCalculator, except fully on device.Example:
void concatenator_example() {
auto data_1 = FixedRateData(1.0);
auto data_2 = FixedRateData(2.0);
// Will return an array with the results of the contained
// objects concatenated - [1,2] in this case
auto concatenated_data = ConcatenatorData(data_1, data_2);
return;
}
Pipeline data#
Dimensionality: the output dimensionality of the final step in the pipeline
Details: Takes as arguments any number of data objects, and passes the outputs from left to right. The data objects must have compatible input/output dimensions (see example). All calculations are performed on device.
Example:
void pipeline_example() {
// Here we will pipe particle velocities into a specular reflection data
// object See the documentation of the specular reflection for more details
// In case we wish to remap the default boundary normal for the specular
// reflection
auto used_map = get_default_map();
auto velocity_data = extract<2>("VELOCITY");
auto specular_reflection = SpecularReflectionData<2>(used_map);
// The pipeline object will first evaluate the velocity extractor
// and will then pipe this into the specular reflection,
// resulting in a specularly reflected velocity (assuming that the surface
// normal is correctly set)
//
// In general, this allows for more flexibility, as we might transform the
// velocity data somehow before passing it to the reflection data
auto pipeline = PipelineData(velocity_data, specular_reflection);
// Alternative syntax
auto pipeline_quick = pipe(velocity_data, specular_reflection);
return;
}
Array transform data#
The following data objects perform unary or binary transformation on arrays, allowing for composition of simple operations.
Unary array transform data#
Dimensionality: varies based on the transformation applied - input dimension generally non-zero, so these objects must be part of a pipeline
Details: Allow taking the result of another data object and applying a unary transformation on it, returning the result. A number of transformations are implemented - see example below
Example:
void unary_array_transform_examples() {
// Each unary array transform data object applies a transformation to
// an input array
// This is just some dummy data for the example
auto position_data = extract<2>("POSITION");
// ------------------------
// PolynomialArrayTransform
// ------------------------
// The following transformation takes an array of polynomial
// coefficients and applies them to the input.
//
// The two template arguments are the expected input dimension and the
// polynomial order (1 less than the dimension of the coefficients)
//
// In this case, the elementwise polynomial will be 2x+1
auto linear_poly =
PolynomialArrayTransform<2, 1>(std::array<REAL, 2>{1.0, 2.0});
// The transform can be wrapped into a reaction data object
auto linear_poly_data = UnaryArrayTransformData(linear_poly);
// And the object can have the input data piped into it
auto pipe_poly = pipe(position_data, linear_poly_data);
// ------------------------
// ScalerArrayTransform
// ------------------------
// The following transform just scales the input elementwise by a number
// The template argument is the expected input dimension
auto scaler = ScalerArrayTransform<2>(2.0);
auto scaler_data = UnaryArrayTransformData(scaler);
// or more succinctly
auto scaler_data_quick = scale_by<2>(2.0);
// And either of the above can then be used in a pipeline.
// ------------------------
// UnaryProjectArrayTransform and UnaryProjectNormalArrayTransform
// ------------------------
//
// These transforms take in a constant direction, and either project the input
// onto that direction, or onto the plane normal to it. Note that if the
// direction vector isn't a unit vector the projection will be scaled by
// the square of its magnitude
//
// The transform expects the same dimensionality of input data as that of the
// direction vector
auto dir = std::array<REAL, 2>{1.0, 0.0};
// Projects onto dir
auto project = UnaryProjectArrayTransform(dir);
// Project onto the plane normal to dir (in this case [0,1])
auto project_normal = UnaryProjectNormalArrayTransform(dir);
// The above can then be wrapped in UnaryArrayTransformData and used in a
// pipeline
return;
}
Binary array transform data#
Dimensionality: varies based on the transformation applied
Details: Allows taking two reaction data objects and applying a binary transformation on their result. A number of transformations are implemented - see example below
Example:
void binary_array_transform_examples() {
// Each binary array transform data object applies a transformation to
// the results of two other data objects
// Unlike the unary data, binary array transforms do not by definition require
// being part of a pipeline
// 1D and 2D reaction data objects for the examples
auto position_data_x = extract<1>("POSITION");
auto position_data_xy = extract<2>("POSITION");
auto velocity_data_xy = extract<2>("VELOCITY");
// ------------------------
// Doing arithmetic with ReactionData objects
// ------------------------
// Binary arithmetic with ReactionData objects is enabled through
// BinaryArrayTransforms
// The following produces a ReactionData object that returns the product of
// the results of the two RHS objects
auto pos_squared = position_data_xy * position_data_xy;
// If one of the ReactionData objects is 1D, the arithmetic implementation
// supports broadcasting onto higher dimensionality objects
// For example, the following ReactionData object would produce a 2D reaction
// data object that returns [pos_x^2,pos_x*pos_y]
auto pos_times_pos_x = position_data_xy * position_data_x;
// All 4 basic arithmetic operations [+,-,*,/] are supported between
// ReactionData objects.
// ------------------------
// BinaryDotArrayTransform
// ------------------------
// Takes the dot product of the results of the contained data objects
auto binary_dot_transform = BinaryDotArrayTransform<2>();
// The following data object, when used, will calculate the dot product
// of the position and velocity vectors (evaluated using the above extractors)
auto dot_transform_data = BinaryArrayTransformData(
binary_dot_transform, position_data_xy, velocity_data_xy);
// Or, more succinctly
auto dot_transform_data_quick =
dot_product(position_data_xy, velocity_data_xy);
// ------------------------
// BinaryProjectArrayTransform and BinaryProjectNormalArrayTransform
// ------------------------
//
// These transforms take in two objects, and either project the first
// input onto the second, or onto the plane normal to it. Note that if the
// direction vector (the second input) isn't a unit vector the projection
// will be scaled by the square of its magnitude
//
// The transform expects the same dimensionality of input data as that of
// the direction vector
auto binary_project_transform = BinaryProjectArrayTransform<2>();
auto binary_normal_project_transform = BinaryProjectNormalArrayTransform<2>();
// Either of the above can be wrapped into a BinaryArrayTransformData object
// acting on two inputs
//
// For example, the result of the following is equivalent to the result of
// dot_product(position_data_xy,velocity_data_xy) * velocity_data_xy
auto binary_project_data = BinaryArrayTransformData(
binary_project_transform, position_data_xy, velocity_data_xy);
return;
}
EXPERIMENTAL Lambda wrappers#
The above unary and binary transform data all rely on built-in transforms or standard operators. There are situations in which those do not allow enough flexibility, so VANTAGE-Reactions offers a wrapper for lambda functions that can be used in binary and unary array transform data objects.
Warning
This is an experimental feature designed to provide device-copyable lambdas. While it works on some common backends, due to the nature of the workaround, there is at least one known issue with the generic adaptivecpp backend! Future work is planned on adressing this, but the wrappers shouldn’t currently be used in production!
void lambda_wrapper_array_transform_examples() {
// 1D and 2D reaction data objects for the examples
auto position_data_x = extract<1>("POSITION");
auto position_data_xy = extract<2>("POSITION");
auto velocity_data_xy = extract<2>("VELOCITY");
// Supported lambdas are binary and unary functions of either
// conforming REAL arrays or of REAL values, allowing for full array
// or elementwise application
// ------------------
// Binary full array
// ------------------
auto binary_lambda_full = [](const std::array<REAL, 2> &a,
const std::array<REAL, 2> &b) {
return std::array<REAL, 2>{a[0] * b[1], b[1]};
};
// The lambda wrapper is templated on the type of the lambda, as well as the
// expected array dimension
auto lambda_wrapper_binary_full =
utils::LambdaWrapper<decltype(binary_lambda_full), 2>{binary_lambda_full};
// batData is a helper function for turning lambda wrappers into
// full array binary transform data object
auto full_binary_transform_data =
batData(lambda_wrapper_binary_full, position_data_xy, velocity_data_xy);
// ------------------
// Binary elementwise
// ------------------
auto binary_lambda_elementwise = [](const REAL &a, const REAL &b) {
return 2 * a + b;
};
auto lambda_wrapper_elementwise =
utils::LambdaWrapper<decltype(binary_lambda_elementwise), 1>(
binary_lambda_elementwise);
// betData is a helper function for turning lambda wrappers into elementwise
// binary transform data objects
auto elementwise_binary_transform_data =
betData(lambda_wrapper_elementwise, position_data_xy, position_data_xy);
// ------------------
// Unary full array
// ------------------
auto unary_lambda_full = [=](const std::array<REAL, 2> &a) {
return std::array<REAL, 2>{a[0] * a[0], a[1] * a[1]};
};
auto unary_lambda_wrapper_full =
utils::LambdaWrapper<decltype(unary_lambda_full), 2>(unary_lambda_full);
// uatData is a helper function templated on the dimensionality of the data
// and the type of the lambda generating unary full array transform data from
// a unary lambda
auto full_unary_transform_data =
uatData<2, decltype(unary_lambda_wrapper_full)>(
unary_lambda_wrapper_full);
// As with regular unary array transforms, these must be in a pipeline to be
// used as standard ReactionData
auto unary_lambda_pipeline =
pipe(position_data_xy, full_unary_transform_data);
// ------------------
// Unary elementwise
// ------------------
auto unary_lambda_elementwise = [](const REAL &a) { return 2 * a; };
auto unary_lambda_wrapper_elementwise =
utils::LambdaWrapper(unary_lambda_elementwise);
// uetData is a helper function generating elementwise unary transform data
// from a lambda
auto elementwise_unary_transform_data =
uetData<2, decltype(unary_lambda_wrapper_elementwise)>(
unary_lambda_wrapper_elementwise);
// Same as above - we need a pipeline in order to make use of unary array
// transform data
auto unary_lambda_pipeline_elementwise =
pipe(position_data_xy, elementwise_unary_transform_data);
return;
}
Surface reaction data#
The data objects in this section are specialised for surface reaction uses, meaning that they require surface interaction data to be present on the particles.
Specular reflection data#
Dimensionality: variable - depends on the velocity dimension (2 or 3). Requires input of the same dimension, representing the ingoing velocities.
Required properties: Simple props: boundary intersection normal; Species props: none
Details: Given a surface normal and an input velocity, reflects the velocity specularly based ont the normal. Should be used as part of a pipeline, allowing for modification of input and output velocities.
Example - see pipeline example
Spherical basis reflection data#
Dimensionality: 3. Requires input of the same dimension, representing the output velocity in (\(v\), \(\theta\), \(\varphi\)), where the first entry is the velocity magnitude, the second the angle with respect to the surface normal, and the third the angle with respect to the initial (pre-reflection) velocity projection onto the surface
Required properties: Simple props: velocity, boundary intersection normal; Species props: none
Details: Given the reflected velocity in spherical coordinates, uses the particle velocity and the surface normal to construct a local basis for reflection. Useful when reflected data is given in spherical coordinates (such as from the TRIM database)
Example:
void spherical_basis_reflection_example() {
// In case we wish to remap the default boundary normal or velocity name
auto used_map = get_default_map();
// v = 2, theta = pi/4, phi = 3*pi/4
std::array<REAL, 3> coords{2.0, M_PI / 4, 3 * M_PI / 4};
// Here we just use a fixed reflection value, but this can
// be calculated using any other reaction data object
auto coord_data = FixedArrayData<3>(coords);
auto spherical_reflection = SphericalBasisReflectionData();
// The pipeline object will first evaluate the reflected coordinate (fixed
// here) and will then pipe this into the reflection kernel.
//
// The reflection kernel will use the particle velocity and the surface normal
// to generate the correct post-reflection basis, where the v, theta, and phi
// coordinate values would be used
//
// NOTE: This will only work with 3-dimensional data
auto pipeline = pipe(coord_data, spherical_reflection);
return;
}
Cartesian basis reflection data#
Dimensionality: 3. Requires input of the same dimension, representing the output velocity in the Cartesian basis defined by the ingoing velocity and the surface normal. The first two components are parallel to the surface (with the first component in the direction determined by the projection of the ingoing particle velocity). The final component is in the direction of the surface normal (directed back into the domain).
Required properties: Simple props: velocity, boundary intersection normal; Species props: none
Details: Given the velocity and surface normal, determines the local basis and sets the outgoing particle velocity based on the input Cartesian components. Useful when the reflected data is given in Cartesian coordinates (such as for thermal reflection)
Example:
void cartesian_basis_reflection_example() {
// In case we wish to remap the default boundary normal or velocity name
auto used_map = get_default_map();
// Reflect in the direction normal to the surface, and back into the domain
std::array<REAL, 3> coords{0, 0, 1};
// Here we just use a fixed reflection value, but this can
// be calculated using any other reaction data object
auto coord_data = FixedArrayData<3>(coords);
auto spherical_reflection = CartesianBasisReflectionData();
// The pipeline object will first evaluate the reflected coordinate (fixed
// here) and will then pipe this into the reflection data object.
//
// The reflection data object will use the particle velocity and the surface
// normal to generate the correct post-reflection basis, where the
// post-reflection velocity is given with local cartesian coordinates
//
// NOTE: This will only work with 3-dimensional data
auto pipeline = pipe(coord_data, spherical_reflection);
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:
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, velocity, 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:
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;
}
General absorption kernels#
Overview: These kernels represent a general absorption process, i.e. anything that removes the weight of particles without creating new particles. Unlike ionisation, it only stores the particle, momentum, and energy sources due to the absorbed particle, and not due to any of the particles that might be interacting with it.
Required properties:
Parent: Simple props: weight, velocity, source_density, source_energy, source_momentum; Species props: N/A
Descendant: N/A; Species props: N/A
Scattering kernel: N/A
Weight kernel: N/A
Transformation kernel: N/A
Feedback kernel: The weight is removed from the parent, and together with the velocity of the particle it determines all three sources. NOTE: Unlike specific kernels, the sources are not species-specific, which means that property remapping is required. This is particularly important when using these kernels to specify surface processes.
Example:
void general_absorption_kernels_example() {
// For this kernel, we do need to remap the sources in general
// For example we might want to use these kernels to represent absorption at a
// surface
auto properties_map = PropertiesMap();
properties_map[VANTAGE::Reactions::default_properties.source_density] =
"SURFACE_SOURCE_DENSITY";
properties_map[VANTAGE::Reactions::default_properties.source_momentum] =
"SURFACE_SOURCE_DENSITY";
properties_map[VANTAGE::Reactions::default_properties.source_energy] =
"SURFACE_SOURCE_DENSITY";
auto target_species = Species("ION", 2.0, 0.0,
-1); // This is absorbed species
auto absorption_kernels =
GeneralAbsorptionKernels<2 // velocity dat dimensionality
>(
target_species, // Absorbed species
properties_map.get_map() // Our remapped sources
);
return;
}
General linear scattering kernels#
Overview: These kernels produce one product, and the velocities of the product are set by the DataCalculator values. Sources (momentum and energy) are then calculated based on the parent and product velocities and the reacted weight.
Required properties:
Parent: Simple props: weight, velocity, source_energy, source_momentum; Species props: N/A
Descendant: internal_state, velocity, weight; Species props: N/A
Scattering kernel: The velocities of the product are set from the values calculated by the DataCalculator of the containing reaction
Weight kernel: All reacted weight is passed onto the product
Transformation kernel: The product internal_state is set from the outgoing particle ID
Feedback kernel: The weight is removed from the parent, and the momentum and energy sources are calculated using the parent and product velocities. NOTE: Unlike specific kernels, the sources are not species-specific, which means that property remapping is required. This is particularly important when using these kernels to specify surface processes. The calculation and writing of sources can be turned off using a template argument (see below).
Example:
void general_linear_scattering_kernels_example(
NP::ParticleGroupSharedPtr particle_group) {
// For this kernel, we do need to remap the sources in general
// For example we might want to use these kernels to represent absorption at a
// surface
//
// In this example we will indeed set up a specular reflection kernel using
// the general linear scattering kernel and specular reflection data
auto properties_map = PropertiesMap();
properties_map[VANTAGE::Reactions::default_properties.source_momentum] =
"SURFACE_SOURCE_DENSITY";
properties_map[VANTAGE::Reactions::default_properties.source_energy] =
"SURFACE_SOURCE_DENSITY";
auto projectile_species = Species("ION", 2.0, 0.0,
-1); // This is the projectile species
// Implementing specular reflection on ingoing particle velocities
auto velocity_data = ExtractorData<2>(NP::Sym<REAL>("VELOCITY"));
auto specular_reflection = SpecularReflectionData<2>();
auto pipeline = pipe(velocity_data, specular_reflection);
// Wrapping the pipeline in a DataCalculator
auto data_calculator = DataCalculator<decltype(pipeline)>(pipeline);
auto scattering_kernel =
LinearScatteringKernels<2, // velocity dat dimensionality
true // Whether sources should be written to
// (remapping not required if false)
>(projectile_species, // Scattered species
properties_map.get_map() // Our remapped sources
);
// The above can then be used to creat a specular reflection reaction
// that can then be applied on particles that have hit a surface
//
// Here we set a constant rate and reflect the particle in the same internal
// state
auto specular_reflection_reaction =
LinearReactionBase<1, FixedRateData, decltype(scattering_kernel),
decltype(data_calculator)>(
particle_group->sycl_target, 0, std::array<int, 1>{0},
FixedRateData(1.0), scattering_kernel, data_calculator);
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.
void electron_impact_ion_example(NP::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:
void recombination_reaction_example(NP::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);
// Recombination with 2 velocity dimensions
auto recombination_reaction =
Recombination<decltype(reaction_rate), decltype(data_calculator), 2>(
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.