Îto diffusion

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

(9)\[\Delta\mathbf{X} = \left(\mathbf{V} + \nabla D\right)\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.

The $\nabla D$ drift correction

Eq. 9 is an Îto-interpretation update, and the $\nabla D$ term in its drift is what makes it transport the same diffusive flux as the fluid solvers. Without that term the update is $\Delta\mathbf{X} = \mathbf{V}\Delta t + \sqrt{2D\Delta t}\mathbf{W}$, whose Fokker-Planck equation is

\[\frac{\partial n}{\partial t} = -\nabla\cdot\left(\mathbf{V}n\right) + \nabla^2\left(Dn\right),\]

i.e. it transports the flux $\mathbf{V}n - n\nabla D - D\nabla n$. CdrSolver integrates

\[\frac{\partial n}{\partial t} = -\nabla\cdot\left(\mathbf{V}n\right) + \nabla\cdot\left(D\nabla n\right),\]

i.e. the flux $\mathbf{V}n - D\nabla n$. The two differ by $\nabla\cdot\left(n\nabla D\right)$ and coincide exactly where $D$ is constant in space, which is why the discrepancy is invisible in any test with a constant diffusion coefficient. Adding $\nabla D$ to the drift, as Eq. 9 does, recovers the fluid form; equivalently, it integrates the stochastic equation in the Hänggi-Klimontovich (anti-Îto) interpretation.

The correction is on by default and can be turned off with

ItoSolver.diffusion_grad_drift = true

in which case the solver reverts to the plain Îto update and transports $\mathbf{V}n - \nabla\left(Dn\right)$. ItoSolver computes $\nabla D$ from the mesh diffusion coefficient and interpolates it to the particle positions; it is the time stepper that adds $\Delta t\nabla D$ to the particle displacement. Setting ItoSolver.plt_vars to include grad_dco writes $\nabla D$ to the plot files.

Note

The correction is defined for a scalar, isotropic $D$. It is also only meaningful when $D$ is a mesh field: a solver that assigns particle diffusion coefficients parametrically from the particle energy has no mesh field whose gradient means anything. Finally, the time step estimates in Computing time steps do not account for this term – see the computeDt documentation for why.

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:

Mobility, diffusion coefficient, energy (not currently used), and a holder for a scratch scalar storage.
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:

Move particles with a microscopic drift-diffusion model.
Compute particle intersection with embedded boundaries and domain edges.
Deposit particles and other particle types on the mesh.
Interpolate velocities and diffusion coefficients to the particle positions.
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:

Mobility.
Diffusion coefficient.
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:

Initial particles specified from a container (ParticleSoA<ItoParticle>) of particles.
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
 * @note The grad(D) drift correction (ItoSolver.diffusion_grad_drift) is not accounted for here. These
 * estimates are built from the particle velocity and diffusion columns, and grad(D) lives on neither -- it is
 * interpolated at the point of use and added to the displacement by the time stepper. The omission is
 * deliberate rather than an oversight: the correction displaces a particle by dt*grad(D), which for a streamer
 * in atmospheric air (D ~ 0.1 m^2/s varying over ~10 um, so |grad D| ~ 1e4 m/s) is about two percent of the
 * diffusion hop sqrt(2*D*dt) at dt = 1 ps, and the hop is what the diffusive limit already bounds. A model
 * that made grad(D) comparable to the drift velocity would make these estimates optimistic.
 */
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
 * @note The grad(D) drift correction (ItoSolver.diffusion_grad_drift) is not accounted for here. These
 * estimates are built from the particle velocity and diffusion columns, and grad(D) lives on neither -- it is
 * interpolated at the point of use and added to the displacement by the time stepper. The omission is
 * deliberate rather than an oversight: the correction displaces a particle by dt*grad(D), which for a streamer
 * in atmospheric air (D ~ 0.1 m^2/s varying over ~10 um, so |grad D| ~ 1e4 m/s) is about two percent of the
 * diffusion hop sqrt(2*D*dt) at dt = 1 ps, and the hop is what the diffusive limit already bounds. A model
 * that made grad(D) comparable to the drift velocity would make these estimates optimistic.
 */
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
 * @note The grad(D) drift correction (ItoSolver.diffusion_grad_drift) is not accounted for here. These
 * estimates are built from the particle velocity and diffusion columns, and grad(D) lives on neither -- it is
 * interpolated at the point of use and added to the displacement by the time stepper. The omission is
 * deliberate rather than an oversight: the correction displaces a particle by dt*grad(D), which for a streamer
 * in atmospheric air (D ~ 0.1 m^2/s varying over ~10 um, so |grad D| ~ 1e4 m/s) is about two percent of the
 * diffusion hop sqrt(2*D*dt) at dt = 1 ps, and the hop is what the diffusive limit already bounds. A model
 * that made grad(D) comparable to the drift velocity would make these estimates optimistic.
 */
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_method; 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. Both kd_cell and reinitialize are of this kind. The exceptions are kd_patch, kd_amr and nn_amr, which are distributed merges dispatched over the whole container rather than cell by cell. kd_amr and nn_amr 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.

The merging algorithm is selected in two steps: ItoSolver.merge_method names the scope the merge groups particles over, and a small set of kd_*/nn_* specifiers say how it groups within that scope. Splitting the two apart is what keeps the selector list short: kd_patch and kd_amr run the same tree build over a different scope, and the three nearest-neighbour backends differ only in the spatial index they search.

ItoSolver.merge_method takes one of:

none - No particle merging/splitting is performed.
kd_cell Partition each grid cell’s particles with a kd tree and reduce every leaf to one particle. Cell-scoped, so a leaf never reaches past the cell it belongs to and the per-cell density is preserved exactly. Reads ItoSolver.kd_partition and ItoSolver.kd_placement.
kd_patch The same tree build, but one tree per patch rather than per cell: nodes are bisected on the longest axis and are never snapped to the grid, so a leaf can straddle a cell face and merge particles across it. How far splitting goes is set by a live per-cell quota – the number of leaves centred in a cell is that cell’s post-merge population, and a split that would push a cell past ItoSolver.particles_per_cell is refused. No particle ghost halo is filled and no particle is ever contested, so each patch reduces only the particles it owns: communication-free, at the cost of no coordination across patch boundaries.
kd_amr The same whole-patch tree build and per-cell quota as kd_patch, but resolved across patch and rank boundaries, so a leaf may draw members owned by a neighbour. ItoSolver.kd_amr_boundary picks how that is done:
- carve – a single, non-iterative “z-buffer carve”: competing leaves from neighbouring patches are ranked by a deterministic key and the tightest one wins each contested particle. There is no drain loop and nothing to tune beyond the tree itself; the ghost width is fixed at 1, and the maximum per-axis extent of a mergeable leaf is fixed at one cell width to match it.
- nn – every leaf holding no ghost particle is committed immediately (holding no ghost is the whole safety condition: such a leaf’s members are all resident in this patch, and any neighbouring patch that draws one of them into a leaf of its own sees it as a ghost and is disqualified by the same test), and the contested remainder – the skin – is drained by the nearest-neighbour pair merge, whose Chebyshev-1 search radius matches the width-1 ghost halo exactly. Boundary exposure is never consulted, so an exposed but uncontested leaf merges locally, which is what keeps the skin small. The two tiers share one per-cell budget: the skin drains each cell to what is left of its target after the interior results are counted against it. Also reads ItoSolver.nn_iterate, ItoSolver.nn_fallback and ItoSolver.nn_max_rounds; ItoSolver.nn_max_cell_dist does not apply.
nn_amr A distributed, MPI-safe nearest-neighbour pair merge that reaches the target particle count over the whole AMR hierarchy. 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_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). ItoSolver.nn_search picks the spatial index the candidate search is backed by:
- tree – one whole-patch PointCloudBVH per patch.
- hash – one whole-patch PointCloudHashGrid (a uniform spatial hash grid) per patch. Identical behaviour to tree; only the backend differs.
- onecell – one PointCloudBVH per occupied grid cell. 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_max_cell_dist does not apply.
Tunable through ItoSolver.nn_iterate, ItoSolver.nn_fallback, ItoSolver.nn_max_rounds and (for tree/hash) ItoSolver.nn_max_cell_dist.
reinitialize Discard each cell’s spatial information entirely and rebuild it: the cell’s physical particles are divided into near-equal integer weights and placed at uniformly random positions in the cell. Cell-scoped. Requires integer weights.
external Use an externally injected particle merging algorithm. In order to use this feature the user must supply one through setParticleCellMerger.

Every kd_* method reads the same two specifiers, because every one of them builds a kd tree and reduces each leaf to one particle – only the scope of the tree differs.

ItoSolver.kd_partition – how a node is divided into two children. weight divides it so the two halves carry as nearly equal a weight as possible, which keeps the merged weights equal; the per-cell build divides a particle across the split when it must, so it can create particles during the build, while the patch and AMR builds never do and so cannot reach exact equality. count divides on particle count and never divides a particle, so it creates none and can never overshoot the target mid-build, but it leaves the merged weights uneven – measurably so under repeated merging, where the effective particle count per cell falls well below the target. hybrid is the two combined by node size: the count median while a node is wider than ItoSolver.kd_hybrid_leaf_dx cell widths, the weight median once it is narrower. That is the rule the patch and AMR scopes want, since a split of a node spanning several cells apportions leaves between cells, which is a question of counts, while a narrow node is where the weight median earns its keep. weight and count are the two limits of the same rule – a crossover above every node’s size and one below it – which is why one setting serves all three scopes. weight_capped is weight with a pre-pass that divides every particle heavier than a leaf’s target weight into pieces before the tree is built; a particle heavier than that target cannot be equalized away by any partition, since it alone sets a floor on whichever leaf it lands in. The pre-pass matters most for kd_patch, whose build may not divide a particle mid-tree at all and so cannot otherwise reach equal weights.

Warning

weight_capped is refused with merge_method = kd_amr. That scope’s boundary tier arbitrates contested particles by id and resolves a winner’s weight and position by looking the id up in a table built before the merge divides anything; the cap’s pieces inherit their parent’s id, so the lookup would return the parent’s full pre-cap weight for every piece and double-count its mass across the patch boundary. Lifting the restriction needs the pieces to carry fresh ids and the id table to be rebuilt after the cap.
ItoSolver.split_placement – where the pieces of a split particle go. A split divides one particle’s weight among several pieces; center puts every piece at the parent’s position, jitter draws each from the local mean interparticle spacing around the parent, and cell draws each uniformly in the owning cell. All three keep every piece inside the parent’s own cell, so the merge still moves no mass across a cell face, and all three fall back to center in a cut cell.

center perturbs nothing spatially, but it makes the pieces indistinguishable, and a placement that inherits member positions (sample) then reproduces them. Over repeated merges the number of distinct positions a cell holds falls well below its particle count – measured at a median of 6 distinct positions per 16 particles after 500 merges. That loss is invisible in the weight statistics, which stay perfectly uniform, while the cell’s effective spatial sampling degrades by the same factor. jitter restores it, at the cost of moving mass by less than the distance between neighbouring particles.

Note

The jitter kernel is truncated to the cell, not clamped and not reflected. Clamping stacks every out-of-range draw onto the wall. Reflecting folds the overhanging part of the kernel back onto the part that stayed, doubling the density within half a bandwidth of each wall – and because the fold happens independently per direction, a parent near a corner is folded once per direction and the corner receives $2^D$ times the density it should. Truncation is the correct conditional distribution and has no directional structure.

ItoSolver.kd_placement – where the particle a leaf reduces to is placed. centroid puts it at the leaf’s weighted centroid, which conserves the leaf’s centre of mass exactly. sample puts it at one of the leaf’s own particles, drawn with probability proportional to weight. random puts it at a uniformly random point in the leaf’s bounding box (falling back to the centroid in cut cells, to keep the particle inside the embedded boundary).

The centroid keeps a leaf’s mean position and discards its spread, so it contracts the sub-cell distribution a little at every merge and, since the leaves of a locally uniform population are near enough equal-volume sub-boxes, drives the particles onto a regular sub-cell lattice that is the same in every grid cell and therefore coherent across the whole domain. sample is distribution-preserving – a weight-proportional draw from a partition of a sample is itself a sample of the same distribution – so neither the lattice nor the contraction occurs, at the cost of conserving the per-leaf centre of mass in expectation rather than exactly. random also avoids the lattice, but draws from the leaf’s bounding box rather than the leaf itself, which is not quite the same distribution. Prefer sample wherever sub-cell position matters – near a refinement boundary, or near an embedded boundary, where a centroid is not a position any particle occupied.

Note

For kd_cell the leaf lies inside one cell, so every placement keeps the leaf’s weight in the cell it came from and the per-cell density is preserved exactly. For kd_patch and kd_amr the leaves are deliberately not snapped to the grid, so a leaf straddling a cell face may place its merged particle on either side whichever placement is chosen – that is what lets those scopes merge across a face in the first place, and it is why their per-cell particle count is a target rather than a guarantee. In all three, a placement that would land inside an embedded boundary falls back to the leaf’s centroid.

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_method = external.
 * @details The built-in merge methods build their own per-cell merger internally; this one is used
 * only when merge_method = 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.
Most are run-time configurable; seed and diffusion_grad_drift are read once at setup.

# ====================================================================================================
# ItoSolver class options
# ====================================================================================================
ItoSolver.verbosity                 = -1             ## Class verbosity
ItoSolver.particles_per_cell        = 32             ## Target computational particles per cell (one value, or one per level)
ItoSolver.merge_method              = kd_cell        ## Merge scope. One of 'kd_cell', 'kd_patch', 'kd_amr', 'nn_amr', 'reinitialize', 'none', or 'external'
ItoSolver.regrid_superparticles     = solver         ## Merge run during regrids: 'solver' (use merge_method), 'none', or any merge_method selector
ItoSolver.kd_partition              = weight         ## kd_cell/kd_patch/kd_amr: split rule. 'weight', 'count', 'hybrid', or 'weight_capped' (not with kd_amr)
ItoSolver.kd_hybrid_leaf_dx         = 1.0            ## kd_partition = hybrid: crossover node size, in cell widths
ItoSolver.kd_placement              = centroid       ## kd_cell/kd_patch/kd_amr: leaf placement. 'centroid', 'sample', or 'random'
ItoSolver.split_placement           = center         ## kd_*: where a split particle's pieces go. 'center', 'jitter', or 'cell'
ItoSolver.kd_amr_boundary           = carve          ## kd_amr: boundary tier. 'carve' or 'nn'
ItoSolver.nn_search                 = tree           ## nn_amr: neighbour search. 'tree', 'hash', or 'onecell'
ItoSolver.nn_iterate                = true           ## nn_amr/kd_amr_boundary = nn: iterate the local tier within each round
ItoSolver.nn_fallback               = 1              ## nn_amr/kd_amr_boundary = nn: fallback candidates per query
ItoSolver.nn_max_cell_dist          = 1              ## nn_amr: max merge distance in cells. Unused by nn_search = onecell
ItoSolver.nn_max_rounds             = 3              ## nn_amr/kd_amr_boundary = nn: max drain rounds per merge
ItoSolver.plt_vars                  = phi vel dco    ## 'phi', 'vel', 'dco', 'grad_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.diffusion_grad_drift      = true           ## Add grad(D) to the drift so the particles transport D*grad(n), as CdrSolver does
ItoSolver.checkpointing             = particles      ## 'particles' or 'numbers'
ItoSolver.ppc_restart               = 32             ## Particles per cell redrawn when restarting from a 'checkpointing = numbers' checkpoint
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_irreg          = mirror         ## Cut-cell deposition: 'native', 'ngp', 'mirror', 'redistribute', 'redistribute_blended'
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.

How the cut cells are treated when depositing is selected with a single flag,

ItoSolver.deposition_irreg, with valid options
- native – deposit as-is, with no cut-cell treatment. A cut cell then holds $\kappa n$.
- ngp – put the particle’s entire cloud in its own cell when that cell is cut.
- mirror – even extension of the density about the embedded boundary. Under this option a cut cell holds $n$, not $\kappa n$, unlike every other option here; consumers of $\phi$ must opt into that change of meaning knowingly. Incompatible with ItoSolver.deposition = ngp.
- redistribute – hybrid divergence, with the cut-cell mass excess redistributed to the neighbours.
- redistribute_blended – as redistribute, but blended with the non-conservative divergence.

These are mutually exclusive answers to one question, which is why they are one selector rather than independent flags: mirror combined with redistribution would apply a $1/\kappa$ correction to a field that no longer needs one.

Note

mirror is experimental. It is verified to add mass only where it should and never to remove any, but its acceptance suite is not finished, and the change of meaning above – a cut cell holding $n$ rather than $\kappa n$ – has not been audited against every consumer of $\phi$. Use it for deliberate experiments, not production runs.

Cut-cell interpolation is separate and remains a boolean:

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.

Îto diffusion

The \(\nabla D\) drift correction

ItoParticle

ItoSolver

ItoSpecies

Supplying initial data

Particle containers

Remapping particles

Particle deposition

AMR synchronization after a deposit

Deposition of other quantities

Particle interpolation

Mobility coefficient interpolation

Diffusion coefficient interpolation

Particle intersections

Computing time steps

Advective time step

Diffusive time step

Advective-diffusive time step

Superparticle management

Example transport kernel

I/O

Plot files

Input options

Plot file variables

Particle-mesh configuration

Checkpoint-restart

Example application(s)