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 below
if (pGrid->nparticle+2 > pGrid->arrsize) particle_realloc(pGrid, pGrid->nparticle+2);
nparticle+2 is the new array size. Note that the array size
arrsizemust be at least one more than
Set the parameters needed to calculate the stopping time. When necessary, also set the particle size (
Grain_Property, see Data Structure) for each particle type.
tsmode=3 (fixed stopping time), one simply sets the array
tstop0for the stopping time of each particle type.
tsmode=2 (Epstein drag law), one needs to set the global array
grrhoafor each particle type, so that the stopping time is calculated by
csare the gas density and sound speed.
tsmode=1 (general drag law), in addition to setting the array
grrhoa, one further needs to set the global variable
alamcoeff, so that the ratio between particle size and gas mean free path is given by
If feedback is included, set the mass for each type of particles (in
grproperty, see Data Structure). Note that particles in the code are “super-particles”, each representing a swarm of real particles. The particle mass is the inertia of each super-particle, 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 user-defined 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
ftis 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
u3are to be set as a function of position