.. Implementation details
   Loic Hausammann, 2020
   Matthieu Schaller, 2020

.. _implementation_details:

Implementation details
======================

This section contains technical information about internals of the code.

Random number generator
~~~~~~~~~~~~~~~~~~~~~~~

Often subgrid models require random numbers, for this purpose 
SWIFT has a random number generator. The random number generator
of SWIFT is based on a combination of the standard `rand_r` and `erand48`
random number generators. Since for some applications in cosmology
we want to be able to sample random numbers with a probability lower than 
:math:`10^{-8}`, we could not simply use the 32-bit `rand_r` due to its cut-off
and spacing of around :math:`2^{-32} \approx 2 \times 10^{-10}`.
For the `erand48` algorithm with 48 bits the spacing and cutoff are 
significantly smaller and around :math:`2^{-48} \approx 3.5 \times 10^{-15}`,
so this is very suitable for our applications. 

Reproducible random numbers
---------------------------

In our simulations we want to be able to reproduce the exact same random 
numbers if we do exactly the same simulation twice. Because of this we 
want to seed the random number generator in a reproducible way. In order to do this
we seed with the particle ID of the current particle and the current 
integer simulation time. 

To produce several different types of random numbers we have an additional
argument called the type of random number which is basically the nth random
number for the specified seed, which is added to the particle ID, thus providing
a distinct state per random number type.

If the user wishes to run a simulation with a different set of random number,
an option during the configuration (``--with-random-seed=INT``) is available.
This option simply flip some bits in the initial number composed of the ID and the
current simulation time through the binary operator XOR.

Implementation
--------------

Our random number generator packs the particle ID (plus the random number type) and
the current simulation time as two 64-bit values, plus a constant 16-bit value,
into a 144-bit buffer. This buffer is treated as an array of 9 `uint16` values.

In a first pass we initialize the seed to 0 and run through the 9 `uint16` values,
xor-ing them with the seed and calling `rand_r` on the seed to advance it. Using
`rand_r` with the thus-generated seed, we generate a sequence of 9 16-bit values
and xor them with the original 144-bit buffer.

The 9 bytes of this modified buffer are then used for three passes of `erand48`,
xor-ing the state in the same way as before. `erand48` is then called one final
time with the last state, producing a random double-precision value with a
48-bit mantissa.

What to do if we break the random number generator?
---------------------------------------------------

The most likely case is that the RNG is not strong enough for our application,
in this case we could simply do multiple passes of both shuffling the state and
generating the final value from the state. This increases the computational cost but
will make the random number generator stronger. 

An other case is that we need probabilities that are lower than :math:`1 \times 10^{-17}`,
in this case we simply cannot use our random number generator and for example
need to generate two random numbers in order to probe these low probabilities. 


------------------------------------------------------

Generating new unique IDs
~~~~~~~~~~~~~~~~~~~~~~~~~

When spawning new particles (not converting them) for star formation or other
similar processes, the code needs to create new unique particle IDs. This is
implemented in the file ``space_unique_id.c`` and can be switched on/off in the
star formation file ``star_formation_struct.h`` by setting the variable
``star_formation_need_unique_id`` to 1 or 0.

The generation of new IDs is done by computing the maximal ID present in the
initial condition (across all particle types) and then attributing two batches
of new, unused IDs to each MPI rank.  The size of each batch is computed in the
same way as the count of extra particles in order to ensure that we will have
enough available IDs between two tree rebuilds (where the extra particles are
regenerated).

When a new ID is requested, the next available ID in the first batch is
returned. If the last available ID in the first batch is requested, we switch to
the next batch of IDs and flag it for regeneration at the next rebuild time.  If
the second batch is also fully used, the code will exit with an error message
[#f1]_. At each tree-rebuild steps, the ranks will request a new batch if
required and make sure the batches are unique across all MPI ranks.

As we are using the maximal ID from the ICs, nothing can be done against the user
providing the maximal integer possible as an ID (that can for instance be the
case in some of the EAGLE ICs as the ID encode their Lagrangian position on a
Peano-Hilbert curve). 


.. [#f1] Thanks to the size of the fresh ID batches, the code should run out of
	 extra particles before reaching this point and triggered a new rebuild
	 if this is allowed by the star formation scheme.
