Tracer particles

Tracer particles are particles that move along a prescribed velocity field

\[\frac{\partial\mathbf{X}}{\partial t} = \mathbf{V}\]

where \(\mathbf{X}\) is the particle position and \(\mathbf{V}\) is the particle velocity. The velocity is interpolated from a mesh-based field as

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

where \(\mathbf{v}\) is a velocity field defined on the mesh. Such particles are useful, for example, for numerical integration along field lines.

Tip

The chombo-discharge tracer particle functionality resides in $DISCHARGE_HOME/Source/TracerParticles.

TracerParticleSolver

The tracer particle solver is templated as

template <typename P>
class TracerParticleSolver

where P is the payload type used for the solver. The particles are stored Struct-of-Arrays in a ParticleContainer: the position and weight are container-owned columns, while P declares the extra payload columns. The template constraints on P are

  1. It must expose the velocity as SpaceDim scalar payload columns named D_DECL(v_x, v_y, v_z).

  2. It must have a ParticleTraits<P> specialization listing its columns (including v_x/v_y/v_z).

Users are free to provide their own payload provided that it meets these constraints. However, we also define a plug-and-play payload for the tracer-particle stepper, see TracerParticle.

Note

The TracerParticleSolver<P> API is available at https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classTracerParticleSolver.html.

TracerParticle

TracerParticle is the plug-and-play payload used by the tracer-particle stepper. It declares the interpolated velocity together with Runge-Kutta stage scratch, all as per-component SoA columns:

struct TracerParticle
{
  ParticleReal v_x = 0.0; ///< Interpolated velocity, x-component.
  ParticleReal v_y = 0.0; ///< Interpolated velocity, y-component.
#if CH_SPACEDIM == 3
  ParticleReal v_z = 0.0; ///< Interpolated velocity, z-component.
#endif

  double xk_x = 0.0; ///< Position snapshot x^k, x-component (double: position-like).
  double xk_y = 0.0; ///< Position snapshot x^k, y-component (double: position-like).
#if CH_SPACEDIM == 3
  double xk_z = 0.0; ///< Position snapshot x^k, z-component (double: position-like).
#endif

  ParticleReal k1_x = 0.0; ///< Runge-Kutta stage k1, x-component.
  ParticleReal k1_y = 0.0; ///< Runge-Kutta stage k1, y-component.
#if CH_SPACEDIM == 3
  ParticleReal k1_z = 0.0; ///< Runge-Kutta stage k1, z-component.
#endif

  ParticleReal k2_x = 0.0; ///< Runge-Kutta stage k2, x-component.
  ParticleReal k2_y = 0.0; ///< Runge-Kutta stage k2, y-component.
#if CH_SPACEDIM == 3
  ParticleReal k2_z = 0.0; ///< Runge-Kutta stage k2, z-component.
#endif

  ParticleReal k3_x = 0.0; ///< Runge-Kutta stage k3, x-component.
  ParticleReal k3_y = 0.0; ///< Runge-Kutta stage k3, y-component.
#if CH_SPACEDIM == 3
  ParticleReal k3_z = 0.0; ///< Runge-Kutta stage k3, z-component.
#endif
};

The velocity columns D_DECL(v_x, v_y, v_z) are the ones required by the solver; the remaining columns hold intermediate Runge-Kutta integration state. Weight and position are container-owned and are therefore not payload members.

Initialization

To initialize the solver, one can use the full constructor

/**
 * @brief Full constructor
 * @param[in] a_amr      Handle to AmrMesh.
 * @param[in] a_compGeom Computational geometry.
 */
TracerParticleSolver(const RefCountedPtr<AmrMesh>& a_amr, const RefCountedPtr<ComputationalGeometry>& a_compGeom);

Getting the particles

To obtain the solver particles, simply call

/**
 * @brief Get all particles.
 * @return m_particles
 */
virtual ParticleContainer<P>&
getParticles();

/**
 * @brief Get all particles. Const version.
 * @return m_particles
 */
virtual const ParticleContainer<P>&
getParticles() const;

This returns the ParticleContainer<P> holding the particles.

Setting \(\mathbf{v}\)

To set the velocity field on the mesh, use

/**
 * @brief Set the tracer particle velocity field.
 * @param[in] a_velocityField Velocity field.
 */
virtual void
setVelocity(const EBAMRCellData& a_velocityField);

This will associate the input velocity a_velocityField with \(\mathbf{v}\).

Interpolating velocities

To compute \(\mathbf{V} = \mathbf{v}\left(\mathbf{X}\right)\) for all particles that reside in the solver, use

/**
 * @brief Interpolate particles velocities.
 */
virtual void
interpolateVelocities();

This will interpolate the velocities to the particle positions using the user-defined interpolation method (see Input options).

If desirable, one can also interpolate a scalar field defined on the mesh onto the particle weight by calling

/**
 * @brief Interpolate a scalar field onto the particle weight.
 * @param[in] a_scalar Scalar field to interpolate from.
 */
virtual void
interpolateWeight(const EBAMRCellData& a_scalar) noexcept;

The interpolation function is set by the user, see Input options. See Particle-mesh for further details.

Deposit particles

To deposit the particles, call

/**
 * @brief Deposit particle weight on mesh.
 * @param[out] a_phi Deposited weight.
 */
virtual void
deposit(EBAMRCellData& a_phi) const noexcept;

This will deposit the particle weights onto the input data holder.

The deposition function is set by the user, see Input options. Complete details regarding how the deposition functions work are available in Particle-mesh.

Input options

Available input options for the tracer particle solver are given in the listing below.

Listing 31 List of configuration options for TracerParticleSolver<P>. All options are run-time configurable.
# ====================================================================================================
# TracerParticleSolver class options
# ====================================================================================================
TracerParticleSolver.verbosity         = -1            ## Solver verbosity level.
TracerParticleSolver.deposition        = cic           ## Deposition method. Must be 'ngp' or 'cic'
TracerParticleSolver.interpolation     = cic           ## Interpolation method. Must be 'ngp' or 'cic'
TracerParticleSolver.deposition_cf     = transition    ## 'interp', 'halo', 'halo_ngp', 'transition'.
TracerParticleSolver.plot_weight       = true          ## Turn on/off plotting of the particle weight.
TracerParticleSolver.plot_velocity     = true          ## Turn on/off plotting of the particle velocities.
TracerParticleSolver.volume_scale      = false         ## If true, depositions yield density * volume instead of just volume