Î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 derives from GenericParticle as follows:

class ItoParticle : public GenericParticle<5, 3>

From the signature one can see that ItoParticle contains a number of extra class Real and RealVect class members. These extra fields are used for storing the following information in the particle:

  1. Particle weight, mobility, diffusion coefficient, energy (not currently used), and a holder for a scratch 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 the 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 positons.

  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 to simply to improve flexibility when assigned 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 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 indicates 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 list (List<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
*/
List<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
*/
List<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 then 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.
*/
virtual ParticleContainer<ItoParticle>&
getParticles(const 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(const 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 Generic deposition function which deposits a particle field onto the mesh using a specified deposition method. 
  @details The template parameters indicate the particle type and quantity to be deposited. The second template parameter must be a pointer to a member function in 
  the particle class with signature 'const Real& P::function() const'. E.g. 'const Real& P::mass() const' which is the default quantity to be deposited. 
  @param[out] a_phi                  Mesh data -- must have exactly one component. 
  @param[in]  a_particles            Particles to be deposited
  @param[in]  a_deposition           Deposition method
  @param[in]  a_coarseFineDeposition Coarse-fine deposition strategy
*/
template <class P, class Ret, Ret (P ::*MemberFunc)() const>
void
depositParticles(EBAMRCellData&             a_phi,
                 ParticleContainer<P>&      a_particles,
                 const DepositionType       a_deposition,
                 const CoarseFineDeposition a_coarseFineDeposition) const;

This version permits the user to select any particle container a_particles and deposit them onto some pre-allocated mesh storage a_phi. Note that the template type P does not need to be ItoParticle, although this is the most common use case.

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 onto mesh. 
  @details This will deposit the mass (i.e., computational weight) "bulk" particles into m_phi. 
  @note Calls the other version with a_container = WhichContainer::Bulk
*/
virtual void
depositParticles();

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.

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
*/
virtual void
depositConductivity(EBAMRCellData& a_phi, ParticleContainer<ItoParticle>& a_particles) const;

/*!
  @brief Deposit diffusivity (i.e. mass*D/volume)
  @details This deposits mass*diffusion (not multiplied by charge)
  @param[out] a_phi       Mesh data
  @param[in]  a_particles Particle data
  @note Calls the other versions with a_deposition = m_deposition
*/
virtual void
depositDiffusivity(EBAMRCellData& a_phi, ParticleContainer<ItoParticle>& a_particles) 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 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 signatures 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

/*!
  @brief Do boundary intersection tests.
  @details This will intersect the particles in the "bulk" particles data holder with the domain faces and EBs. If a particle crossed the EB it will
  be put into the "EB" particle data holder and likewise for the particles that crossed the domain side. 
  @param[in] a_ebIntersection      Enum for switching between various types of intersection tests. 
  @param[in] a_deleteParticles     If true, the origin particle will also be removed from the bulk particle data holder.
  @param[in] a_nonDeletionModifier Optional input argument for letting the user manipulate particles that were intersected but not deleted
  @note This will call the other version with WhichContainer::Bulk, WhichContainer::EB, and WhichContainer::Domain.
*/
virtual void
intersectParticles(
  const EBIntersection                    a_ebIntersection,
  const bool                              a_deleteParticles,
  const std::function<void(ItoParticle&)> a_nonDeletionModifier = [](ItoParticle&) -> void {
    return;
  });

Here, EbIntersection is a 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. 
*/
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
*/
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). 
*/
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 can currently handle superparticles through several internal functions, and is also equipped with an interface in which the user can inject an external particle-handling routine. The function for splitting and merging the particles 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(const WhichContainer a_container, const int a_particlesPerCell);

Calling this function will merge/split the particles.

Important

Particle merging is currently performed within each grid cell, and particles must therefore be sorted by their cell index before calling the merging routine.

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.

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

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 particle merger. This will get called when merging particles using makeSuperparticles.
  @param[in] a_particleMerger Particle merger
*/
virtual void
setParticleMerger(const ParticleManagement::ParticleMerger<ItoParticle>& a_particleMerger) noexcept;

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

/*!
  @brief Concept for splitting/merging particles
  @param[inout] a_particles Particles to be merged/split
  @param[in] a_cellInfo Cell info
  @param[in] a_numTargetParticles Number of target particles
*/
template <class P>
using ParticleMerger = std::function<
  void(List<P>& a_particles, const CellInfo& a_cellInfo, const int a_numTargetParticles)>;

Tip

ItoSolver uses the kD-node implementation from Superparticles 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:

List<ItoParticle> particles;

for (ListIterator<ItoParticle>& lit(particles); lit.ok(); ++lit) {
   ItoParticle& p = lit();

   p.oldPosition() = p.position();
   p.position()   += p.velocity() * a_dt + sqrt(2.0*p.diffusion()*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        = equal_weight_kd    ## Particle merging algorithm. Either 'reinitialize', 'equal_weight_kd', or 'reinitialize_bvh'
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                 ## Maximum number of computational particles to generate for restarts.
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) is 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