20.3 Using the Particles Unit

The Particles unit encompasses nearly all aspects of Lagrangian particles. The exceptions are input/output the movement of related data structures between different blocks as the particles move from one block to another, and mapping the particle attributes to and from the grid.

Beginning with release of version 4 it is possible to add particles to a simulation during evolution, a new function Particles_addNew has been added to the unit's API for this purpose. It has been possible to include multiple different types of particles in the same simulation since release FLASH3.1. Particle types must be specified in the Config file of the Simulations unit setup directory for the application, and the syntax is explained in Sec:FlashHparticles. At configuration time, the setup script parses the PARTICLETYPE specifications in the Config files, and generates an F90 file Particles_specifyMethods.F90 that populates a data structure gr_ptTypeInfo. This data structure contains information about the method of initialization and interpolation methods for mapping the particle attributes to and from the grid for each included particle type. Different time integration schemes are applied to active and passive particles. However, in one simulation, all active particles are integrated using the same scheme, regardless of how many active types exists. Similarly, only one passive integration scheme is used. The added complexity of multiple particle types allows different methods to be used for initialization of particles positions and their mapping to and from the grid quantities. Because several different implementations of each type of functionality can co-exist in one simulation, there are no defaults in the Particles unit Config files. These various functionalities are organized into different subunits; a brief description of each subunit is included below and further expanded in subsections in this chapter.

• The ParticlesInitialization subunit distributes a given set of particles through the spatial domain at the simulation startup. Some type of spatial initialization is always required; the functionality is provided by Particles_initPositions. The users of active particles typically have their own custom initialization. The following two implementations of initialization techniques are included in the FLASH4 distribution (they are more likely to used with the passive tracer particles):

  Lattice

  distributes particles regularly along the axes directions throughout a subsection of the physical grid.

  WithDensity

  distributes particles randomly, with particle density being proportional to the grid gas density.

  Users have two options for implementing custom initialization methods. The two files involved in the process are: Particles_initPositions and pt_initPositions. The former does some housekeeping such as allowing for inclusion of one of the available methods along with the user specified one, and assigning tags at the end. A user wishing to add one custom method with no constraints on tags etc is advised to implement a custom version of the latter. This approach allows the user to focus the implementation on the placement of particles only. Users desirous of refining the grid based on particles count during initialiation should see the setup PoisParticles for an example implementation of the Particles_initPositions routine. If more than one implementation of pt_initPositions is desired in the same simulation then it is necessary to implement each one separately with different names (as we do for tracer particles: pt_initPositionsLattice and pt_initPositionsWithDensity) in their simulation setup directory. In addition, a modified copy of Particles_initPostions, which calls these new routines in the loop over types, must also be placed in the same directory.

• The ParticlesMain subunit contains the various time-integration options for both active and passive particles. A detailed overview of the different schemes is given in Sec:Particles Integration.

• The ParticlesMapping subunit controls the mapping of particle properties to and from the grid. FLASH currently supplies the following mapping schemes:

  Cloud-in-cell

  (ParticlesMapping/meshWeighting/CIC), which weights values at nearby grid cells; and

  Quadratic

  (ParticlesMapping/Quadratic), which performs quadratic interpolation.

  Some form of mapping must always be included when running a simulation with particles. As mentioned in Sec:Particles Mapping the quadratic mapping scheme is only available to map from the grid quantities to the corresponding particle attributes. Since active particles require the same mapping scheme to be used in mapping to and from the mesh, they cannot use the quadratic mapping scheme as currently implemented in FLASH4. The CIC scheme may be used by both the active and passive particles.

  For active particles, we use the mapping routines to assign particles' mass to the particle density grid-based solution variable (PDEN_VAR). This mapping is the initial step in the particle-mesh (PM) technique for evaluating the long range gravitational force between all particles. Here, we use the particle mapping routine Particles_mapToMeshOneBlk to “smear” the particles' attribute over the cells of a temporary array. The temporary array is an input argument which is passed from the grid mapping routine Grid_mapParticlesToMesh. This encapsulation means that the particle mapping routine is independent of the current state of the grid, and is not tied to a particular Grid implementation. For details about the task of mapping the temporary array values to the cells of the appropriate block(s), please see Sec:GridParticlesMapToMesh. New schemes can be created that “smear” the particle across many more cells to give a more accurate PDEN_VAR distribution, and thus a higher quality force approximation between particles. Any new scheme should implement a customized version of the pt_assignWeights routine, so that it can be used by the Particles_mapToMeshOneBlk routine during the map.

• The ParticlesForces subunit implements the long and short range forces described in Equation (20.3) in the following directories:
  • longRange collects different long-range force laws (requiring elliptic solvers or the like and dependent upon all other particles);
  • shortRange collects different short-range force laws (directly summed or dependent upon nearest neighbors only).

  Currently, only one long-range force law (gravitation) with one force method (particle-mesh) is included with FLASH. Long-range force laws are contained in the Particles/ParticlesForces/longRange, which requires that the Gravity unit be included in the code. In the current release, no shortRange implementation of ParticlesForces is supplied with FLASH. However, note that the sink particle implementation described below in Sec:Particles Sink includes directly computed particle-particle forces.

After particles are moved during time integration or by forces, they may end up on other blocks within or without the current processor. The redistribution of particles among processors is handled by the GridParticles subunit, as the algorithms required vary considerably between the grid implementations. The boundary conditions are also implemented by the GridParticles unit. See Sec:GridParticles for more details of these redistribution algorithms. The user should include the option -with-unit=Grid/GridParticles on the setup line, or REQUIRES Grid/GridParticles in the Config file.

In addition, the input-output routines for the Particles unit are contained in a subunit IOParticles. Particles are written to the main checkpoint files. If the user desires, a separate output file can be created which contains only the particle information. See Sec:Particles IO below as well as Sec:IOParticles for more details. The user should include the option -with-unit=IO/IOParticles on the setup line, or REQUIRES IO/IOParticles in the Config file.

In FLASH4, the initial particle positions can be used to construct an appropriately refined grid, i.e. more refined in places where there is a clustering of particles. To use this feature the flash.par file must include: refine_on_particle_count=.true. and max_particles_per_blk=[some value]. Please be aware that FLASH will abort if the criterion is too demanding. To overcome the abort, specify a less demanding criterion, or increase the value of lrefine_max.

20.3.1 Particles Runtime Parameters

There are several general runtime parameters applicable to the Particles unit, which affect every implementation. The variable useParticles obviously must be set equal to .true. to utilize the Particles unit. If the variable initializeParticleAtRestart is .true., then particle data is initilized from scratch instead of a checkpoint file when restarting a simulation. The time stepping is controlled with pt_dtFactor; a value less than one ensures that particles will not step farther than one entire cell in any given time interval. The Lattice initialization routines have additional parameters. The number of evenly spaced particles is controlled in each direction by pt_numX and similar variables in $Y$ and $Z$. The physical range of initialization is controlled by pt_initialXMin and the like. Finally, note that the output of particle properties to special particle files is controlled by runtime parameters found in the IO unit. See Sec:Particle files for more details.

20.3.2 Particle Attributes

By default, particles are defined to have eight real properties or attributes: 3 positions in x,y,z; 3 velocities in x,y,z; the current block identification number; and a tag which uniquely identifies the particle. Additional properties can be defined for each particle. For example, active particles usually have the additional properties of mass and acceleration (needed for the integration routines, see Table Table 20.1). Depending upon the simulation, the user can define particle properties in a manner similar to that used for mesh-based solution variables. To define a particle attribute, add to a Config file a line of the form

PARTICLEPROP property-name

For attributes that are meant to merely sample and record the state of certain mesh variables along trajectories, FLASH can automatically invoke interpolation (or, in general, some map method) to generate attribute values from the appropriate grid quantities. (For passive tracer particles, these are typically the only attributes beyond the default set of eight mentioned above.) The routine Particles_updateAttributes is invoked by FLASH at appropriate times to effect this mapping, namely before writing particle data to checkpoint and particle plot files. To direct the default implementation of Particles_updateAttributes to act as desired for tracer attributes, the user must define the association of the particle attribute with the appropriate mesh variable by including the following line in the Config file:

PARTICLEMAP TO property-name FROM VARIABLE variable-name

These particle attributes are carried along in the simulation and output in the checkpoint files. At runtime, the user can specify the attributes to output through runtime parameters
particle_attribute_1, particle_attribute_2, etc. These specified attributes are collected in an array by the Particles_init routine. This array in turn is used by Particles_updateAttributes to calculate the values of the specified attributes from the corresponding mesh quantities before they are output.

20.3.3 Particle I/O

Particle data are written to and read from checkpoint files by the I/O modules (Sec:FLASH output formats). For more information on the format of particle data written to output files, see Sec:HDF5 and Sec:PnetCDF IO.

Particle data can also be written out to the flash.dat file. The user should include a local copy of IO_writeIntegralQuantities in their Simulation directory. The Orbit test problem supplies an example IO_writeIntegralQuantities routine that is useful for writing individual particle trajectories to disk at every timestep.

There is also a utility routine Particles_dump which can be used to dump particle output to a plain text file. An example of usage can be found in Particles_unitTest. Output from this routine can be read using the fidlr routine particles_dump.pro.

20.3.4 Unit Tests

The unit tests provided for Particles exercise the Particles_advance methods for tracer particles. Tests under Simulation/SimulationMain/unitTest/ParticlesAdvance can be used to examine and compare convergence behavior of various time integration schemes. The tests compare numerical and analytic solutions for a problem (with a given velocity field) where analytic solutions can be computed.

Currently only one ParticlesAdvance test is provided. It is designed to be easily modified by replacing a few source files that contain implementations of the equation and the analytic solution. The use the test, configure it with a command like  

./setup -auto -1d unitTest/ParticlesAdvance/HomologousPassive    		-unit=Particles/ParticlesMain/passive/EstiMidpoint2

and replace EstiMidpoint2 with one of the other available methods (or omit the option to get the default method), see Sec:Passive Partices Integration. Add other options as desired.

For unitTest/ParticlesAdvance/HomologousPassive,  

./setup -auto -1d unitTest/ParticlesAdvance/HomologousPassive +ug -nxb=80

is recommended to get started.

When varying the test, the following runtime parameters defined for
Simulation/SimulationMain/unitTest/ParticlesAdvance will probably need to be adjusted:

PARAMETER sim_schemeOrder INTEGER 2 --

The order of the integration scheme. This should probably always be either 1 or 2.

PARAMETER sim_maxTolCoeff0 REAL 1.0e-8 --

Zero-th order error coefficient $C_0$, used for convergence criterion if sim_schemeOrder$=0$.

PARAMETER sim_maxTolCoeff1 REAL 0.0001 --

First order error coefficient $C_1$, used for convergence criterion if sim_schemeOrder$=1$.

PARAMETER sim_maxTolCoeff2 REAL 0.01 --

Second order error coefficient $C_2$, used for convergence criterion if sim_schemeOrder$=2$.

A test for order $k$ is considered successful if the following criterion is satisfied:

$${\tt maxError} \leq C_k \times {\tt maxActualDt}^k \;,$$

where maxError is the maximum absolute error between numerical and analytic solution for any particle that was encountered during a simulation run, and maxActualDt is the maximum time step $\Delta t$ used in the run. The appropriate runtime parameters of various units, in particular Driver, Particles, and Grid, should be used to control the desired simulation run. In particular, it is recommended to vary dtmax by several orders of magnitude (over a range where it directly determines maxActualDt) for a given test in order to examine convergence behavior.