Commit 319025da authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'line_of_sight' into 'master'

Line-of-sight output along the z-axis

See merge request !1050
parents b9af415a 9a4032cc
......@@ -40,6 +40,7 @@ Parameters:
-u, --fof Run Friends-of-Friends algorithm and
black holes seeding.
-x, --velociraptor Run with structure finding.
--line-of-sight Run with line-of-sight outputs.
--limiter Run with time-step limiter.
--sync Run with time-step synchronization
of particles hit by feedback events.
......
......@@ -127,6 +127,7 @@ Parameters:
-u, --fof Run Friends-of-Friends algorithm to
perform black hole seeding.
-x, --velociraptor Run with structure finding.
--line-of-sight Run with line-of-sight outputs.
--limiter Run with time-step limiter.
--sync Run with time-step synchronization
of particles hit by feedback events.
......
......@@ -37,6 +37,7 @@ can be found by typing ``./swift -h``:
-u, --fof Run Friends-of-Friends algorithm to
perform black hole seeding.
-x, --velociraptor Run with structure finding.
--line-of-sight Run with line-of-sight outputs.
--limiter Run with time-step limiter.
--sync Run with time-step synchronization
of particles hit by feedback events.
......
......@@ -7,8 +7,75 @@
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
......
.. Snapshots
Matthieu Schaller, 23rd May 2020
.. _line_of_sight:
Line-of-sights outputs
======================
The line-of-sight outputs are designed to be processed by the
``SpecWizard`` tool (`Theuns et al. 1998 <https://ui.adsabs.harvard.edu/abs/1998MNRAS.301..478T/>`_,
`Tepper-Garcia et al. 2011
<https://ui.adsabs.harvard.edu/abs/2011MNRAS.413..190T/>`_).
TO BE DONE.
......@@ -742,7 +742,8 @@ full section would be:
delta_time: 1.02
invoke_stf: 1
Showing all the parameters for a basic hydro test-case, one would have:
Showing all the parameters for a basic non-cosmological hydro test-case, one
would have:
.. code:: YAML
......@@ -767,6 +768,36 @@ following pages:
* :ref:`Output_list_label` (to have snapshots not evenly spaced in time),
* :ref:`Output_selection_label` (to select what particle fields to write).
.. _Parameters_line_of_sight:
Line-of-sight outputs
---------------------
The ``LineOfSight`` section of the parameter file contains all the options related to
the dump of simulation outputs in the form of HDF5 :ref:`line_of_sight` data to
be processed by the ``SpecWizard`` tool
(See `Theuns et al. 1998 <https://ui.adsabs.harvard.edu/abs/1998MNRAS.301..478T/>`_,
`Tepper-Garcia et al. 2011
<https://ui.adsabs.harvard.edu/abs/2011MNRAS.413..190T/>`_). The parameters are:
.. code:: YAML
LineOfSight:
basename: los
scale_factor_first: 0.02 # Only used when running in cosmological mode
delta_time: 1.02
time_first: 0.01 # Only used when running in non-cosmological mode
output_list_on: 0 # Overwrite the regular output times with a list of output times
num_along_x: 0
num_along_y: 0
num_along_z: 100
allowed_los_range_x: [0, 100.] # Range along the x-axis where LoS along Y or Z are allowed
allowed_los_range_y: [0, 100.] # Range along the y-axis where LoS along X or Z are allowed
allowed_los_range_z: [0, 100.] # Range along the z-axis where LoS along X or Y are allowed
range_when_shooting_down_x: 100. # Range along the x-axis of LoS along x
range_when_shooting_down_y: 100. # Range along the y-axis of LoS along y
range_when_shooting_down_z: 100. # Range along the z-axis of LoS along z
.. _Parameters_fof:
Friends-Of-Friends (FOF)
......
......@@ -23,14 +23,14 @@ difference is the parameter file that will need to be adapted for SWIFT.
HydroSchemes/index
TimeStepping/index
SubgridModels/index
random
Planetary/index
FriendsOfFriends/index
VELOCIraptorInterface/index
LineOfSights/index
EquationOfState/index
ExternalPotentials/index
NewOption/index
Task/index
VELOCIraptorInterface/index
AnalysisTools/index
Logger/index
ImplementationDetails/index
.. Random number generator
Folkert Nobels, 11th of July 2019
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.
......@@ -91,6 +91,15 @@ InitialConditions:
generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
# Parameters of the line-of-sight outputs
LineOfSight:
basename: eagle_los
num_along_x: 0
num_along_y: 0
num_along_z: 100
scale_factor_first: 0.1
delta_time: 1.1
# Impose primoridal metallicity
EAGLEChemistry:
init_abundance_metal: 0.
......@@ -104,6 +113,7 @@ EAGLEChemistry:
init_abundance_Silicon: 0.0
init_abundance_Iron: 0.0
# EAGLE cooling parameters
EAGLECooling:
dir_name: ./coolingtables/
H_reion_z: 7.5 # Planck 2018
......
......@@ -11,7 +11,7 @@ fi
if [ ! -e yieldtables ]
then
echo "Fetching EAGLE yield tables..."
../getEagleYieldtable.sh
../getEagleYieldTable.sh
fi
if [ ! -e coolingtables ]
......
......@@ -91,6 +91,15 @@ InitialConditions:
generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
# Parameters of the line-of-sight outputs
LineOfSight:
basename: eagle_los
num_along_x: 0
num_along_y: 0
num_along_z: 100
scale_factor_first: 0.1
delta_time: 1.1
# Impose primoridal metallicity
EAGLEChemistry:
init_abundance_metal: 0.
......
......@@ -11,7 +11,7 @@ fi
if [ ! -e yieldtables ]
then
echo "Fetching EAGLE yield tables..."
../getEagleYieldtable.sh
../getEagleYieldTable.sh
fi
if [ ! -e coolingtables ]
......
......@@ -89,6 +89,15 @@ InitialConditions:
generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
# Parameters of the line-of-sight outputs
LineOfSight:
basename: eagle_los
num_along_x: 0
num_along_y: 0
num_along_z: 100
scale_factor_first: 0.1
delta_time: 1.1
# Impose primoridal metallicity
EAGLEChemistry:
init_abundance_metal: 0.
......
......@@ -11,7 +11,7 @@ fi
if [ ! -e yieldtables ]
then
echo "Fetching EAGLE yield tables..."
../getEagleYieldtable.sh
../getEagleYieldTable.sh
fi
if [ ! -e coolingtables ]
......
......@@ -88,6 +88,15 @@ InitialConditions:
generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
# Parameters of the line-of-sight outputs
LineOfSight:
basename: eagle_los
num_along_x: 0
num_along_y: 0
num_along_z: 100
scale_factor_first: 0.1
delta_time: 1.1
# Impose primoridal metallicity
EAGLEChemistry:
init_abundance_metal: 0.
......
......@@ -11,7 +11,7 @@ fi
if [ ! -e yieldtables ]
then
echo "Fetching EAGLE yield tables..."
../getEagleYieldtable.sh
../getEagleYieldTable.sh
fi
if [ ! -e coolingtables ]
......
......@@ -62,7 +62,7 @@ Scheduler:
Restarts:
onexit: 1
delta_hours: 2.0
delta_hours: 6.0
# Parameters related to the initial conditions
InitialConditions:
......@@ -73,6 +73,15 @@ InitialConditions:
generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
# Parameters of the line-of-sight outputs
LineOfSight:
basename: los
num_along_x: 0
num_along_y: 0
num_along_z: 100
scale_factor_first: 0.1
delta_time: 1.1
# Quick Lyman-alpha cooling (EAGLE with fixed primoridal Z)
QLACooling:
dir_name: ./coolingtables/
......
......@@ -106,6 +106,7 @@ int main(int argc, char *argv[]) {
struct spart *sparts = NULL;
struct bpart *bparts = NULL;
struct unit_system us;
struct los_props los_properties;
int nr_nodes = 1, myrank = 0;
......@@ -173,6 +174,7 @@ int main(int argc, char *argv[]) {
int with_qla = 0;
int with_eagle = 0;
int with_gear = 0;
int with_line_of_sight = 0;
int verbose = 0;
int nr_threads = 1;
int with_verbose_timers = 0;
......@@ -224,6 +226,8 @@ int main(int argc, char *argv[]) {
NULL, 0, 0),
OPT_BOOLEAN('x', "velociraptor", &with_structure_finding,
"Run with structure finding.", NULL, 0, 0),
OPT_BOOLEAN(0, "line-of-sight", &with_line_of_sight,
"Run with line-of-sight outputs.", NULL, 0, 0),
OPT_BOOLEAN(0, "limiter", &with_timestep_limiter,
"Run with time-step limiter.", NULL, 0, 0),
OPT_BOOLEAN(0, "sync", &with_timestep_sync,
......@@ -514,6 +518,16 @@ int main(int argc, char *argv[]) {
return 1;
}
if (!with_hydro && with_line_of_sight) {
if (myrank == 0) {
argparse_usage(&argparse);
printf(
"\nError: Cannot use line-of-sight outputs without gas, --hydro must "
"be chosen.\n");
}
return 1;
}
/* Let's pin the main thread, now we know if affinity will be used. */
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
if (with_aff &&
......@@ -656,6 +670,18 @@ int main(int argc, char *argv[]) {
}
}
/* Check that we can write the line of sight files by testing if the
* output directory exists and is searchable and writable. */
if (with_line_of_sight) {
char losbasename[PARSER_MAX_LINE_SIZE];
parser_get_param_string(params, "LineOfSight:basename", losbasename);
const char *losdirp = dirname(losbasename);
if (access(losdirp, W_OK | X_OK) != 0) {
error("Cannot write line of sight in directory %s (%s)", losdirp,
strerror(errno));
}
}
/* Prepare the domain decomposition scheme */
struct repartition reparttype;
#ifdef WITH_MPI
......@@ -1063,6 +1089,9 @@ int main(int argc, char *argv[]) {
with_hydro, with_self_gravity, with_star_formation,
with_DM_background_particles, talking, dry_run, nr_nodes);
/* Initialise the line of sight properties. */
if (with_line_of_sight) los_init(s.dim, &los_properties, params);
if (myrank == 0) {
clocks_gettime(&toc);
message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
......@@ -1186,6 +1215,7 @@ int main(int argc, char *argv[]) {
engine_policies |= engine_policy_structure_finding;
if (with_fof) engine_policies |= engine_policy_fof;
if (with_logger) engine_policies |= engine_policy_logger;
if (with_line_of_sight) engine_policies |= engine_policy_line_of_sight;
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
......@@ -1196,7 +1226,7 @@ int main(int argc, char *argv[]) {
&reparttype, &us, &prog_const, &cosmo, &hydro_properties,
&entropy_floor, &gravity_properties, &stars_properties,
&black_holes_properties, &feedback_properties, &mesh, &potential,
&cooling_func, &starform, &chemistry, &fof_properties);
&cooling_func, &starform, &chemistry, &fof_properties, &los_properties);
engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
nr_threads, with_aff, talking, restart_file);
......
......@@ -588,8 +588,8 @@ int main(int argc, char *argv[]) {
/*hydro_properties=*/NULL, /*entropy_floor=*/NULL, &gravity_properties,
/*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
/*feedback_properties=*/NULL, &mesh, /*potential=*/NULL,
/*cooling_func=*/NULL,
/*starform=*/NULL, /*chemistry=*/NULL, &fof_properties);
/*cooling_func=*/NULL, /*starform=*/NULL, /*chemistry=*/NULL,
&fof_properties, /*los_properties=*/NULL);
engine_config(/*restart=*/0, /*fof=*/1, &e, params, nr_nodes, myrank,
nr_threads, with_aff, talking, NULL);
......
......@@ -212,12 +212,29 @@ StructureFinding:
config_file_name: stf_input.cfg # Name of the STF config file.
basename: haloes # Common part of the name of output files.
subdir_per_output: stf # (Optional) Sub-directory (within Snapshots:subdir) to use for each invocation of STF. Defaults to "" (i.e. no sub-directory)
scale_factor_first: 0.92 # (Optional) Scale-factor of the first snaphot (cosmological run)
scale_factor_first: 0.92 # (Optional) Scale-factor of the first structure finding (cosmological run)
time_first: 0.01 # (Optional) Time of the first structure finding output (in internal units).
delta_time: 1.10 # (Optional) Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
output_list_on: 0 # (Optional) Enable the output list
output_list_on: 0 # (Optional) Enable the use of an output list
output_list: stflist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
# Parameters related to the Line-Of-Sight (SpecWizard) outputs
LineOfSight:
basename: los # Basename of the files
scale_factor_first: 0.02 # (Optional) Scale-factor of the first line-of-sight (cosmological run)
time_first: 0.01 # (Optional) Time of the first line-of-sight output (in internal units).
delta_time: 1.02 # (Optional) Time difference between consecutive line-of-sight outputs (in internal units) in simulation time intervals.
output_list_on: 0 # (Optional) Enable the use of an output list
num_along_x: 0 # Number of sight-lines along the x-axis
num_along_y: 0 # Number of sight-lines along the y-axis
num_along_z: 100 # Number of sight-lines along the z-axis
allowed_los_range_x: [0, 100.] # (Optional) Range along the x-axis where LoS along Y or Z are allowed (Defaults to the box size).
allowed_los_range_y: [0, 100.] # (Optional) Range along the y-axis where LoS along X or Z are allowed (Defaults to the box size).
allowed_los_range_z: [0, 100.] # (Optional) Range along the z-axis where LoS along X or Y are allowed (Defaults to the box size).
range_when_shooting_down_x: 100. # (Optional) Range along the x-axis of LoS along x (Defaults to the box size).
range_when_shooting_down_y: 100. # (Optional) Range along the y-axis of LoS along y (Defaults to the box size).
range_when_shooting_down_z: 100. # (Optional) Range along the z-axis of LoS along z (Defaults to the box size).
# Parameters related to the equation of state ------------------------------------------
EoS:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment