Îto diffusion

The Îto diffusion model advances computational particles as drifting Brownian walkers

(9)\[\Delta\mathbf{X} = \mathbf{V}\Delta t + \sqrt{2D\Delta t}\mathbf{W}\]

where \(\mathbf{X}\) is the spatial position of a particle, \(\mathbf{V}\) the particle drift velocity, and \(D\) is the diffusion coefficient in the continuum limit. The vector term \(\mathbf{W}\) indicates a random number sampled from a Gaussian distribution with mean value of 0 and standard deviation of 1.

Tip

The code for Îto diffusion is given in /Source/ItoDiffusion.

ItoParticle

The ItoParticle is used as the underlying particle type for running the Ito drift-diffusion solvers. It is a Struct-of-Arrays payload (see ParticleSoA) whose columns are

struct ItoParticle
{
  ParticleReal mobility  = 0.0; ///< Mobility coefficient.
  ParticleReal diffusion = 0.0; ///< Diffusion coefficient.
  ParticleReal energy    = 0.0; ///< Average particle energy.
  ParticleReal scratch   = 0.0; ///< Scratch scalar storage.

  double old_x = 0.0; ///< Previous position, x-component.
  double old_y = 0.0; ///< Previous position, y-component.
#if CH_SPACEDIM == 3
  double old_z = 0.0; ///< Previous position, z-component.
#endif

  ParticleReal vx = 0.0; ///< Interpolated velocity, x-component.
  ParticleReal vy = 0.0; ///< Interpolated velocity, y-component.
#if CH_SPACEDIM == 3
  ParticleReal vz = 0.0; ///< Interpolated velocity, z-component.
#endif

  ParticleReal scratch_x = 0.0; ///< Scratch vector storage, x-component.
  ParticleReal scratch_y = 0.0; ///< Scratch vector storage, y-component.
#if CH_SPACEDIM == 3
  ParticleReal scratch_z = 0.0; ///< Scratch vector storage, z-component.
#endif
};

In addition to the container-owned position and weight, ItoParticle stores the payload columns above. These extra fields are used for storing the following information in the particle:

1. Mobility, diffusion coefficient, energy (not currently used), and a holder for a scratch scalar storage.

2. The previous particle position, the velocity, and a holder for a RealVect scratch storage.

Tip

Several member functions are available for obtaining the particle properties. See the full ItoParticle C++ API

ItoSolver

The ItoSolver class encapsulates the implementation of Eq. 9 in chombo-discharge. This class can advance a set of computational particles (see ItoParticle) with the following functionality:

1. Move particles with a microscopic drift-diffusion model.

2. Compute particle intersection with embedded boundaries and domain edges.

3. Deposit particles and other particle types on the mesh.

4. Interpolate velocities and diffusion coefficients to the particle positions.

5. Manage superparticle splitting and merging.

Internally, ItoSolver stores its particles in various ParticleContainer<ItoParticle> containers. Although the particle velocities and diffusion coefficients can be manually assigned, they can also be interpolated from the mesh. ItoSolver stores the following properties on the mesh:

1. Mobility.

2. Diffusion coefficient.

3. Velocity function.

The reason for storing both the mobility and velocity function is simply to improve flexibility when assigning the particle velocity \(\mathbf{V}\). Note that the velocity function does not have to represent the particle velocity. When using both the mobility and velocity function, one can compute the particle velocity as \(\mathbf{V} = \mu\mathbf{v}\), where \(\mathbf{v}\) is a velocity field. This is typically done for discharge simulations where for simplicity we assign \(\mathbf{v}\) to be the electric field, and \(\mu\) to the field-dependent mobility. Additional information is available in Particle interpolation.

ItoSpecies

ItoSpecies is a class for parsing solver information into ItoSolver, e.g., whether or not the particle type is mobile or not. The constructor for the ItoSpecies class is

/**
 * @brief Full constructor
 * @param[in] a_name         Species name
 * @param[in] a_chargeNumber Charge number
 * @param[in] a_mobile       Mobile species or not
 * @param[in] a_diffusive    Diffusive species or not
 */
ItoSpecies(const std::string& a_name, const int a_chargeNumber, const bool a_mobile, const bool a_diffusive);

Here, a_name indicates a variable name for the solver. This variable will be used in, e.g., error messages and I/O functionality. a_chargeNumber indicates the charge number of the species and the two booleans a_mobile and a_diffusive indicate whether or not the solver is mobile or diffusive.

Note

The C++ ItoSpecies API is available at https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classItoSpecies.html.

Supplying initial data

Initial data for the ItoSolver is provided through ItoSpecies by providing it with the following:

1. Initial particles specified from a container (ParticleSoA<ItoParticle>) of particles.

2. Provide a density description from which initial particles are stochastically sampled within each grid cell.

In particular, there are two data members that must be populated:

/**
 * @brief Initial particles
 */
ParticleSoA<ItoParticle> m_initialParticles;

/**
 * @brief Initial density, in case the user wants to generate particles from a density distribution
 */
std::function<Real(const RealVect& x, const Real& t)> m_initialDensity;

These can either be populated during construction, or explicitly supplied via the following set functions:

/**
 * @brief Set the initial species density
 * @param[in] a_initialDensity Initial density.
 */
virtual void
setInitialDensity(const std::function<Real(const RealVect& x, const Real& t)>& a_initialDensity);

/**
 * @brief Get initial particles -- this is called by ItoSolver when filling the solver with initial particles.
 * @return Returns m_initialParticles
 */
ParticleSoA<ItoParticle>&
getInitialParticles();

When ItoSolver initializes the data in the solver, it will copy the particle list m_initialParticles from the species and into the solver.

Tip

When using MPI, the user must ensure that each MPI rank does not provide duplicate particles. The ParticleOps class contains lots of supporting functionality for sampling particles with MPI, see the ParticleOps C++ API

When sampling particles from a mesh-based density, the solver will generate the particles so that the specified density is approximately reached within each grid cell. If the density that is supplied does not lead to an integer number of particles in the grid cell (which is virtually always the case), the evaluation of the number of particles is stochastically evaluated. E.g., if the density is \(\phi\) and the grid cell volume is \(\Delta V\), and \(\phi\Delta V = 1.2\), then there is a 20% chance that there will be generated two particles within the grid cell, and 80% chance that only one particle will be generated.

Tip

The number of initially sampled particles is set through ItoSolver.ppc_restart.

Particle containers

Internally, ItoSolver contains several ParticleContainer<ItoParticle> for storing various categories of particles. These categories exist because the transport kernel will almost always lead to particles that leave the domain or intersect the EB. Chemistry models that use ItoSolver for tracking particles might also require new particles to be added into the domain.

ItoSolver defines an enum WhichContainer for classification of ParticleContainer<ItoParticle> data holders for holding particles that live on:

• Main particles (WhichContainer::Bulk).

• The embedded boundary (WhichContainer::EB).

• On the domain edges/faces (WhichContainer::Domain).

• Representing ‘’source particles’’ (WhichContainer::Source).

• Particles that live inside the EB (WhichContainer::Covered).

The particles are available from the solver through the function

/**
 * @brief Get a general particle container
 * @param[in] a_container Which container to fetch.
 * @return Particles
 */
virtual ParticleContainer<ItoParticle>&
getParticles(WhichContainer a_container);

Usually, ItoSolver will perform a drift-diffusion advance and the user will then check if some of the particles crossed into the EB. The solver can then automatically fill the boundary particles containers, see Particle intersection.

Remapping particles

ItoSolver has two functions for remapping particles:

/**
 * @brief Remap the bulk particle container.
 */
virtual void
remap();

/**
 * @brief Remap all particles in the input container
 * @param[in] a_container Particle container
 */
virtual void
remap(WhichContainer a_container);

The bottom function lets the user remap any ParticleContainer<ItoParticle> that lives in the solver. Here, a_container indicates which particle container to remap.

Particle deposition

ItoSolver contains several member functions for depositing various particle properties onto the mesh. The most general version is given below:

/**
 * @brief Deposit a gathered per-particle quantity on the mesh (kappa-conservative + redistribution).
 * @param[out] a_phi                  Mesh data -- must have exactly one component.
 * @param[in]  a_particles            SoA particles to be deposited.
 * @param[in]  a_deposition           Deposition method.
 * @param[in]  a_coarseFineDeposition Coarse-fine deposition strategy.
 * @param[in]  a_gather               Per-particle value gatherer (leaf, index) -> Real.
 * @tparam     Gather                 Callable (const ParticleSoA<ItoParticle>&, std::size_t) -> Real.
 * @note This leaves coarse levels un-averaged and ghost cells stale. Call coarsenAndFillGhosts() afterwards unless
 * the result is a term in a sum that the caller synchronizes itself.
 */
template <typename Gather>
void
depositGathered(EBAMRCellData&                        a_phi,
                const ParticleContainer<ItoParticle>& a_particles,
                DepositionType                        a_deposition,
                CoarseFineDeposition                  a_coarseFineDeposition,
                Gather                                a_gather) const;

This version permits the user to deposit an arbitrary per-particle quantity from a particle container a_particles onto some pre-allocated mesh storage a_phi. The quantity to be deposited is supplied through the a_gather callable, which returns a Real value for each particle in the SoA container.

Important

The ItoSolver deposition methods are specified in the input script, see Input options. Both the base deposition scheme (e.g., NGP or CIC) must be specified, as well as the handling near refinement boundaries.

A simpler version that deposits the bulk particles as a density on the mesh is

/**
 * @brief Deposit particles on to mesh.
 * @param[in] a_container Which container to deposit.
 * @details This will deposit mass (i.e., computational weight) of the the input particle container particles onto the
 * classes member 'm_phi'.
 * @note Calls the general version with arguments: m_phi, m_particles.at(a_container), m_deposition.
 */
virtual void
depositParticles(WhichContainer a_container);

The particles are deposited into the class member m_phi, which stores the particle density on the mesh. This data can then be fetched with

/**
 * @brief Get the mesh data.
 * @return Returns m_phi
 */
virtual EBAMRCellData&
getPhi();

For the full list of available deposition functions, see the ItoSolver C++ API https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classItoSolver.html.

AMR synchronization after a deposit

The deposition functions put the particles on the mesh and, if redistribution is enabled, redistribute the cut-cell mass. They do not average the result down onto the coarser levels, and they do not fill its ghost cells. Whether that is wanted depends on what happens next, which only the caller knows: a deposit that is the final state of a field needs it, whereas a deposit that is merely one term of a larger sum does not, because only the assembled sum needs synchronizing. Callers that need it should therefore call

/**
 * @brief Coarsen the input data and interpolate its ghost cells.
 * @details This is the AMR synchronization that a deposit does *not* do for itself. Whether it is wanted depends on
 * what the caller does next: a deposit that is the final state of a field needs it, whereas a deposit that is one
 * term in a sum does not -- there, only the assembled sum needs synchronizing, and doing it per term produces values
 * that are immediately overwritten. Deposits therefore leave it to the caller.
 * @param[in,out] a_phi Cell-centered mesh data (one component).
 */
void
coarsenAndFillGhosts(EBAMRCellData& a_phi) const;

after depositing. Note that ItoSolver::depositParticles already does this internally, so m_phi (i.e. the data returned by getPhi()) is always synchronized.

Deposition of other quantities

One can also deposit the following quantities on the mesh:

• Conductivity, which deposits \(\mu W\).

• Diffusivity, which deposits \(D W\).

Here, \(W\) is the particle weight, \(\mu\) is the particle mobility, \(D\) is the particle diffusion coefficient. It is up to the user to first interpolate or directly set the particle mobilities and diffusion coefficients before depositing the conductivity onto the mesh.

Functionality for the above deposited quantities exist as the following functions:

/**
 * @brief Deposit conductivities (i.e. mass*mobility / volume)
 * @details This deposits mass*mobility (not multiplied by charge)
 * @param[out] a_phi                  Mesh data
 * @param[in]  a_particles            Particle data
 * @param[in]  a_deposition           Deposition method
 * @param[in]  a_coarseFineDeposition Coarse-fine deposition method.
 */
virtual void
depositConductivity(EBAMRCellData&                  a_phi,
                    ParticleContainer<ItoParticle>& a_particles,
                    DepositionType                  a_deposition,
                    CoarseFineDeposition            a_coarseFineDeposition) const;
/**
 * @brief Deposit diffusivity (i.e. mass*D/volume)
 * @details This deposits mass*mobility (not multiplied by charge)
 * @param[out] a_phi                  Mesh data
 * @param[in]  a_particles            Particle data
 * @param[in]  a_deposition           Deposition method
 * @param[in]  a_coarseFineDeposition Coarse-fine deposition method.
 */
virtual void
depositDiffusivity(EBAMRCellData&                  a_phi,
                   ParticleContainer<ItoParticle>& a_particles,
                   DepositionType                  a_deposition,
                   CoarseFineDeposition            a_coarseFineDeposition) const;

Particle interpolation

Interpolating particle velocities for ItoSolver is done by interpolating the mobility and particle velocities to the mesh,

\[\mathbf{V} = \mu\left(\mathbf{X}\right) \mathbf{v}\left(\mathbf{X}\right).\]

There is, however, some freedom in choosing how the mobility coefficient is calculated, which is discussed below. In either case, there is some interpolation from a mesh-based variable onto the particle position \(\mathbf{X}\). This interpolation method is always parsed from an options file, and is usually an NGP or CIC scheme.

Important

When interpolating particle properties from the mesh, the user must first ensure that ghost cells are properly updated.

The separation into a mobility function and a velocity field is motivated by the introduction of an electric conductivity that permits a rather simple velocity relation as \(\mathbf{v} = \mu\mathbf{E}\), where \(\mathbf{E}\) is the electric field. Complete interpolation of the particle velocity consists of calling two functions:

/**
 * @brief Interpolate mobilities
 * @details This will switch between the two ways of computing the particle mobility.
 */
virtual void
interpolateMobilities();
/**
 * @brief Interpolate the particle velocities.
 * @details This will compute the particle velocities as v = mu * V(Xp) where mu is the particle mobility and V(Xp) is
 * the interpolation of m_velocityFunction to the particle position.
 */
virtual void
interpolateVelocities();

Here, the calling sequence is such that the mobilities must be interpolated first, and then the velocity fields.

Mobility coefficient interpolation

The mobility coefficient of a particle is usually interpolated directly, i.e.,

\[\mu = \mu\left(\mathbf{X}\right).\]

The other option is to compute the mobility as

\[\mu = \frac{\left(\mu\left|\mathbf{v}\right|\right)\left(\mathbf{X}\right)}{\left|\mathbf{v}\left(\mathbf{X}\right)\right|}.\]

This method ensures that the particle velocity becomes \(\mathbf{V} = \left(\mu\mathbf{v}\right)\left(\mathbf{X}\right)\).

Tip

One can switch between the two interpolation methods in the ItoSolver run-time input options.

Diffusion coefficient interpolation

Interpolation of the diffusion coefficient is always done using an interpolation method

\[D = D\left(\mathbf{X}\right).\]

The function signature is

/**
 * @brief Interpolate the diffusion field to the particle positions.
 * @details This computes D_p = Df(X_p) where Df is the diffusion field on the mesh.
 */
virtual void
interpolateDiffusion();

Particle intersections

It will happen that particles occasionally hit the embedded boundary or leave through the domain sides. In this case one might want to keep the particles in separate data holders rather than discard them. ItoSolver supplies several functions for transferring the particles to separate data containers when they intersect the EB or domain. The most relevant function is

virtual void
intersectParticles(
  const EBIntersection                                               a_ebIntersection,
  const bool                                                         a_deleteParticles,
  const std::function<void(ParticleSoA<ItoParticle>&, std::size_t)>& a_nonDeletionModifier =
    [](ParticleSoA<ItoParticle>&, std::size_t) -> void {
    return;
  });

Here, EBIntersection is just an enum for putting logic into how the intersection is computed. Valid options are EBIntersection::Bisection and EBIntersection::Raycast. These algorithms are discussed in Boundary interaction. The flag a_deleteParticles specifies if the original particles should be deleted when populating the other particle containers (again, see Boundary interaction).

After calling intersectParticles, the particles that crossed the EB or domain walls are available through the getParticles routine, see ItoSolver and can then be parsed separately by user code.

Computing time steps

While ItoSolver has no fundamental requirement on the time steps that can be used, several functions are available for computing various types of drift and diffusion related time steps.

Important

All time step calculations below are imposed on the particles and not on the mesh variables.

Advective time step

The drift time step routines are implemented such that one restricts the time step such that the fastest particle does not move more than a specified number of grid cells. This routine is implemented as

/**
 * @brief Compute advection time step dt = dx/vMax where vMax is the largest velocity component of the particle.
 * @return Computed advective dt
 */
virtual Real
computeAdvectiveDt() const;

which returns a CFL-like condition

\[\Delta t = \frac{\Delta x}{\textrm{max}(\left|v_x\right|, \left|v_y\right|, \left|v_z\right|)}.\]

Diffusive time step

The signatures for the diffusion time step are similar to the ones for drift:

/**
 * @brief Compute the diffusive dt. This computes dt = dx*dx/(2*SpaceDim*D) for all particles
 * @return Computed diffusive dt
 */
virtual Real
computeDiffusiveDt() const;

which returns a CFL-like condition

\[\Delta t = \frac{\Delta x^2}{2dD},\]

where \(d\) is the spatial dimension and \(D\) is the particle diffusion coefficient.

Advective-diffusive time step

A combination of the advection and diffusion time step routines also exists as

/**
 * @brief Compute a time step for the advance -- this calls the level function.
 * @details This computes the time step differently whether or not diffusion and advection are active. The Ito
 * particle model does not have a fundamental time step limitation, so these limits "replicate" the time step
 * selections in a 1D fluid model. If we only use advection advection the time step is computed as dt = dx/sum(|V_i|)
 * = dtA. If only diffusion is active the time step is computed as dt = (dx*dx)/(2*SpaceDim*D) = dtD. If both
 * advection and diffusion are active the time step is computed as dt = 1/(1/dtA + 1/dtD).
 * @return Computed dt
 */
virtual Real
computeDt() const;

This time step limitation is inspired by fully explicit and non-split fluid models, and is calculated as

\[\Delta t = \frac{1}{\frac{\Delta x}{\left|v_x\right| + \left|v_y\right| + \left|v_z\right|} + \frac{\Delta x^2}{2dD}}.\]

Superparticle management

It can occasionally be necessary to merge or split computational particles. This occurs in, e.g., plasma simulations where chemical reactions lead to exponential growth of particles. ItoSolver handles superparticles via a configurable merger functor selected at parse time through ItoSolver.merge_algorithm; the user can also supply a custom functor through setParticleCellMerger. The entry point for splitting and merging is in all cases

/**
 * @brief Make superparticles for a full container -- this is the AMR version that users will usually call.
 * @param[in] a_container        Which container to repartition into new superparticles
 * @param[in] a_particlesPerCell Target number of particles per cell
 */
virtual void
makeSuperparticles(WhichContainer a_container, int a_particlesPerCell);

Calling this function will merge/split the particles.

Important

Most merging algorithms are performed within each grid cell, and particles must therefore be sorted by their cell index (organizeParticlesByCell) before calling the merging routine. The exceptions are nn_pair_tree, nn_pair_onecell, nn_pair_hash, kd_carve, kd_patch, and kd_skin_nn, which are distributed AMR-level merges dispatched over the whole container rather than cell by cell. All of these except kd_patch match particles across patch and rank boundaries and therefore require that a particle ghost halo has been filled; kd_patch is patch-local and instead requires that no ghost halo is present.

In order to specify the merging algorithm the user must set the ItoSolver.merge_algorithm to one of the following:

• none - No particle merging/splitting is performed.

• equal_weight_kd Use a kd-tree with bounding volume hierarchies to partition and split/merge the particles. This conserves the particle center-of-mass.

• reinitialize Re-initialize the particles in each grid cell, ensuring that weights are as uniform as possible.

• reinitialize_bvh Re-initialize the particles in each node of a kd-tree. Weights are as uniform as possible.

• nn_sfc Reach the target particle count by space-filling-curve nearest-neighbour clustering: when there are more particles than the target the nearest neighbours (along a Hilbert curve) are merged until the target count is reached, and when there are fewer the highest-weight particles are split. This gives spatially tight groups but does not equalize the weights.

• nn_pair_tree A distributed, MPI-safe nearest-neighbour pair merge that reaches the target particle count over the whole AMR hierarchy, searching for candidates via one whole-patch PointCloudBVH per patch. Over-full cells are drained by matching each over-crowded particle with its true nearest neighbour across patch and rank boundaries (a propose/judge/verdict protocol over a particle ghost halo) and merging the pair to its weighted centroid; because a single round merges pairs, the round is repeated until every cell reaches the target, or until ItoSolver.nn_pair_max_rounds rounds have run. Under-full cells are then brought up to the target by splitting the heaviest particle into two co-located daughters (floor/ceil weights, so integer weights stay integer). Tunable through ItoSolver.nn_pair_iterate, ItoSolver.nn_pair_fallback, ItoSolver.nn_pair_max_rounds and ItoSolver.nn_pair_max_cell_dist.

• nn_pair_onecell The same distributed nearest-neighbour pair merge and drain/split protocol as nn_pair_tree, but candidates are found via one PointCloudBVH per occupied grid cell instead of one per patch: a query only ever searches its own cell and its Moore-adjacent neighbours, so the merge distance is structurally fixed at Chebyshev cell distance 1 and ItoSolver.nn_pair_max_cell_dist does not apply. Tunable through ItoSolver.nn_pair_iterate, ItoSolver.nn_pair_fallback and ItoSolver.nn_pair_max_rounds.

• nn_pair_hash The same distributed nearest-neighbour pair merge and drain/split protocol as nn_pair_tree, but candidates are found via one PointCloudHashGrid (a uniform spatial hash grid) per patch instead of a PointCloudBVH. Identical tunables and behaviour to nn_pair_tree (ItoSolver.nn_pair_iterate, ItoSolver.nn_pair_fallback, ItoSolver.nn_pair_max_rounds, ItoSolver.nn_pair_max_cell_dist); only the per-patch spatial-index backend differs.

• kd_carve A distributed, MPI-safe whole-patch merge built around a spatial partition (“kd tree”) rather than a nearest-neighbour graph: each patch splits its local-plus-ghost particles purely by position, bisecting the longest axis and never snapping to the grid, so particles near a cell face can merge across it – at the count median while a node is still larger than ItoSolver.kd_split_weight_leaf_dx cell widths, and at the weight median once it is smaller, so that the resulting super-particles come out with comparable weights. How far splitting goes is set by a live per-cell quota: every group becomes exactly one super-particle at its weighted centroid, so the number of groups centred in a cell is that cell’s post-merge population, and any split that would push a cell past the target count is refused. A group entirely clear of any patch/rank boundary merges immediately with no communication; a group that touches one is resolved by a single, non-iterative “z-buffer carve” – competing groups from neighbouring patches are ranked by a deterministic key and the tightest one wins each contested particle. Unlike the nn_pair_* family there is no drain loop and no nn_pair_iterate/nn_pair_fallback/nn_pair_max_cell_dist equivalent to tune; the ghost width is fixed at 1. The maximum per-axis extent of a mergeable group is fixed at one cell width – a physical safety bound applied as a filter on the finished partition, not a tunable, and matching the hardcoded width-1 ghost halo. Tunable through ItoSolver.kd_split_weight_leaf_dx (count-median/weight-median crossover, in cell widths; 0 disables the weight median).

• kd_patch The same whole-patch kd-tree build and per-cell quota as kd_carve, but with no boundary tier: no particle ghost halo is filled, every group merges regardless of boundary exposure, and no particle is ever contested, so each patch reduces only the particles it owns. Cheaper and communication-free, at the cost of no coordination across patch boundaries.

• kd_skin_nn The same whole-patch kd-tree build and per-cell quota as kd_carve, but the boundary tier is the nearest-neighbour pair merge instead of the carve arbitration. The particles are split across two containers: every group holding no ghost particle is committed into one of them (holding no ghost is the whole safety condition – such a group’s members are all resident in this patch, and any neighbouring patch that draws one of them into a group of its own sees it as a ghost and is disqualified by the same test), and the contested remainder – the skin – is drained from the other by the same algorithm as nn_pair_onecell, whose Chebyshev-1 search radius matches the width-1 ghost halo exactly. Unlike kd_carve boundary exposure is never consulted, so an exposed but uncontested group merges locally, which is what makes the skin small. The skin tier drains each cell to what is left of its target after the already-merged interior results are counted against it, so the two tiers share one budget. Tunable through ItoSolver.kd_split_weight_leaf_dx and, for the skin tier, ItoSolver.nn_pair_iterate, ItoSolver.nn_pair_fallback and ItoSolver.nn_pair_max_rounds; ItoSolver.nn_pair_max_cell_dist does not apply.

• external Use an externally injected particle merging algorithm. In order to use this feature the user must supply one through setParticleCellMerger.

The user can set the merging algorithm through the input script (see Input options), or supply one externally by setting the merge algorithm to external. In addition, the user must first supply a particle merging function:

/**
 * @brief Set the user-supplied per-cell particle merger used by merge_algorithm = external.
 * @details The built-in merge methods build their own per-cell merger internally; this one is used
 * only when merge_algorithm = external, applied cell-by-cell by makeSuperparticles().
 * @note The merger runs on ItoMergeParticle, not ItoParticle: every merge in this solver operates on
 * the reduced particle, and a per-cell merger sees exactly the columns a merge is entitled to change
 * (position, weight, energy). Out-of-tree mergers written against ItoParticle break at compile time.
 * @param[in] a_particleCellMerger Per-cell particle merger.
 */
virtual void
setParticleCellMerger(const ParticleManagement::ParticleMerger<ItoMergeParticle>& a_particleCellMerger) noexcept;

In the code above, ParticleManagement::ParticleMerger<P> is an alias:

template <class P, class Traits = ParticleTraits<P>>
using ParticleMerger = std::function<
  void(ParticleSoA<P, Traits>& a_particles, const CellInfo& a_cellInfo, const int a_numTargetParticles)>;

Tip

ItoSolver uses the kd-tree implementation from Merging and splitting particles and partitioners for splitting the particles into two subsets with equal weights.

Example transport kernel

Transport kernels for the particles within ItoSolver will typically be imposed externally by the user through a TimeStepper subclass that advances the particles. For completeness, we here include a simple transport kernel for the ItoSolver which simply consists of a drift-diffusion kick. ItoParticle is a Struct-of-Arrays payload (see Particles), so the kernel operates on a ParticleSoA<ItoParticle> leaf and addresses each particle by index; the position is a container-owned column accessed through position(i)/setPosition(i, ...), while the interpolated velocity, diffusion coefficient, and old position are payload columns accessed through get<...>(i). The loop below is the per-patch inner kernel and is run inside the usual level/patch iteration (see Particles):

// One grid patch. The velocity columns (vx/vy/vz) have already been filled
// by ItoSolver::interpolateVelocities().
ParticleSoA<ItoParticle>& leaf = particles[lvl][dit()];

for (std::size_t i = 0; i < leaf.size(); i++) {
   const RealVect      x = leaf.position(i);
   const RealVect      v = RealVect(D_DECL(leaf.get<&ItoParticle::vx>(i),
                                           leaf.get<&ItoParticle::vy>(i),
                                           leaf.get<&ItoParticle::vz>(i)));
   const ParticleReal& D = leaf.get<&ItoParticle::diffusion>(i);

   // Store the old position in the payload's old-position columns.
   D_TERM(leaf.get<&ItoParticle::old_x>(i) = x[0];,
          leaf.get<&ItoParticle::old_y>(i) = x[1];,
          leaf.get<&ItoParticle::old_z>(i) = x[2];);

   // Drift-diffusion kick.
   leaf.setPosition(i, x + v * a_dt + sqrt(2.0 * D * a_dt) * this->randomGaussian());
}

The function randomGaussian implements a diffusion hopping and returns a 2D/3D dimensional vector with values drawn from a normal distribution with standard width of one and mean value of zero. The implementation uses the random number generators in Random numbers.

I/O

Plot files

For a complete list of available plot variables, see Input options.

Input options

Several input options are available for configuring the run-time configuration of ItoSolver, which are listed in Listing 27.

Listing 27 Input options for the ItoSolver class. All options are run-time configurable.

# ====================================================================================================
# ItoSolver class options
# ====================================================================================================
ItoSolver.verbosity                   = -1             ## Class verbosity
ItoSolver.merge_algorithm             = kd_carve       ## Particle merging algorithm. One of 'reinitialize', 'equal_weight_kd', 'reinitialize_bvh', 'nn_sfc', 'nn_pair_tree', 'nn_pair_onecell', 'nn_pair_hash', 'kd_carve', 'kd_patch', 'kd_skin_nn', 'none', or 'external'
ItoSolver.particles_per_cell          = 32             ## Target computational particles per cell (one value, or one per level)
ItoSolver.regrid_superparticles       = solver         ## Merge run during regrids: 'solver' (use merge_algorithm), 'none', or any merge_algorithm selector
ItoSolver.nn_pair_iterate             = true           ## nn_pair_tree/nn_pair_onecell/nn_pair_hash: iterate the local tier to convergence within each round
ItoSolver.nn_pair_fallback            = 1              ## nn_pair_tree/nn_pair_onecell/nn_pair_hash: fallback candidates per query
ItoSolver.nn_pair_max_cell_dist       = 1              ## nn_pair_tree/nn_pair_hash: max merge distance in cells (Chebyshev); unused by nn_pair_onecell (fixed at 1)
ItoSolver.nn_pair_max_rounds          = 3              ## nn_pair_tree/nn_pair_onecell/nn_pair_hash: max drain rounds per merge (stops early on convergence)
ItoSolver.kd_split_weight_leaf_dx     = 1.0            ## kd_carve/kd_patch: leaf size (cell widths) at/below which splitting switches from count-median to weight-median (0 = never)
ItoSolver.plt_vars                    = phi vel dco    ## 'phi', 'vel', 'dco', 'part', 'eb_part', 'dom_part', 'src_part', 'energy_density', 'energy'
ItoSolver.intersection_alg            = bisection      ## Intersection algorithm for EB-particle intersections.
ItoSolver.bisect_step                 = 1.E-4          ## Bisection step length for intersection tests
ItoSolver.normal_max                  = 5.0            ## Maximum value (absolute) that can be drawn from the exponential distribution.
ItoSolver.redistribute                = false          ## Turn on/off redistribution.
ItoSolver.blend_conservation          = false          ## Turn on/off blending with nonconservative divergenceo
ItoSolver.checkpointing               = particles      ## 'particles' or 'numbers'
ItoSolver.ppc_restart                 = 32             ## Particles per cell redrawn when restarting from a 'checkpointing = numbers' checkpoint
ItoSolver.irr_ngp_deposition          = true           ## Force irregular deposition in cut cells or not
ItoSolver.irr_ngp_interp              = true           ## Force irregular interpolation in cut cells or not
ItoSolver.mobility_interp             = direct         ## How to interpolate mobility, 'direct' or 'velocity', i.e. either mu_p = mu(X_p) or mu_p = (mu*E)(X_p)/E(X_p)
ItoSolver.plot_deposition             = cic            ## Cloud-in-cell for plotting particles.
ItoSolver.deposition                  = cic            ## Deposition type.
ItoSolver.deposition_cf               = transition     ## 'interp', 'halo', 'halo_ngp', 'transition'.

Plot file variables

Plot variables are specified using ItoSolver.plt_vars, see Plot files. To add a variable to HDF5 output files, one can modify the ItoSolver.plt_vars input variable to include, e.g., the following variables:

• \(\phi\), i.e. the deposited particle weights (ItoSolver.plt_vars = phi)

• \(\mathbf{v}\), the advection field (ItoSolver.plt_vars = vel).

• \(D\), the diffusion coefficient (ItoSolver.plt_vars = dco).

Particle-mesh configuration

To specify the mobility interpolation, use ItoSolver.mobility_interp. Valid options are direct and velocity, see Particle interpolation.

Deposition and coarse-fine deposition (see Particle-mesh) are controlled using the flags

• ItoSolver.deposition for the base deposition scheme. Valid options are ngp, cic, and tsc.

• ItoSolver.deposition_cf for the coarse-fine deposition strategy. Valid options are interp, halo, or halo_ngp.

To modify the deposition scheme in cut-cells, one can enforce NGP interpolation and deposition through

• ItoSolver.irr_ngp_deposition for enforcing NGP deposition. Valid options are true or false.

• ItoSolver.irr_ngp_interp for enforcing NGP interpolation. Valid options are true or false.

Checkpoint-restart

Available input options for the ItoSolver are listed below:

Example application(s)

Example applications that use ItoSolver are found in

• $DISCHARGE_HOME/Physics/BrownianWalker, see Brownian walker.

• $DISCHARGE_HOME/Physics/ItoKMC, see Îto-KMC plasma model.