Please contact Jamie Lombardi (jamie.lombardi@allegheny.edu) with any questions.
OBTAINING AND COMPILING
Get the code from the repository at https://github.com/jalombar/starsmasher/.
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:
make
(The final thing that make does is ./makelib.sh 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.
make
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 http://ciera.northwestern.edu/StarCrash/). 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 http://adsabs.harvard.edu/abs/2010MNRAS.402..105G. The AV scheme is described in Ponce et al. (2011).
EXAMPLE RUNS
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:
make
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
cd 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
make
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/sph.in* .
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.
BIMPACT=4.d0,
semimajoraxis=118.57d0,
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
make
To run the code, get to the same directory as the executable:
cd ..
Then,
cp ../sph.pbs.kfs ./sph.pbs
qsub -q debug sph.pbs
RESTARTING RUNS
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.
DATA VISUALIZATION
I recommend Price's SPLASH software. It will need to be customized to read in our output file type.