Running a simulation

Two files control a run. sph.init chooses the kind of calculation, and sph.input sets everything else. Both must be present.

The minimum sph.init is three lines, for example:

&INITT
INAME='erg'
&END

INAME is a three-letter code naming the setup routine. Two of them cover most work:

erg

Build a star from a stellar-evolution profile, a MESA model for instance, and relax it into an SPH star. This is where nearly every project starts, and it is the subject of Relaxing a star.

hyp

Put two already-relaxed bodies on a Keplerian orbit and let them meet: a collision or fly-by between two stars, or between a star and a compact object. See A collision.

The usual sequence is therefore two runs: an erg run per star to make it, then one hyp run to collide them.

The other thirteen initialization scripts

Code	Sets up
1es	a single polytrope, from starmass and starradius
1mc	a single polytrope with a compact core
2cr	a corotating binary
bps	a binary plus a single star
bph	a binary plus a black hole
hbs	a binary encountering a single star
meq	a star whose particle-mass scheme varies with position
res	a rescaling of an existing model
tri	a triple system
bhe	an encounter with a supermassive black hole
grs	a general-relativistic setup
rin	a restart from an existing dump
txt	particles laid out from an ASCII image

sph.init gives the same list with the routine each one calls.

Each code selects one initialization script, and they live in files named initialize_*.f in parallel_bleeding_edge/src: initialize_parent.f for erg, initialize_hyperbolic.f for hyp, and so on.

Which script you choose also decides which other files the run needs: a stellar-evolution profile for erg, a relaxed star or two for hyp. See Input files for what those are and where they come from.

The important thing about these routines is how little they do. A setup routine runs once, before the first step, and its only job is to decide where the particles start: their positions, masses, velocities and internal energies. It then hands that state to the integrator and takes no further part. Everything after that is the same code regardless of which INAME you chose: the hydrodynamics, the gravity, the timestep control, and what gets written to disk.

Two consequences are worth keeping in mind. A setting that controls the evolution means the same thing for every INAME, so advice about cn1 or nnopt carries across unchanged. A setting that describes the starting state is read only by the routine that needs it, so e0 does nothing in a relaxation and equalmass does nothing in a collision, and neither will warn you.

See sph.input for every setting and its default.

Choosing the number of ranks

mpirun -np 8 ./my_run_gpu_sph

Worth testing on your own machine, but sixteen ranks is often a sensible ceiling. Beyond that, communication costs eat into the gains, and leaving cores free lets you run several simulations at once.

ngravprocs sets how many ranks compute gravity. It cannot exceed the number of GPUs present: ask for more and it is clamped to what you have. Leave it at 0 to have it detected automatically.

GPU or CPU

make builds the GPU version and make cpu the CPU one. The GPU takes over only the gravity: neighbour finding and the SPH sums run on the CPU either way. That is why the two are indistinguishable on a small model and diverge sharply on a large one. Gravity is the part that grows fastest with the number of particles, and it is the only part the card touches.

The card computes gravity in single precision. Positions are sent to it relative to the gas centroid weighted by \(1/h^2\), with the subtraction done in double precision, so a star far from the origin, as near a supermassive black hole, keeps its close separations. That point sits on a dense core, and midway between two identical stars. Gas in two dense places far apart, such as two stars many thousands of smoothing lengths from each other, still loses precision on the card, and is better run with make cpu.

What was measured

Relaxing a one-solar-mass polytrope, INAME='erg' aside, with:

&input
n=N, nnopt=22, tf=1.0, dtout=1000.0, nrelax=1, treloff=0,
starmass=1.0, starradius=1.0,
&end

on four MPI ranks, on one machine: an RTX 5080 with the CPU build using all four ranks for gravity. Both builds were checked to produce the same answer, with identical W, U and total energy to every digit written, so the two columns are the same work, not merely similar work.

\(N\)	Iterations	GPU, per step	CPU, per step	CPU/GPU	CPU, whole run
1 000	23	82 ms	75 ms	0.9	1.7 s
4 000	30	68 ms	69 ms	1.0	2.1 s
16 000	43	68 ms	171 ms	2.5	7.4 s
64 000	64	129 ms	1817 ms	14	116 s
128 000	80	225 ms	7104 ms	32	568 s

The per-step columns are the comparison worth making. The iteration count is not a constant of the experiment: every run above stops at the same tf, but a bigger model takes shorter timesteps and so needs more of them, 23 steps at \(N=1000\) against 80 at \(N=128\,000\). A table of total wall times alone would mix that in with the cost of a step and show neither clearly.

What this means in practice

Below about \(N=4000\) it makes no difference. The GPU is fractionally slower, because there is not enough work to cover the cost of moving it to the card. By \(N=16\,000\) the CPU is taking two and a half times as long, and by \(N=128\,000\) it is thirty times.

So a relaxation of a few thousand particles is perfectly comfortable without a card, and the tutorials in this documentation are sized to run either way. A production model of \(10^5\) particles is not: the same job that occupies a GPU for a few minutes will occupy the CPU build for hours.

These figures are one machine and one kind of run. The protocol is written above so you can repeat it on yours rather than trusting the numbers, and it is worth doing before committing to a long job.

Warning

If you have built both versions, ./*_sph matches two files and the shell hands you the alphabetically first, which is the CPU build. It runs, just more slowly, and nothing says so. The give-away is the first lines of output: the GPU build reports SPHgrav found 1 CUDA devices and is running on ... with gpu 0, and the CPU build reports neither.

Block timesteps

By default every particle advances with one shared timestep, the smallest that any particle needs. In an encounter with a black hole, a handful of gas particles close to the black hole can need steps thousands of times shorter than the rest of the star, and the whole star then pays for them.

Setting nblock=1 gives each particle its own step instead, a power-of-2 fraction of the largest step dtmaxblk (dtout when dtmaxblk is zero or less). The smallest allowed step is dtmaxblk/2**nbinmax. Each particle is integrated with its own kick-drift-kick leapfrog. Particles that are not active in a substep are drifted to the current time, and their smoothing lengths and densities are predicted rather than solved. For the force calculation their positions also include the half acceleration term, so that they sit on the parabola of their step rather than its chord. Without it, a star falling in a strong external field (a black hole’s) sees its inactive particles lag behind by up to a dt**2/8, and momentum and energy drift at a rate the timestep limits do not control. Two safeguards keep neighbours consistent:

• nblimit (default 2) keeps every particle’s step within a factor 2**nblimit of its neighbours’ steps, and a particle whose neighbour suddenly needs a much shorter step is woken up early (Saitoh & Makino 2009).

• Any inactive particle with an active particle inside its kernel is refreshed: its smoothing length, density and related sums are re-solved at the current positions, so that the active particles see up-to-date values. The hydrodynamic forces then come from the same routine as with shared steps, using the neighbour lists of the active and refreshed particles.

When only a few particles are active, the work of a substep is arranged so that its cost depends on those particles rather than on the whole star. Their neighbours are found in a search tree that is kept across substeps and rebuilt only when particles have moved too far for it to be trusted, with a search radius enlarged to guarantee that every current neighbour is found. The active and refreshed particles are dealt out evenly across MPI ranks, whatever their particle numbers, and the black hole to gas timestep check is shared the same way.

The black hole to gas timestep limits from cn5, cn6 and cn7 are applied to both members of each pair. Energies are written only when all particles are in step, every dtmaxblk, and so is restartrad.sph: at the first such time at least 1000 substeps after the previous one. Output files and restarts work the same way in both modes, so a run can be restarted with or without block timesteps.

Block timesteps need a dynamical run (nrelax=0) with nintvar=2, ncooling=0, nkernel=2 and nselfgravity=1, and the code stops with a message otherwise. They have been tested with the CPU gravity build only.

Block timesteps work with automatic jumps ahead (rjump or rjumpfactor). A jump waits for the next full synchronization, when every particle is at the same time, and the block stepper then starts afresh from the state the jump leaves, as it does on a restart.

In a star and black hole encounter, block timesteps should use a few times less CPU time than shared steps. They may need cn5 and cn6 tightened together (say, cn5=cn6=0.0025 rather than 0.01) to conserve energy as well as the shared-step run.

nblockfull and dtforce exist for testing only. nblockfull=1 recomputes every particle’s smoothing length and hydrodynamics at every substep, and dtforce replaces the shared timestep with a fixed value.

A run that starts from rest at t=0 across a discontinuity, such as a shock tube, has \(\dot u=0\) at first, so the cn3 limit cannot act on the first step and that step can be several times too long. nstartramp=1 makes the first step 1/64 of what the limits allow and lets each later step at most double until the limits take over. By default (nstartramp=-1) it is on only for a new shock-tube run, that is a periodic box read with iname='pbx' and nrelax=0, and nstartramp=0 turns it off there too. A restart never ramps.

Restarting

A restartrad.sph file is written every few iterations, overwriting the previous one. If one is present when a run starts, it is used automatically. You can also restart from any snapshot by renaming it to restartrad.sph.

A restart is not only a way to recover from a crash. It is also how settings are changed part-way through a calculation, since sph.input is reread every time: artificial viscosity can be turned off, the output interval changed, or an orbital jump armed once you have seen how the first passage went. See Jumping ahead on a wide orbit.

Note

Output files are numbered per stage, not per rank. A fresh run writes log0.sph and energy0.sph, the first resume writes log1.sph and energy1.sph, and so on. A resumed run does not append to energy0.sph, which matters if you are parsing an energy trace across a restart.