A Smoothed Particle Hydrodynamics code for smashing stars (and planets)

View project on GitHub

Please contact Jamie Lombardi ( with any questions.


Get the code from the repository at

The main directory contains at least three subdirectories:

example_input: contains example input files for various cases
misc: miscellaneous files such as assorted makefiles and the StarCrash manual
parallel_bleeding_edge: latest version of code (this is probably what you want)

For runs on a gpu-enabled cluster, the directory you care about is parallel_bleeding_edge. Let's try to compile the code in there:

cd parallel_bleeding_edge/src/SPHgrav_lib

The SPHgrav_lib subdirectory contains code written by Evghenii Gaburov (and somewhat modified by Jamie Lombardi and Sam Knarr) for calculating softened gravitational forces and potentials on NVIDIA GPUs. If you are running on something other than keeneland, you will need to customize the Makefile for your system: any path that contains "cuda" or "CUDA" probably needs to be changed. To make the library, simply issue the following command:


(The final thing that make does is ./ which is what will make the library.) If it doesn't work, and you want to have a fresh start, first use make clean. If all goes well, the library by the name libGPUsph_gatherscatter.a will be generated.

If you type the command below, you will be moved up a level in the directory structure. The SPH source code is in the ".." directory. This is a parallel code. Execute the following command to move into the src subdirectory.

cd ..

Now, make a copy of the makefile.

cp makefile.keeneland makefile

To compile, type the command below.


You will likely need to customize the makefile for your computer. For example, you may need to change the compiler or update paths to point to the appropriate locations. The executable name will end with _sph and will depend on the name of the directory that src appears within. The executable will automatically be moved up a level in the directory structure (that is, to the ".." directory). You probably want to remove the executable for now so that you don't later copy it to other places where you don't want it:

cd ..
rm *_sph

Although this Starsmasher code is not well documented, there is documentation for Starcrash, which is a previous version of this code (see Although Starcrash is not the same as Starsmasher, the variable names, input files, output files, and parallelization strategy are very similar. The Starcrash documentation is available as the file usersmanual.pdf in the misc subdirectory and may be somewhat helpful to peruse. You will be most interested in pages 31 - 36, which talks about input and output. Remember, however, that this Starsmasher SPH code is different in many ways. The algorithms of Starsmasher are mostly described in The AV scheme is described in Ponce et al. (2011).


There are two steps to modelling a collision or fly-by of two stars. First, you have to model the individual star: a relaxation run. Second, you simulate the interaction: a dynamical calculation.

Step 1: Relaxation runs

Step 1(a): Relaxing a n=1.5 polytrope

Let's have you do your own simulation for practice. We'll use a small number of particles N so that the code runs very quickly. Let's try an equation of state governed by the adiabatic index gamma=5/3, and a structure specified by the polytropic index n=1.5.

First, let's get the directories ready. We'll be using parallel_bleeding_edge as a template directory that can be duplicated as needed. So assuming you are still in the directory parallel_bleeding_edge, do a
cd ..
(or maybe cd ~/starsmasher depending on where your installation is).
cp -r parallel_bleeding_edge GAM1.667_n1.5
cd GAM1.667_n1.5
cp ../example_input/relaxation_preMS/sph.init .
cp ../example_input/relaxation_preMS/sph.input .

The file sph.init controls what type of calculation we will do.
Confirm that the three letter code word inside sph.init is '1es' for this polytrope relaxation (and update the text if necessary):
emacs -nw sph.init

You should look at sph.input:
emacs -nw sph.input
[or use some other editor]

Change line 4 of sph.input to be N=2000 and save the file. There are several other variables that could be changed, but the defaults should be fine for now. See the Starcrash documentation for an explanation of what some of the other variables are.

Now let's compile the source code:
cd src
cat starsmasher.h

Confirm that the last line sets the adiabatic index GAM=5.d0/3.d0.
emacs -nw initialize_polyes.f

Confirm that the lines right after the declaration section will set the polytropic constant npoly to be 1.5.
The radius and mass are set on near here. These values are in physical units of runit and munit, which are set at the end of sph.input and are typically the radius and mass of the sun in cgs units. Save initialize_polyes.f when you are done and then try to compile:

Once everything goes well with the compilation, the last line you should see is "mv GAM1.667_n1.5_sph ..".
Get positioned to run the code:
cd ..

And get a good batch script if you don't have one:
cp ../sph.pbs.kfs sph.pbs

The default sph.pbs file asks for only 1 minute, so you'll need to change that for a real run. To run the code,
qsub -q debug sph.pbs ### or, for long jobs, "qsub sph.pbs"

This creates, among other things, a file log0.sph that collects ascii tex about the run.

Step 1(b): Relaxing a TWIN model

Get the directories ready:
cd ~/starsmasher ### or cd to wherever your main installation is
cp -r parallel_bleeding_edge TWIN
cp ../example_input/relaxation_TWIN_model/* .

To relax a TWIN model instead, the three letter code word in sph.init should be 'erg'. The SPH code will look for a model named 'eg.last1.muse_s2mm'.
You can create a symbolic link from an appropriate output file of TWIN if you like
(e.g., "cp -s m_8_t_5.7.last1.muse_s2mm eg.last1.muse_s2mm").
The file m_8_t_5.7.last1.muse_s2mm is the default in example_input/relaxation_TWIN_model. It represents an 8M_sun MS star at t=5.7Myr.

The sph.input file, also available within example_input/relaxation_TWIN_model/, should be similar to the sph.input file used above for the polytrope (probably just change N and nothing else unless necessary). The file on the svn site does have TRELOFF to 200 in this case, to make sure there was enough time for oscillations to die off, and TF is 300 so the star could be monitored for 100 time units after the relaxation was turned off.

Run the code like before:
cd src
cd ..
cp ../sph.pbs.kfs ./sph.pbs
qsub -q debug sph.pbs

Step 2: Dynamical calculation

If the model from the relaxation run at t=TRELOFF looks good, then you can use that to start a dynamical calculation. Start by get the directories ready:
cd ~/starsmasher ### or cd to wherever your main installation is
cp -r parallel_bleeding_edge/ collision_rp3.9_a119
cd collision_rp3.9_a119
cp ../example_input/collision/* .
emacs -nw sph.init

Make sure line 2 of sph.init says INAME='hyp'.

The new sph.input file is set up for doing dynamical calculations. Take a look at it:
emacs -nw sph.input

While you're in the file, check that DTOUT=1 on line 3, which means that an out*.sph file will be dumped every 1 time unit. Also, on line 16 or so, make BIMPACT=3.9d0, which means the periastron separation of the initial orbit is 3.9 stellar radii.

In sph.input you can use any of the following pairs of data to initialize a collision in hyperbolic.f: (e0,vinf2), (e0,bimpact), (semimajoraxis,bimpact), (semimajoraxis,e0), (bimpact,vinf2). Here, e0 is the eccentricity of the initial orbit, and bimpact is the periastron separation (not really the impact parameter). For example, sph.input could contain the following.


and e0 and vinf2 could be unspecified. Or you could specify BIMPACT and e0, and the same code works, etc. The code in initialize_hyperbolic.f figures out what it needs to solve for.

Whether the encounter is hyperbolic, parabolic, or elliptical is controlled by the velocity at infinity squared, vinf2, or the semi-major axis a. The value of vinf2 is negative for elliptical encounters, zero for parabolic encounters, and positive for hyperbolic encounters. Note that the orbital energy GMmu/(-2a) just equals muvinf2/2.
Therefore, there is a simple relation between the semi-major axis a and vinf2, namely. a=-GM/vinf2, where M is the total mass of the two stars. For example, if you are colliding a 0.3 Msun and 8 Msun star and if vinf2=-0.07, then the semi-major axis a=(0.3+8)/0.07=118.57, where I assume the units chosen are the usual G=Msun=Rsun=1.

Note that ngravprocs=-6 and ppn=16... this means there are abs(ngravprocs)=6 GPU units per 16 CPU cores, as on one node of the supercomputer forge. You should change these values to be appropriate for your machine. For example, on Keeneland KFS, there are 3 gpus and 12 cores per node.

Let's get the star models:
cp ../TWIN//out0200.sph ./sph.start1u
cp ../GAM1.667_n1.5/
/out0100.sph ./sph.start2u

Where we assume the 8 solar mass ZAMS star was relaxed in the directory ../TWIN/.

Now let's compile the code:
cd src

To run the code, get to the same directory as the executable:
cd ..

cp ../sph.pbs.kfs ./sph.pbs
qsub -q debug sph.pbs


Every few iterations, a restartrad.sph file is dumped (overwriting any previously existing restartrad.sph file). If a restartrad.sph file exists in the directory when a new run is launched, then the code will automatically use that file to initiate the calculation. This is useful for restarting in the middle of a long simulation. You can also restart from any out*.sph file simply by renaming it to restartrad.sph.


I recommend Price's SPLASH software. It will need to be customized to read in our output file type.