next up previous contents
Next: References to Cite and Up: StarCrash: A Parallel Smoothed Previous: Common Block Variables   Contents


Subroutines

For your convenience, we list the subroutines which comprise our code, with a brief description of each one.

Subroutine ADJUST:
Changes the smoothing length HP of each particle to try to force each particle to maintain NNOPT neighbors, subject the rule ${\tt HPMIN} \le {\tt HP} \le {\tt HPMAX}$. In practice, we change the smoothing length of each particle to the arithmetic average of the old value, and the power-law guess at the new value, given by Eq. 3-9 as:
\begin{displaymath} h_{new}=\frac{1}{2}h_{old}\times \left[1+\left(\frac{N_{opt}}{N_N}\right)^{\frac{1}{3}}\right]. \end{displaymath} (45)

Subroutine ADOTS:
Selects which AV routine to use to calculate the entropy evolution equation, Eq. 3-16
\begin{displaymath} \frac{dA_i}{dt}=\frac{\Gamma-1}{2\rho_i^{\Gamma-1}}\sum_j m_j \Pi_{ij} (\vec{v}_i-\vec{v}_j)\cdot \nabla W_{ij}. \end{displaymath} (46)

The choice of AV subroutine is set by the parameter NAV. Valid choices are: Other choices of NAV will produce an error.

Subroutine ADVANCE:
Evolves all the hydro quantities each timestep. First, the positions are advanced half a timestep, and used to calculate the rate of change of the entropic variable for each particle (ADOTS), before being advanced to the end of the full timestep. Next, velocities are advanced half a timestep, accelerations are calculated (VDOTS), and then the velocities are advanced using the new acceleration a full timestep from their original values.

Subroutine BALADOTS:
Calculates the change in the entropic variable $A_i$ when the Balsara AV prescription is used (Balsara, 1995). See Eqs. 3-16, 3-23, 3-24, and 3-25. See Sec. 3.5.4 for details about parallelization.

Subroutine BALVDOTS:
Calculates the AV acceleration term when the Balsara AV prescription is used (Balsara, 1995). See Eqs. 3-15, 3-23, 3-24, and 3-25. See Sec. 3.5.4 for details about parallelization.

Subroutine BIOUT:
Writes file ``biout.sph'' with binary-related quantities

Subroutine CALCDIVV:
Calculates the SPH form of the divergence of the velocity at each particle position. This is given by Eq. 3-22 as
\begin{displaymath} (\vec{\nabla}\cdot \vec{v})_i={1 \over \rho_i}\sum_j m_j (\vec{v}_j-\vec{v}_i)\cdot\vec{\nabla} W_{ij}. \end{displaymath} (47)

This routine is parallelized like other hydro loops, see Sec. 3.1.3.

Subroutine CHECKPT:
Writes a binary file called ``restart.sph'' every NITCH iterations.

Subroutine CLAADOTS:
Calculates the change in the entropic variable $A_i$ when the Monaghan AV prescription is used (Monaghan, 1989). See Eqs. 3-16, 3-17, 3-18, and 3-19. See Sec. 3.5.4 for details about parallelization.

Subroutine CLAVDOTS:
Calculates the AV acceleration term when the Monaghan AV prescription is used (Monaghan, 1989). See Eqs. 3-15, 3-17, 3-18, and 3-19. See Sec. 3.5.4 for details about parallelization.

Subroutine CMADJ:
Adjusts the system center of mass back to the origin when relaxation is on.

Subroutine crfft:
This C-language routine performs a 3-d Real $\rightarrow$ Complex inverse transform using the FFTW libraries.

Subroutine DENSPLOT:
Generate a compact data file at every iteration for use in animations and other displays. The SPH density is calculated on a grid, and stored in binary form as an integer in the range 0-255, converted to a character.

Subroutine DOQGRAV:
Calculates the gravitational force as well as the radiation reaction force on each particle, when radiation reaction is turned on by setting NGRAVRAD=1. Our scheme for radiation reaction is taken from Blanchet et al. (1990). See Eqs. 3-27, 3-28, 3-29, and 3-30. Note that the speed of light SOL, in code units, must be set for this routine to work.

Subroutine DUMP:
Writes to the logical unit passed an integer input variable (i.e., a previously opened file) the unformatted data which can be used to restart the run. MAKE SURE THE FORMAT AGREES WITH SUBROUTINE INIT, OR YOU WILL NOT BE ABLE TO RESTART RUNS!

Subroutine DUOUT:
Writes a binary file of the form ``outxxx.sph'', where ``xxx'' represents a three-digit integer, every time an additional time interval of ``DTOUT'' passes.

Fuction DW:
Calculates the derivative of the SPH smoothing kernel function W, which itself is given by Eq. 3-2.

Subroutine ENOUT:
Writes file ``energy.sph'' with various energies and other thermodynamic global quantities.

Subroutine GETCURLV:
Calculates the SPH form of the curl of the velocity at each particle position. This is given by Eq. 3-26 as
\begin{displaymath} (\vec{\nabla}\times \vec{v})_i={1 \over \rho_i}\sum_j m_j (\vec{v}_i-\vec{v}_j)\times\vec{\nabla} W_{ij}. \end{displaymath} (48)

This routine is parallelized like other hydro loops, see Sec. 3.1.3.

Subroutine GRAVQUANT:
Calculates n_lower and n_upper, which determine the range of particles whose neighbor lists are stored on a specific processor, since N, the number of particles, is often set at runtime, and can differ from the user-specified value found in sph.input. It also calculates the grid ranges used in the gravity routines, and runs consistency checks on several other numbers, to ensure the runs will not crash. GRAVQUANT MUST BE CALLED AFTER N IS SET, AND BEFORE NENE IS RUN, OR THE CODE WILL CRASH AND DIE IN A HORRIFIC FASHION!

Subroutine GWOUT:
Writes file ``gwdata.sph'' with gravity wave-related quantities

Subroutine INIT:
Initializes the run. If a file named restart.sph exists, it is used to continue a run which had been started previously. If not, it reads in a 3-character code from sph.init and launches the proper subroutine. See Sec. 5.1.2 for the complete list.

After the chosen subroutine finishes, it runs the neighbor finding routine and writes out the parameters for the run.

Subroutine LFSTART:
initializes the leapfrog timing algorithm.

Program MAIN:
Initializes MPI, initializes the hydro routines, opens up the output files, launches the iteration loop, and at the end, closes all files and shuts down MPI.

Subroutine MAINIT:
Runs a single iteration of the SPH code.

Subroutine makeplans:
This C-language subroutine initializes the FFTW library FFT routines, and returns the memory addresses of structures containing plans for calculating the actual Fourier Transforms.

Subroutine NENE:
Computes the nearest-neighbor list for each particle, and saves the information in an array which is distributed over the processors.

We compute the neighbor list using a head-of-cell/linked list method, which assigns all particles to a 3-d array, then loops over the neighboring cell to each particle to construct the neighbor list. All neighbors must lie within two smoothing lengths of the particle.

Subroutine NEWADOTS:
Calculates the change in the entropic variable $A_i$ when the Hernquist and Katz AV prescription is used (Hernquist & Katz, 1989). See Eqs. 3-16, 3-20, and 3-21. See Sec. 3.5.4 for details about parallelization.

Subroutine NEWVDOTS:
Calculates the AV acceleration term when the Hernquist and Katz AV prescription is used (Hernquist & Katz, 1989). See Eqs. 3-15, 3-20, and 3-21. See Sec. 3.5.4 for details about parallelization.

Subroutine OPTHP2:
Calculates the approximately optimal particle smoothing length for runs in which equal-mass particles are laid down Monte Carlo style. The equation is given by
\begin{displaymath} h_i=0.9 R\left(\frac{3}{32\pi}\frac{N_{opt} M}{N \rho_i}\right)^{\frac{1}{2}}. \end{displaymath} (49)

where $M$ and $R$ are the star's mass and radius, $\rho_i$ is the ideal density at that particle's position, and $N$ and $N_{opt}$ are the total number of particles in the star and the optimal number of neighbors, respectively.

Subroutine OPTHP3:
Calculates the approximately optimal particle smoothing length for the secondary star for unequal-mass stars with equal-mass particles (subroutine SETUP2QM), using Eq. 7-5, with $R$ and $M$ appropriate for the secondary.

Subroutine OUTPUT:
Calls out all the subroutines which write data to files, having first synchronized the velocities to the same time values as the positions according to the leapfrog algorithm (see Sec. 3.1.5).

Subroutine POLY:
Calculates the value of the entropic variable $A_i$ and the density profile for a spherical polytrope with specified index, mass and radius.

This is one of the very few subroutines which has values passed to it during the subroutine call. They are, in order, the adiabatic index $n\equiv 1/(\Gamma-1)$, the stellar mass, the stellar radius, and the size of the vector containing the density profile. It returns the density profile as a vector, and the vale of the entropic variable for the star. Note that this routine does NOT include the spha.h header file, and thus has no access to common block variables. Note also that it calls subroutines DERIVS, RKTAB, RK4, which are the standard Numerical Recipes routine for doing a fourth order Runge-Kutta integration (Press et al., 1992).

Subroutine QGRAV:
Calculates the gravitational force on each particle by means of an FFT convolution algorithm. We use the MPI-parallelized routines of the FFTW software package, available free at http://www.fftw.org.

The Green's function for the transform requires a grid twice as large as the grid which holds the data, a method known as zero-padding the data. Essentially, we need to handle independently the effect of particles located at both higher and lower x, y, and z, and one can check that the range of possible values of $\vec{r}_i-\vec{r}_j$ is twice as large in every direction as the data itself.

Densities are uploaded to the grid by a cloud-in-cell algorithm, i.e. densities are assigned with weights to the 8 grid points around any particle. We then calculate the FFT of the data, convolve with the $1/r$ Green's function, and inverse FFT. As this is a parallel FFT, with 3-d arrays distributed along the z-direction, we need to take additional steps to finite difference in the z-direction when we are near the edge of each processor's data segment.

Subroutine QGRAV2:
Calculates the gravitational force on each particle by means of an FFT convolution algorithm, as well as all second derivatives of the gravitational potential, which are needed for calculating the radiation reaction force. Other than that, this routine is an exact replica of subroutine QGRAV

Subroutine QGRAVRAD:
Calculates the radiation reaction potential given by Eq. 3-30, in the notation of the BDS radiation reaction formalism (Blanchet et al., 1990).

Function RAN1:
Random number generator, returning a value $0 < {\tt RAN1} < 1$. Use IDUM=-2391 as a seed. Note that the generator is fully deterministic, and thus your runs are repeatable.

Subroutine rcfft:
This C-language routine performs a 3-d Complex $\rightarrow$ Real forward transform using the FFTW libraries.

Subroutine RELAX:
Adds the relaxation drag terms to the equations of motion. For single stars, the drag force is given by
\begin{displaymath} \frac{d\vec{v}}{dt}=\vec{a}-\frac{\vec{v}}{t_{relax}}, \end{displaymath} (50)

where the user sets $t_{relax}$. For binary stars, the calculation is performed in the corotating frame, and the relaxation force takes the form
\begin{displaymath} \frac{d\vec{v}}{dt}=\vec{a}+\Omega^2 \vec{r}-\frac{\vec{v}}{t_{relax}}, \end{displaymath} (51)

where $\Omega$ is calculated from the inward force on each star.

Subroutine RHOS:
Calculates the SPH expression for the particle density, given by
\begin{displaymath} \rho_i=\sum_j m_j W_{ij}, \end{displaymath} (52)

where the sum is taken over all the neighbors of a particle. We use a ``gather-scatter'' algorithm, i.e., for each particle-neighbor pair, we count half of the density contribution to the particle, and half to the neighbor. This routine uses MPI to distribute the computation.

Subroutine SETUP1EM:
Constructs a single star with equal-mass particles, laid down Monte-Carlo style weighted by the density at a given radius according to a polytropic fit. Called by using the 3-letter code ``1em'' in sph.init.

Subroutine SETUP1ES:
Constructs a single star with equally-spaced particles of varying mass, which yield a density profile very near to a polytropic fit. Called by using the 3-letter code ``1es'' in sph.init.

Subroutine SETUP2CM:
Constructs a corotating binary containing equal-mass stars, using equal-mass particles, laid down Monte-Carlo style weighted by the density at a given radius according to a polytropic fit. Called by using the 3-letter code ``2cm'' in sph.init.

Subroutine SETUP2CR:
Constructs a corotating binary containing equal-mass stars, using the output file from a single star run named ``pos.sph''. Called by using the 3-letter code ``2cr'' in sph.init.

Subroutine SETUP2CS:
Constructs a corotating binary containing equal-mass stars, using equally-spaced particles of varying mass, which yield a density profile very near to a polytropic fit. Called by using the 3-letter code ``2cs'' in sph.init.

Subroutine SETUP2I1:
Constructs an irrotational binary containing equal-mass stars, using the output file from a single star run named ``pos.sph'', and axis ratios specified by a field called ``axes1.sph''. Called by using the 3-letter code ``2i1'' in sph.init.

Subroutine SETUP2IQ:
Constructs an irrotational binary containing unequal-mass stars, using the output filed from single star run named ``pos.sph'' and ``pos2.sph'', and axis ratios specified by a field called ``axes1.sph'' and ``axes2.sph''. Called by using the 3-letter code ``2iq'' in sph.init.

Subroutine SETUP2QM:
Constructs a corotating binary containing unequal-mass stars, using equal-mass particles, laid down Monte-Carlo style weighted by the density at a given radius according to a polytropic fit. Called by using the 3-letter code ``2qm'' in sph.init.

Subroutine SETUP2QR:
Constructs a corotating binary containing unequal-mass stars, using the output files from single star runs named ``pos.sph'' and ``pos2.sph''. Called by using the 3-letter code ``2qr'' in sph.init.

Subroutine SETUP2QS:
Constructs a corotating binary containing unequal-mass stars, using equally-spaced particles of varying mass, which yield a density profile very near to a polytropic fit. Called by using the 3-letter code ``2qs'' in sph.init.

Subroutine SETUPHYP1:
Constructs a binary on a hyperbolic orbit containing equal-mass stars, using the output file from a single star run named ``pos.sph''. Called by using the 3-letter code ``hy1'' in sph.init.

Subroutine SETUPHYPQ:
Constructs a binary on a hyperbolic orbit containing unequal-mass stars, using the output files from single star runs named ``pos.sph'' and ``pos2.sph''. Called by using the 3-letter code ``hyq'' in sph.init.

Subroutine TABULINIT:
calculates the SPH smoothing kernel W, used to calculate all hydrodynamic quantities, and its derivative. The kernel function we use is given by
\begin{displaymath} W(r,h)={1\over\pi h^3}\cases{1-{3\over2}\left({r\over h}\rig... ...\right]^3,& $1\le{r\over h}<2$,\cr 0,& ${r\over h}\ge2$.\cr} \end{displaymath} (53)

Subroutine TSTEP:
Calculates the time step used in the code.

Subroutine VDOTS:
Selects which AV routine to use to calculate the all forces on the particle, including the AV acceleration, Eq. 3-15. The choice of AV subroutine is set by the parameter NAV. Valid choices are: Other choices of NAV will produce an error.

Function W:
Calculates the SPH smoothing kernel function W, which is given by Eq. 3-2.

next up previous contents
Next: References to Cite and Up: StarCrash: A Parallel Smoothed Previous: Common Block Variables   Contents
Joshua Faber 2003-06-28