Reactions controllers#
Introduction#
While one or two reactions can be applied manually, once multiple species and multiple processes are treated, complexity is sure to explode. To address this, VANTAGE-Reactions offers reaction controller objects, that bundle reactions and transformations on particles, while respecting certain constraints.
For example, let us assume we wish to apply two reactions with an Euler timestep of \(dt=0.1\). Let us assume that the rate of the first reaction is \(K_1=1.0\) and the rate of the second reaction is \(K_2=1.0\). Now let us assume we are applying the reactions to a particle with weight \(w=0.01\). It is easy to see that, if reactions are applied naively and sequentially, all of the weight of the particle would be consumed by whichever reaction is applied first, since both \(K_1 dt\) and \(K_2 dt\) are greater than \(w\). Instead, what we would prefer to happen isthat half of the particle weight is spent on the first and half on the second reaction. This accounting is automatically performed by the reactions, and is enabled by the controllers knowing the order of operations.
Another example of potential complications is the filtering of particles based on their species/internal state, since reactions have a predefined ingoing state. This is also done automatically for each ID detected by the reaction controller amongst the reactions it is responsible for.
Finally, we might wish to perform operations on the products before they are added to the particle group, or on the reactants/parents after reactions are applied. Examples include small particle merging or removal, as well as the projection of particle properties onto the grid.
Reaction controllers provide a centralised interface for managing the above workload.
Reaction controller modes#
Different behaviour can be obtained by passing different mode flags to the ReactionController object.
Currently implemented flags are:
standard_mode - leads to deterministic application, where all reactions are applied to all particles in accordance to their rates
semi_dsmc_mode - semi-deterministic Direct Simulation Monte Carlo (DSMC) mode application, where particles going through reactions are determined through random sampling, and all reactions are applied to them, completely consuming them
Below are the explanations of how each mode is applied.
Deterministic reaction application#
When the standard_mode flag (or no flag) is passed to the controller, it applies all reactions to all eligible particles, splitting the reacting weight amongst the reactions proportional to their rate. The following is the sequence of actions taken by the controller, given an Euler timestep \(dt\):
Filter the passed group based on particle species/internal_state
Run the rate loop on all reactions, passing in the subgroups containing their respective ingoing states
Run the product loop with the requested step \(dt\), storing all the produced particles into a separate child group
Perform any transformations on the child group, making sure they are species-wise
Perform any transformations on the post-reaction parents
Add the child group to the parent group, completing the application of reactions
The above is done cell-block-wise, so that buffers do not overflow when many reactions are applied. The user can control the maximum number of particles per cell (on average) and the number of cells per block. The defaults are intentionally greedy, and should be reduced in case of memory issues.
ReactionController and using it with multiple reactions and transforms#inline void reaction_controller_example(ParticleGroupSharedPtr particle_group) {
auto particle_spec = particle_group->get_particle_spec();
auto prop_map = get_default_map();
// We shall generate a CX and two ionisation reactions quickly, in order to
// populate the controller with some reactions
auto electron_species = Species("ELECTRON", 5.5e-4, -1.0);
auto ion_species_1 = Species("ION", 1.2, 0.0, 0);
auto ion_species_2 = Species("ION2", 2.0, 0.0, 1);
// We'll use dummy rate data to simplify the setup
auto rate_data = FixedRateData(1.0);
// CX beam data
auto vx_beam_data = FixedRateData(1.0);
auto vy_beam_data = FixedRateData(-1.0);
auto data_calculator =
DataCalculator<FixedRateData, FixedRateData>(vx_beam_data, vy_beam_data);
auto cx_kernel = CXReactionKernels<2>(ion_species_1, ion_species_2, prop_map);
// CX reaction
auto cx_reaction = std::make_shared<
LinearReactionBase<1, FixedRateData, CXReactionKernels<2>,
DataCalculator<FixedRateData, FixedRateData>>>(
particle_group->sycl_target, ion_species_1.get_id(),
std::array<int, 1>{static_cast<int>(ion_species_2.get_id())}, rate_data,
cx_kernel, data_calculator);
// Ionisation reactions
auto ionise_reaction_1 =
std::make_shared<ElectronImpactIonisation<FixedRateData, FixedRateData>>(
particle_group->sycl_target, rate_data, rate_data, ion_species_1,
electron_species);
auto ionise_reaction_2 =
std::make_shared<ElectronImpactIonisation<FixedRateData, FixedRateData>>(
particle_group->sycl_target, rate_data, rate_data, ion_species_2,
electron_species);
// We can now initialise a reaction controller and populate it with the above
// reactions We start off by creating transformations for the children and
// parents
//
// We will remove any children/parents below some set weight using the
// following transform
auto remove_wrapper = std::make_shared<TransformationWrapper>(
std::vector<std::shared_ptr<MarkingStrategy>>{
make_marking_strategy<ComparisonMarkerSingle<REAL, LessThanComp>>(
Sym<REAL>(prop_map[default_properties.weight]), 1e-6)},
make_transformation_strategy<SimpleRemovalTransformationStrategy>());
// But will first try merge any children/parents below a higher weight
// threshold
auto merge_transform =
make_transformation_strategy<MergeTransformationStrategy<2>>();
auto merge_wrapper = std::make_shared<TransformationWrapper>(
std::vector<std::shared_ptr<MarkingStrategy>>{
make_marking_strategy<ComparisonMarkerSingle<REAL, LessThanComp>>(
Sym<REAL>(prop_map[default_properties.weight]), 1e-2)},
merge_transform);
auto reaction_controller = ReactionController(
std::vector{
merge_wrapper,
remove_wrapper}, // the order matters! this will first merge parents,
// then remove any remaining small particles
std::vector{merge_wrapper, remove_wrapper}
// this will do the same to the children
// before merging them into the parents
);
reaction_controller.set_cell_block_size(
256); // This is the greedy default value, reduce this if memory issues
// are found
reaction_controller.set_max_particles_per_cell(
16384); // This is the default maximum (average) number of particles per
// cell for the use in reaction data buffers, modify as needed
// Now we can add the reaction simply
reaction_controller.add_reaction(cx_reaction);
reaction_controller.add_reaction(ionise_reaction_1);
reaction_controller.add_reaction(ionise_reaction_2);
// We can now request an Euler step of 0.01 time units in standard deterministic mode
reaction_controller.apply(particle_group, 0.01, ControllerMode::standard_mode);
return;
}
Semi-DSMC reaction application#
By passing ControllerMode:semi_dsmc_mode to the controller the following behaviour is obtained:
Filter the passed group based on particle species/internal_state
Run the rate loop on all reactions, passing in the subgroups containing their respective ingoing states
Sample uniform random numbers and compare them to \(1-exp(-K_{tot}dt)\), where \(K_{tot}\) is the total reaction rate for any given particle. If the sampled number is lower than that value, the particle is flagged as going through reactions
Run the product loop with the requested step \(dt\), but only applied to the reacted particles, and using all of their weights instead of allowing for some unreacted fraction. This is where the method differes from the standard_mode case.
Perform any transformations on the child group, making sure they are species-wise
Perform any transformations on the post-reaction parents
Add the child group to the parent group, completing the application of reactions
In order to use this mode, the set_rng_kernel() method should be called on the controller before applying it, and the RNG kernel should produce uniform numbers between 0 and 1. For an example of such a kernel see filtered Maxwellian sampling examples.