Marking and transformations#
Introduction#
A common use case when dealing with NESO-Particles ParticleGroup objects is selecting particles based on some condition(s) and applying transformations to them.
For example, assuming a ParticleGroup contains many different species of particles, one might want to select all particles of a given species that also have weight below some threshold, and then merge them to avoid tracking many small particles.
In order to accommodate the above, VANTAGE-Reactions offers a uniform interface for MarkingStrategy, TransformationStrategy, and TransformationWrapper objects.
Marking Strategies#
A marking strategy is the abstract wrapper class for the creation of NESO-Particles particle subgroups. The main marking strategy is the MarkingStrategyDirect, which is effectively a closure for the NESO-Particles subgroup constructor (see below how this enables transformation wrappers).
void direct_marking_example(NP::ParticleGroupSharedPtr particle_group) {
// Here we create a marking strategy marking low weight particles
auto marking_strategy = make_direct_marking_strategy(
"test_strategy", // Name of the strategy for profiling purposes
[](auto w) { return w[0] < 1e-6; }, // Marking kernel
NP::Access::read(NP::Sym<REAL>("WEIGHT")) // Accessors for the kernel
);
// The subgroup can then be created as follows from another subgroup
auto subgroup_low_weight = marking_strategy->make_marker_subgroup(
particle_sub_group(particle_group));
// Contrast the above with the full constructor call from NESO-Particles
auto subgroup_low_weight_from_NP = particle_sub_group(
particle_group, [](auto w) { return w[0] < 1e-6; },
NP::Access::read(NP::Sym<REAL>("WEIGHT")));
return;
}
Transformation Strategies#
Once a suitable subgroup is constructed, we use TransformationStrategy objects to apply any required transformation to those particles. All transformation strategies have the transform method which is applied to particle subgroups. Examples of the various built-in strategies are given below.
Particle Removal Strategy#
Often we wish to remove some particles from the simulation. For example, this might be due to their weight being too low to track. The following code will remove all particles with weight below a given threshold using marking and transformation strategies.
void removal_strategy_example(NP::ParticleGroupSharedPtr particle_group) {
auto subgroup_low_weight = particle_sub_group(
particle_group, [](auto w) { return w[0] < 1e-6; },
NP::Access::read(NP::Sym<REAL>("WEIGHT")));
// The make_transformation_strategy helper function casts concrete
// transformation strategies into std::shared_ptr<TransformationStrategy>.
//
// The simple removal strategy below requires no inputs, and just deletes
// all particles in the passed subgroup
auto removal_strategy =
make_transformation_strategy<SimpleRemovalTransformationStrategy>();
removal_strategy->transform(subgroup_low_weight);
return;
}
ParticleDatZeroer#
VANTAGE-Reactions assumes that particles carry information of their contribution to various sources that need to be projected onto the grid. A common requirement in these situations is to reset the values of sources on the particles after projection. The library offers the ParticleDatZeroer transformation strategy that allows the user to accomplish this.
void zeroer_strategy_example(NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = std::make_shared<NP::ParticleSubGroup>(particle_group);
// A ParticleDatZeroer zeroes INT or REAL particle dats and is
// constructed by passing a vestor of strings with the dat names
auto zeroer = make_transformation_strategy<ParticleDatZeroer<REAL>>(
std::vector<std::string>{"ELECTRON_SOURCE_DENSITY",
"ION_SOURCE_DENSITY"});
zeroer->transform(input_subgroup);
return;
}
Accumulator Strategies#
Another common requirement is the accumulation of particle properties cellwise. This is a requirement for finite volume methods (projection of sources) as well as general particle data analysis (weighted averages of quantities). Three classes of transformation strategies are provided for this use:
CellwiseAccumulator - accumulating one or more properties cellwise
WeightedCellwiseAccumulator - accumulating one or more properties cellwise while weighing them with the particle weights
CellwiseReactionDataAccumulator - accumulating the result of a reaction data object, for use in cases when the first two are two restrictive
void accumulator_strategy_example(NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = particle_sub_group(particle_group);
// CellwiseAccumulators need to access some general data about the particle
// group, so it needs to be available on construction
//
// Similarly to zeroers, accumulators take in the names of the particle dats
// that they should accumulate values for cellwise. Here we use make_shared
// instead of make_transformation_strategy in order to be able to call
// accumulator-specific methods
auto accumulator = std::make_shared<CellwiseAccumulator<REAL>>(
particle_group, std::vector<std::string>{"ELECTRON_SOURCE_DENSITY",
"ION_SOURCE_DENSITY"});
// The accumulator, unlike other transforms, does not modify the particle
// group. Instead, it modifies its own internal state.
accumulator->transform(input_subgroup);
// Upon accumulation, the accumulated data is stored in NESO-Particle
// NP::CellDatConst objects and can be retrieved easily
auto accumulated_electron_source =
accumulator->get_cell_data("ELECTRON_SOURCE_DENSITY");
// Individual cell data buffers can be zeroed
accumulator->zero_buffer("ELECTRON_SOURCE_DENSITY");
// Or all buffers can be zeroed
accumulator->zero_all_buffers();
// A weighted accumulator transform is also available, taking in the REAL
// particle dat to be used as the weight (should be a PartilceDat with 1
// component) in the constructor
auto weighted_accumulator =
std::make_shared<WeightedCellwiseAccumulator<REAL>>(
particle_group, std::vector<std::string>{"VELOCITY", "POSITION"},
"WEIGHT");
weighted_accumulator->transform(input_subgroup);
auto weighted_velocities = weighted_accumulator->get_cell_data("VELOCITY");
weighted_accumulator->zero_buffer("VELOCITY");
// In addition to the standard accumulator methods, the weighted accumulator
// also offers access to the accumulated weight dat
auto accumulated_weight = weighted_accumulator->get_weight_cell_data();
return;
}
void reaction_data_accumulator_strategy_example(
NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = particle_sub_group(particle_group);
// For the reaction data accumulator, we need to define the reaction
// data object whose results should be accumulated cellwise
//
// Here we use extractors and reaction data arithmetic to accumulate
// weighted kinetic energy in each velocity dimension
auto weight = extract<1>("WEIGHT");
auto velocity = extract<2>("VELOCITY");
auto kin_energy = weight * velocity * velocity;
// The accumulator is then defined simply as
auto accumulator =
std::make_shared<CellwiseReactionDataAccumulator<decltype(kin_energy)>>(
particle_group, kin_energy);
// The accumulator can then be called like any other transform
accumulator->transform(input_subgroup);
// Upon accumulation, the accumulated data is stored in NESO-Particle
// NP::CellDatConst objects and can be retrieved easily
auto accumulated_kin_energy = accumulator->get_cell_data();
// The buffer can be zeroed as
accumulator->zero_buffer();
return;
}
Cellwise distributor strategy#
Similar to the requirement for accumulating properties cellwise, there are situations where we want to broadcast one or more property onto all particles in a cell. The CellwiseDistributor transformation strategy offers this, working like the inverse of the CellwiseAccumulator.
void distributor_strategy_example(NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = particle_sub_group(particle_group);
// CellwiseDistributors need to access some general data about the particle
// group, so it needs to be available on construction
//
// Similarly to zeroers, distributors take in the names of the particle dats
// that they should distribute values for cellwise. Here we use make_shared
// instead of make_transformation_strategy in order to be able to call
// distributor-specific methods
auto distributor = std::make_shared<CellwiseAccumulator<REAL>>(
particle_group, std::vector<std::string>{"ELECTRON_SOURCE_DENSITY",
"ION_SOURCE_DENSITY"});
// To set the values for distribution we can do the following
auto buffer = distributor->get_cell_data("ELECTRON_SOURCE_DENSITY");
auto num_cells = particle_group->domain->mesh->get_cell_count();
for (int cellx = 0; cellx < num_cells; cellx++) {
buffer[cellx]->at(0, 0) = 1.0; // at 0th component and 0th layer (see
// NESO-Particles documentation)
}
distributor->set_cell_data("ELECTRON_SOURCE_DENSITY", buffer);
distributor->transform(input_subgroup);
// Distributor buffer zeroing, if needed, can be done in the same way as for
// the accumulator
distributor->zero_buffer("ELECTRON_SOURCE_DENSITY");
distributor->zero_all_buffers();
return;
}
Composite Strategy#
Sometimes multuple strategies need to be applied in order. It is possible to compose transformation strategies by adding them to a composite strategy, allowing all of them to be applied with one transform call. This is particularly useful in the construction of TransformationWrapper objects (see below), where there is a hook left for a single transformation strategy.
void composite_strategy_example(NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = std::make_shared<NP::ParticleSubGroup>(particle_group);
// We wish to compose an accumulator and a particle dat zeroer to accumulate
// some sources and then to reset the particle data that stored them
//
// The sources we wish to accumulate
auto source_names =
std::vector<std::string>{"ELECTRON_SOURCE_DENSITY", "ION_SOURCE_DENSITY"};
// In order to have later access to the accumulator, we construct it using
// make_shared
auto accumulator =
std::make_shared<CellwiseAccumulator<REAL>>(particle_group, source_names);
// A composite transform can be constructed by passing a vector of
// TransformationStrategy objects, so if we wish to include the accumulator,
// it must be dynamically cast
auto composite = std::make_shared<CompositeTransform>(
std::vector<std::shared_ptr<TransformationStrategy>>{
std::dynamic_pointer_cast<TransformationStrategy>(accumulator)});
// We can also directly add TransformationStrategy objects to the composite to
// be applied in sequence (in order of addition)
composite->add_transformation(
make_transformation_strategy<ParticleDatZeroer<REAL>>(source_names));
// The composite can then be applied as one transformation
composite->transform(input_subgroup);
// Since the accumulator was added to the composite via a pointer
// its interface is still accessible
auto accumulated_electron_source =
accumulator->get_cell_data("ELECTRON_SOURCE_DENSITY");
accumulator->zero_all_buffers();
return;
}
Downsampling Strategies#
A common problem in weighted particle methods is ensemble management, and in particular downsampling. VANTAGE-Reactions offers a framework for building downsampling transformation strategies, with a small number of them supplied through helper functions.
In general, downsampling strategies consist of the following steps:
Reduction - where a number of quantities across the particle ensemble are reduced, i.e. moments or other quantities are calculated
Downsampling - where the properties of the particles are modified in such a way that some of them are effectively marked for removal, while the remaining particles are modified according to the downsampling algorithm.
Removal - particles effectively “marked” for removal are removed
The above can be applied on multiple downsampling groups separately, and the downsampling strategies that assume grouping expect that it has been prepared beforehand, by default using the grouping_index integer property. See uniform velocity binning below as an example of a binning transform.
Vranic Merging Strategy#
VANTAGE-Reactions implements a version of the merging algorithim from [VRANIC2015]. It assumes that all particles being merged are of the same species (i.e. have the same mass) and that they are non-relativistic. In the original paper, the authors merge particles within momentum space cells, while we merge all particles in the downsampling group, which could be a momentum/velocity space cell, but doesn’t have to be. Correspondingly, in 3D we use the momentum space bounding box in 3D to determine the plane in which the merged particle momenta lie.
Particles are merged cell-wise and downsampling-group-wise into 2 particles. The properties modified by the merging algorithm are the weights and momenta/velocities. Other properties are taken from 2 other particles in the passed subgroups, i.e. properties like cell and species IDs should be copied consistently, but all other properties should be considered undefined. As such, merging should only be invoked once all particle properties have been used for their respective purposes, such as recording sources. Note that the above means that the first two particles’ positions in the downsampling group will be used as the merged particle positions.
void vranic_merging_strategy_example(
NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = std::make_shared<NP::ParticleSubGroup>(particle_group);
auto prop_map = get_default_map();
auto merging_strat =
make_vranic_merging_strategy<2 // The velocity space dimensionality
>(particle_group, // The particle group
1, // The number of downsampling groups
// (e.g. velocity bins, etc.)
prop_map // Property map used for
// remapping the grouping index, weight,
// linear index, and velocity properties
);
merging_strat->transform(input_subgroup);
return;
}
Simple Thinning Strategy#
A classic alternative to merging is particle thinning, i.e. removing some particles randomly while modifying the properties of the rest. The simple thinning strategy keeps particles with some probability - the thinning_ratio, scaling their weights with the inverse of that probability, while removing the rest of the particles. This procedure conserves particle weight on average only.
void simple_thinning_strategy_example(
NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = std::make_shared<NP::ParticleSubGroup>(particle_group);
auto prop_map = get_default_map();
// Below is a placeholder rng kernel, in practice this would be a uniformly
// sampled value - here we just use a constant number
auto rng_lambda = [&]() -> REAL { return 0.01; };
auto rng_kernel = NP::host_per_particle_block_rng<REAL>(rng_lambda, 1);
auto thinning_strat = make_simple_thinning_strategy(
particle_group, // The particle group
0.1, // The thinning ratio - will on average keep 10% of the particles
rng_kernel, // The uniform random variate rng kernel
prop_map // Property map used for
// remapping the weight, and the panic flag
);
thinning_strat->transform(input_subgroup);
return;
}
DEPRECATED Legacy Merging Strategy#
This is a purely cell-wise version of the above Vranic merging strategy, which also merges particles into the centre of mass. Both of these are drawbacks and the transformation is likely to be removed in the future.
The implementation of of the legacy MergeTransformationStrategy is available in 2D and 3D:
void merging_strategy_example(NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = std::make_shared<NP::ParticleSubGroup>(particle_group);
auto merging_strat =
make_transformation_strategy<MergeTransformationStrategy<2>>();
merging_strat->transform(input_subgroup);
return;
}
Direct transformation strategies#
In cases where the user wants to apply a custom lambda function to a particle subgroup or call a particular ParticleLoop as a transformation strategy VANTAGE-Reactions supplies TransformationStrategyDirect and TransformationStrategyLambda.
void direct_transformation_example(NP::ParticleGroupSharedPtr particle_group) {
auto subgroup_low_weight = particle_sub_group(
particle_group, [](auto w) { return w[0] < 1e-6; },
NP::Access::read(NP::Sym<REAL>("WEIGHT")));
// We can recreate the built-in removal strategy using a direct transformation
// strategy
auto removal_strategy_direct = make_lambda_transformation_strategy(
"removal_lambda", // The name for this transformation strategy
[](auto target) {
target->get_particle_group()->remove_particles(target);
} // The lambda to be applied to any passed subgroup
);
removal_strategy_direct->transform(subgroup_low_weight);
// Or we can apply a particle loop as a transformation strategy
auto set_id_strategy = make_direct_transformation_strategy(
"set_id_0", // Name of the strategy
[](auto id) { id.at(0) = 0; }, // The particle loop kernel to be applied
NP::Access::write(NP::Sym<INT>("ID")) // Accessors for the loop
);
// Apply to the particle_group
set_id_strategy->transform(particle_sub_group(particle_group));
return;
}
Uniform velocity space binning strategy#
A simple strategy that bins particles into uniform velocity cells is provided with VANTAGE-Reactions. It splits each of the velocity space dimensions into some number of uniform cells, up to a given extent, and in addition adds guard cells used to bin particles that might be outside of the extents. The resulting linear bin index is recorded in an integer particle property.
For example, the below strategy will bin particles into a core binning region spanning \((-1.5,1.5] \times (-1.5,1.5]\) split into 10 by 10 cells, and with an outer layer of guard cells capturing any particles with velocity components outside of the binning region - resulting in a total of 144 binning cells.
void uniform_velocity_binning_example(
NP::ParticleGroupSharedPtr particle_group) {
auto input_subgroup = std::make_shared<NP::ParticleSubGroup>(particle_group);
auto prop_map = get_default_map();
auto velocity_bin =
uniform_velocity_bin_transform<2 // Velocity space dimensionality
>(
std::array<REAL, 2>{3.0,
3.0}, // Extent in each dimension of
// the main binning cells - corresponding to
// ranges (-1.5,1.5]x(-1.5,1.5]
std::array<INT, 2>{10, 10}, // Number of main binning cells in each
// dimension - will result in 12 x 12
// total cells, accounting for guard cells
NP::Sym<INT>("REACTIONS_GROUPING_INDEX"), // The linear velocity
// indexing NP::Sym
NP::Sym<REAL>("VELOCITY") // The velocity NP::Sym
);
velocity_bin->transform(input_subgroup);
return;
}
Transformation Wrappers#
Often we wish to encapsulate both some marking conditions as well as the transformation we wish to perform into a single object.
A common example is performing merging on multiple different species, but only on particles where weight is less than some threshold. In other words, the transformation we wish to perform is fixed, but only some of the marking conditions are fixed, while others vary. In this case the library offers the TransformationWrapper class, which wraps marking conditions and a transformation strategy, and applies them to a ParticleGroup. We can then fix the instruction “Merge all particles with weight < threshold”, and can extend it with other conditions, such as “and with species ID = 1”.
TransformationWrapper being used to remove particles with low weights for two different ID values#void transformation_wrapper_example(NP::ParticleGroupSharedPtr particle_group) {
// A transformation wrapper can be constructed with a vector of marking
// strategies or they can be added later.
//
// However, it always requires a transformation strategy at construction
//
// The wrapper below encapsulates the instruction "remove all particles with
// weights < 1e6"
auto wrapper = std::make_shared<TransformationWrapper>(
std::vector<std::shared_ptr<MarkingStrategy>>{
make_direct_marking_strategy(
"low_weight", [](auto w) { return w[0] < 1e-6; },
NP::Access::read(NP::Sym<REAL>("WEIGHT")))},
make_transformation_strategy<SimpleRemovalTransformationStrategy>());
// Wrappers can be copied and/or extended
//
// For example, maybe we want to remove both particles with ID=0 and
// ID=1
//
// For this we create 2 wrappers using the above as a base
auto wrapper_ID0 = std::make_shared<TransformationWrapper>(*wrapper);
auto wrapper_ID1 = std::make_shared<TransformationWrapper>(*wrapper);
// We can then add different further marking conditions to them
wrapper_ID0->add_marking_strategy(
VANTAGE::Reactions::make_direct_marking_strategy(
"id0", [](auto id) { return id[0] == 0; },
NP::Access::read(NP::Sym<INT>("ID"))));
wrapper_ID1->add_marking_strategy(
VANTAGE::Reactions::make_direct_marking_strategy(
"id1", [](auto id) { return id[0] == 1; },
NP::Access::read(NP::Sym<INT>("ID"))));
// Wrappers, unlike strategies, can act on both groups or subgroups
//
// Subselection is assumed to be performed by the successive
// application of MarkingStrategies
wrapper_ID0->transform(particle_group);
wrapper_ID1->transform(particle_group);
return;
}
Vranic et al. - Particle merging algorithm for PIC codes https://www.sciencedirect.com/science/article/pii/S0010465515000405