Problem Generators with Particles
Documentation/UserGuide/Problem Generator with Particles
To write a problem generator, one must take the following steps.

Set the total number of particles (
nparticle
), and reallocate the particle array if necessary. For the grid pointer pGrid, one example is given belowif (pGrid>nparticle+2 > pGrid>arrsize) particle_realloc(pGrid, pGrid>nparticle+2);
where
nparticle
+2 is the new array size. Note that the array sizearrsize
must be at least one more thannparticle
. 
Set the parameters needed to calculate the stopping time. When necessary, also set the particle size (
rad
defined inGrain_Property
, see Data Structure) for each particle type.For
tsmode
=3 (fixed stopping time), one simply sets the arraytstop0[]
for the stopping time of each particle type.For
tsmode
=2 (Epstein drag law), one needs to set the global arraygrrhoa[]
for each particle type, so that the stopping time is calculated bygrrhoa/(d*cs)
, whered
andcs
are the gas density and sound speed.For
tsmode
=1 (general drag law), in addition to setting the arraygrrhoa[]
, one further needs to set the global variablealamcoeff
, so that the ratio between particle size and gas mean free path is given byalamcoeff*rad*d
. 
If feedback is included, set the mass for each type of particles (in
grproperty
, see Data Structure). Note that particles in the code are “superparticles”, each representing a swarm of real particles. The particle mass is the inertia of each superparticle, and has the same unit of gas density. 
Specify the initial positions and velocities of all the particles, and assign a unique id for each particle.

The user has the freedom of adding userdefined forces to the particles in the following function
void Userforce_particle(Real3Vect *ft, const Real x1, const Real x2, const Real x3, const Real v1, const Real v2, const Real v3);
The force vector
ft
is to be given as a function of particle position (e.g., in a static gravitational potential). Be warned: if the force vector is a function of velocity in general this will require modifications to the integrator. 
The user can also study the motion of particles in a fake gas velocity field. This can be achieved by letting the gas be in hydrostatic equilibrium (all the time), while creating a fake velocity field using the following function in the problem generator
void gasvshift(Real x1, Real x2, Real x3, Real *u1, Real *u2, Real *u3);
Here the fake gas velocities
u1
,u2
andu3
are to be set as a function of positionx1
,x2
,x3
.