From c549d4cb978a6aafa071e1ba5189adf436f5598f Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 8 Jun 2020 12:20:45 +0200
Subject: [PATCH] Revert "Revert "Merge branch 'master' into 'updated_MAC'""
This reverts commit 6f6c3499f88dd583de1e01fe709d95a95a142e75.
---
README | 1 +
README.md | 1 +
configure.ac | 11 +
doc/RTD/source/CommandLineOptions/index.rst | 1 +
.../source/ImplementationDetails/index.rst | 69 +-
doc/RTD/source/LineOfSights/index.rst | 16 +
.../ParameterFiles/parameter_description.rst | 33 +-
doc/RTD/source/index.rst | 4 +-
doc/RTD/source/random.rst | 65 --
examples/EAGLE_ICs/EAGLE_12/eagle_12.yml | 42 +-
examples/EAGLE_ICs/EAGLE_12/run.sh | 2 +-
examples/EAGLE_ICs/EAGLE_25/eagle_25.yml | 41 +-
examples/EAGLE_ICs/EAGLE_25/run.sh | 2 +-
examples/EAGLE_ICs/EAGLE_50/eagle_50.yml | 41 +-
examples/EAGLE_ICs/EAGLE_50/run.sh | 2 +-
examples/EAGLE_ICs/EAGLE_50_low_res/README | 3 +
.../EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml | 209 ++++
examples/EAGLE_ICs/EAGLE_50_low_res/getIC.sh | 2 +
.../EAGLE_50_low_res/output_list.txt | 38 +
examples/EAGLE_ICs/EAGLE_50_low_res/run.sh | 35 +
.../vrconfig_3dfof_subhalos_SO_hydro.cfg | 191 ++++
examples/EAGLE_ICs/README | 5 +-
examples/EAGLE_low_z/EAGLE_12/eagle_12.yml | 32 +-
examples/EAGLE_low_z/EAGLE_25/eagle_25.yml | 32 +-
examples/EAGLE_low_z/EAGLE_50/eagle_50.yml | 32 +-
examples/EAGLE_low_z/EAGLE_6/eagle_6.yml | 32 +-
examples/QuickLymanAlpha/L050N0752/qla_50.yml | 11 +-
examples/main.c | 32 +-
examples/main_fof.c | 4 +-
examples/parameter_example.yml | 57 +-
src/Makefile.am | 6 +-
src/black_holes/Default/black_holes.h | 3 +-
src/black_holes/Default/black_holes_iact.h | 6 +-
src/black_holes/EAGLE/black_holes.h | 216 ++--
src/black_holes/EAGLE/black_holes_iact.h | 188 ++-
src/black_holes/EAGLE/black_holes_io.h | 79 +-
.../EAGLE/black_holes_parameters.h | 2 +
src/black_holes/EAGLE/black_holes_part.h | 33 +-
.../EAGLE/black_holes_properties.h | 130 ++-
src/cooling/EAGLE/cooling_io.h | 1 +
src/cooling/const_du/cooling.h | 1 +
src/cooling/const_du/cooling_io.h | 1 +
src/cooling/const_lambda/cooling.h | 1 +
src/cooling/const_lambda/cooling_io.h | 1 +
src/cooling/none/cooling_io.h | 1 +
src/distributed_io.c | 1 +
src/engine.c | 157 ++-
src/engine.h | 18 +-
src/hydro/AnarchyPU/hydro.h | 37 +-
src/hydro/Default/hydro.h | 27 +-
src/hydro/Gadget2/hydro.h | 34 +-
src/hydro/Minimal/hydro.h | 32 +-
src/hydro/PressureEnergy/hydro.h | 37 +-
.../PressureEnergyMorrisMonaghanAV/hydro.h | 37 +-
src/hydro/PressureEntropy/hydro.h | 22 +-
src/hydro/SPHENIX/hydro.h | 52 +-
src/hydro/SPHENIX/hydro_parameters.h | 2 +-
src/line_of_sight.c | 1006 +++++++++++++++++
src/line_of_sight.h | 122 ++
src/parallel_io.c | 1 +
src/parser.h | 2 +-
src/part.c | 2 +
src/part.h | 1 +
src/runner_black_holes.c | 3 -
src/runner_doiact_functions_black_holes.h | 16 +-
src/runner_ghost.c | 8 +-
src/serial_io.c | 1 +
src/sincos.h | 84 ++
src/single_io.c | 1 +
src/swift.h | 1 +
src/task.c | 4 +
tools/analyse_runtime.py | 1 +
72 files changed, 2950 insertions(+), 474 deletions(-)
create mode 100644 doc/RTD/source/LineOfSights/index.rst
delete mode 100644 doc/RTD/source/random.rst
create mode 100644 examples/EAGLE_ICs/EAGLE_50_low_res/README
create mode 100644 examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
create mode 100755 examples/EAGLE_ICs/EAGLE_50_low_res/getIC.sh
create mode 100644 examples/EAGLE_ICs/EAGLE_50_low_res/output_list.txt
create mode 100755 examples/EAGLE_ICs/EAGLE_50_low_res/run.sh
create mode 100644 examples/EAGLE_ICs/EAGLE_50_low_res/vrconfig_3dfof_subhalos_SO_hydro.cfg
create mode 100644 src/line_of_sight.c
create mode 100644 src/line_of_sight.h
create mode 100644 src/sincos.h
diff --git a/README b/README
index 97e3d48276..8587a13b35 100644
--- a/README
+++ b/README
@@ -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.
diff --git a/README.md b/README.md
index 99af7ba4f7..7223f4719a 100644
--- a/README.md
+++ b/README.md
@@ -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.
diff --git a/configure.ac b/configure.ac
index 2181f2167b..46d9fb2a8a 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1383,6 +1383,17 @@ if test "$ax_cv_c_compiler_vendor" = "clang"; then
AC_CHECK_LIB([m],[__exp10f], [AC_DEFINE([HAVE___EXP10F],1,[The __exp10f function is present.])])
fi
+# Check if we have native sincos and sincosf functions. If not failback to our
+# implementations. On Apple/CLANG we have __sincos, so also check for that
+# if the compiler is clang.
+AC_CHECK_LIB([m],[sincos], [AC_DEFINE([HAVE_SINCOS],1,[The sincos function is present.])])
+AC_CHECK_LIB([m],[sincosf], [AC_DEFINE([HAVE_SINCOSF],1,[The sincosf function is present.])])
+if test "$ax_cv_c_compiler_vendor" = "clang"; then
+ AC_CHECK_LIB([m],[__sincos], [AC_DEFINE([HAVE___SINCOS],1,[The __sincos function is present.])])
+ AC_CHECK_LIB([m],[__sincosf], [AC_DEFINE([HAVE___SINCOSF],1,[The __sincosf function is present.])])
+fi
+
+
# Add warning flags by default, if these can be used. Option =error adds
# -Werror to GCC, clang and Intel. Note do this last as compiler tests may
# become errors, if that's an issue don't use CFLAGS for these, use an AC_SUBST().
diff --git a/doc/RTD/source/CommandLineOptions/index.rst b/doc/RTD/source/CommandLineOptions/index.rst
index a5f9e7eff5..2de6c6bb25 100644
--- a/doc/RTD/source/CommandLineOptions/index.rst
+++ b/doc/RTD/source/CommandLineOptions/index.rst
@@ -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.
diff --git a/doc/RTD/source/ImplementationDetails/index.rst b/doc/RTD/source/ImplementationDetails/index.rst
index dda984293e..f56126f1ae 100644
--- a/doc/RTD/source/ImplementationDetails/index.rst
+++ b/doc/RTD/source/ImplementationDetails/index.rst
@@ -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
diff --git a/doc/RTD/source/LineOfSights/index.rst b/doc/RTD/source/LineOfSights/index.rst
new file mode 100644
index 0000000000..42558cee4e
--- /dev/null
+++ b/doc/RTD/source/LineOfSights/index.rst
@@ -0,0 +1,16 @@
+.. 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.
+
diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index ee854dcc89..61cd9c68bc 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -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)
diff --git a/doc/RTD/source/index.rst b/doc/RTD/source/index.rst
index ae75332e7a..381f857728 100644
--- a/doc/RTD/source/index.rst
+++ b/doc/RTD/source/index.rst
@@ -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
diff --git a/doc/RTD/source/random.rst b/doc/RTD/source/random.rst
deleted file mode 100644
index 4ca7e0b5bf..0000000000
--- a/doc/RTD/source/random.rst
+++ /dev/null
@@ -1,65 +0,0 @@
-.. 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.
-
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 7945009d7d..ca539e3dcf 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -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
@@ -180,15 +190,23 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1.0 # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/examples/EAGLE_ICs/EAGLE_12/run.sh b/examples/EAGLE_ICs/EAGLE_12/run.sh
index e908b552de..78dc4692fa 100755
--- a/examples/EAGLE_ICs/EAGLE_12/run.sh
+++ b/examples/EAGLE_ICs/EAGLE_12/run.sh
@@ -11,7 +11,7 @@ fi
if [ ! -e yieldtables ]
then
echo "Fetching EAGLE yield tables..."
- ../getEagleYieldtable.sh
+ ../getEagleYieldTable.sh
fi
if [ ! -e coolingtables ]
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 0ca0f37bcf..e598e96436 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -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.
@@ -181,15 +190,23 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1.0 # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/examples/EAGLE_ICs/EAGLE_25/run.sh b/examples/EAGLE_ICs/EAGLE_25/run.sh
index 17dcb98069..5b5ca5a4e9 100755
--- a/examples/EAGLE_ICs/EAGLE_25/run.sh
+++ b/examples/EAGLE_ICs/EAGLE_25/run.sh
@@ -11,7 +11,7 @@ fi
if [ ! -e yieldtables ]
then
echo "Fetching EAGLE yield tables..."
- ../getEagleYieldtable.sh
+ ../getEagleYieldTable.sh
fi
if [ ! -e coolingtables ]
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index b4aa82a79c..825669a8e2 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -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.
@@ -179,15 +188,23 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1.0 # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
\ No newline at end of file
diff --git a/examples/EAGLE_ICs/EAGLE_50/run.sh b/examples/EAGLE_ICs/EAGLE_50/run.sh
index 18e8f5cfce..2afd0de07b 100755
--- a/examples/EAGLE_ICs/EAGLE_50/run.sh
+++ b/examples/EAGLE_ICs/EAGLE_50/run.sh
@@ -11,7 +11,7 @@ fi
if [ ! -e yieldtables ]
then
echo "Fetching EAGLE yield tables..."
- ../getEagleYieldtable.sh
+ ../getEagleYieldTable.sh
fi
if [ ! -e coolingtables ]
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/README b/examples/EAGLE_ICs/EAGLE_50_low_res/README
new file mode 100644
index 0000000000..3b68fbf9ca
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/README
@@ -0,0 +1,3 @@
+Initial conditions corresponding to the 50 Mpc volume
+of the EAGLE suite at 8x lower resolution.
+The ICs only contain DM particles. The gas particles will be generated in SWIFT.
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
new file mode 100644
index 0000000000..8e0f6b6114
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
@@ -0,0 +1,209 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: EAGLE-L0050N0376-Ref
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98841e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ a_begin: 0.0078125 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_m: 0.307 # Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: eagle # Common part of the name of output files
+ output_list_on: 1
+ output_list: ./output_list.txt
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.02
+ scale_factor_first: 0.05
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ mesh_side_length: 128
+ comoving_DM_softening: 0.006640 # Comoving softening for DM (6.67 ckpc)
+ max_physical_DM_softening: 0.002600 # Physical softening for DM (2.60 pkpc)
+ comoving_baryon_softening: 0.003580 # Comoving softening for baryons (3.58 ckpc)
+ max_physical_baryon_softening: 0.001400 # Physical softening for baryons (1.40 pkpc)
+ dithering: 0
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing length in units of softening.
+ h_max: 0.5 # Maximal smoothing length in co-moving internal units.
+ CFL_condition: 0.2 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100.0 # (internal units)
+ initial_temperature: 268.7 # (internal units)
+ particle_splitting: 1 # Particle splitting is ON
+ particle_splitting_mass_threshold: 5.6e-3 # (internal units, i.e. 5.6e7 Msun ~ 4x initial gas particle mass)
+
+# Parameters of the stars neighbour search
+Stars:
+ resolution_eta: 1.1642 # Target smoothing length in units of the mean inter-particle separation
+ h_tolerance: 7e-3
+
+# Parameters for the Friends-Of-Friends algorithm
+FOF:
+ basename: fof_output # Filename for the FOF outputs.
+ min_group_size: 32 # The minimum no. of particles required for a group.
+ linking_length_ratio: 0.2 # Linking length in units of the main inter-particle separation.
+ black_hole_seed_halo_mass_Msun: 1.5e10 # Minimal halo mass in which to seed a black hole (in solar masses).
+ scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
+ delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+
+Scheduler:
+ max_top_level_cells: 16
+ cell_split_size: 200
+
+Restarts:
+ onexit: 1
+ delta_hours: 6.0
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: EAGLE_L0050N0376_ICs.hdf5
+ periodic: 1
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
+ cleanup_velocity_factors: 1 # Remove the sqrt(a) factor in the velocities inherited from Gadget
+ 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.
+ init_abundance_Hydrogen: 0.752
+ init_abundance_Helium: 0.248
+ init_abundance_Carbon: 0.0
+ init_abundance_Nitrogen: 0.0
+ init_abundance_Oxygen: 0.0
+ init_abundance_Neon: 0.0
+ init_abundance_Magnesium: 0.0
+ init_abundance_Silicon: 0.0
+ init_abundance_Iron: 0.0
+
+# EAGLE cooling parameters
+EAGLECooling:
+ dir_name: ./coolingtables/
+ H_reion_z: 7.5 # Planck 2018
+ H_reion_eV_p_H: 2.0
+ He_reion_z_centre: 3.5
+ He_reion_z_sigma: 0.5
+ He_reion_eV_p_H: 2.0
+
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+ EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+ EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ min_over_density: 57.7 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ EOS_entropy_margin_dex: 0.5 # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+ threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+
+# EAGLE feedback model
+EAGLEFeedback:
+ use_SNII_feedback: 1 # Global switch for SNII thermal (stochastic) feedback.
+ use_SNIa_feedback: 1 # Global switch for SNIa thermal (continuous) feedback.
+ use_AGB_enrichment: 1 # Global switch for enrichement from AGB stars.
+ use_SNII_enrichment: 1 # Global switch for enrichement from SNII stars.
+ use_SNIa_enrichment: 1 # Global switch for enrichement from SNIa stars.
+ filename: ./yieldtables/ # Path to the directory containing the EAGLE yield tables.
+ IMF_min_mass_Msun: 0.1 # Minimal stellar mass considered for the Chabrier IMF in solar masses.
+ IMF_max_mass_Msun: 100.0 # Maximal stellar mass considered for the Chabrier IMF in solar masses.
+ SNII_min_mass_Msun: 8.0 # Minimal mass considered for SNII stars in solar masses.
+ SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII stars in solar masses.
+ SNII_sampled_delay: 1 # Sample the SNII lifetimes to do feedback.
+ SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
+ SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
+ SNII_energy_fraction_n_0_H_p_cm3: 1.4588 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
+ SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
+ SNII_energy_fraction_n_n: 0.8686 # Power-law for the birth density dependance of the SNII energy fraction.
+ SNIa_DTD: Exponential # Functional form of the SNIa delay time distribution.
+ SNIa_DTD_delay_Gyr: 0.04 # Stellar age after which SNIa start in Gyr (40 Myr corresponds to stars ~ 8 Msun).
+ SNIa_DTD_exp_timescale_Gyr: 2.0 # Time-scale of the exponential decay of the SNIa rates in Gyr.
+ SNIa_DTD_exp_norm_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
+ SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
+ AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
+ SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
+ SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
+ SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
+ SNII_yield_factor_Nitrogen: 1.0 # (Optional) Correction factor to apply to the Nitrogen yield from the SNII channel.
+ SNII_yield_factor_Oxygen: 1.0 # (Optional) Correction factor to apply to the Oxygen yield from the SNII channel.
+ SNII_yield_factor_Neon: 1.0 # (Optional) Correction factor to apply to the Neon yield from the SNII channel.
+ SNII_yield_factor_Magnesium: 2.0 # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel.
+ SNII_yield_factor_Silicon: 1.0 # (Optional) Correction factor to apply to the Silicon yield from the SNII channel.
+ SNII_yield_factor_Iron: 0.5 # (Optional) Correction factor to apply to the Iron yield from the SNII channel.
+
+# EAGLE AGN model
+EAGLEAGN:
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/getIC.sh b/examples/EAGLE_ICs/EAGLE_50_low_res/getIC.sh
new file mode 100755
index 0000000000..fd1f7120ba
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/getIC.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/EAGLE_ICs/EAGLE_L0050N0376_ICs.hdf5
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/output_list.txt b/examples/EAGLE_ICs/EAGLE_50_low_res/output_list.txt
new file mode 100644
index 0000000000..592ab8483d
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/output_list.txt
@@ -0,0 +1,38 @@
+# Redshift
+18.08
+15.28
+13.06
+11.26
+9.79
+8.57
+7.54
+6.67
+5.92
+5.28
+4.72
+4.24
+3.81
+3.43
+3.09
+2.79
+2.52
+2.28
+2.06
+1.86
+1.68
+1.51
+1.36
+1.21
+1.08
+0.96
+0.85
+0.74
+0.64
+0.55
+0.46
+0.37
+0.29
+0.21
+0.14
+0.07
+0.00
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/run.sh b/examples/EAGLE_ICs/EAGLE_50_low_res/run.sh
new file mode 100755
index 0000000000..3022d6f44f
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/run.sh
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e EAGLE_L0050N0376_ICs.hdf5 ]
+then
+ echo "Fetching initial conditions for the EAGLE 50Mpc low-res. example..."
+ ./getIC.sh
+fi
+
+# Grab the cooling and yield tables if they are not present.
+if [ ! -e yieldtables ]
+then
+ echo "Fetching EAGLE yield tables..."
+ ../getEagleYieldTable.sh
+fi
+
+if [ ! -e coolingtables ]
+then
+ echo "Fetching EAGLE cooling tables..."
+ ../getEagleCoolingTable.sh
+fi
+
+# The following run-time options are broken down by line as:
+# Basic run-time options
+# Create and run with stars
+# Radiative options - run with cooling and stellar feedback
+# Run with the time-step limiter required to capture feedback
+# Run with black holes - fof is needed for the seeding
+# Threading options - run with threads and pinning (latter not required but improves performance)
+# The corresponding parameter file for this run
+
+../../swift \
+ --cosmology --eagle \
+ --threads=16 --pin \
+ eagle_50.yml
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/vrconfig_3dfof_subhalos_SO_hydro.cfg b/examples/EAGLE_ICs/EAGLE_50_low_res/vrconfig_3dfof_subhalos_SO_hydro.cfg
new file mode 100644
index 0000000000..8590cbf5bc
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/vrconfig_3dfof_subhalos_SO_hydro.cfg
@@ -0,0 +1,191 @@
+#Configuration file for analysing Hydro
+#runs 3DFOF + substructure algorithm, demands subhalos and FOF halos be self-bound, calculates many properties
+#Units currently set to take in as input, Mpc, 1e10 solar masses, km/s, output in same units
+#To set temporally unique halo ids, alter Snapshot_value=SNAP to appropriate value. Ie: for snapshot 12, change SNAP to 12
+
+################################
+#input options
+#set up to use SWIFT HDF input, load gas, star, bh and dark matter
+################################
+HDF_name_convention=6 #HDF SWIFT naming convention
+Input_includes_dm_particle=1 #include dark matter particles in hydro input
+Input_includes_gas_particle=1 #include gas particles in hydro input
+Input_includes_star_particle=1 #include star particles in hydro input
+Input_includes_bh_particle=1 #include bh particles in hydro input
+Input_includes_wind_particle=0 #include wind particles in hydro input (used by Illustris and moves particle type 0 to particle type 3 when decoupled from hydro forces). Here shown as example
+Input_includes_tracer_particle=0 #include tracer particles in hydro input (used by Illustris). Here shown as example
+Input_includes_extradm_particle=0 #include extra dm particles stored in particle type 2 and type 3, useful for zooms
+
+Halo_core_phase_merge_dist=0.25 #merge substructures if difference in dispersion normalised distance is < this value
+Apply_phase_merge_to_host=1 #merge substructures with background if centrally located and phase-distance is small
+
+#units conversion from input input to desired internal unit
+Length_input_unit_conversion_to_output_unit=1.0 #default code unit,
+Velocity_input_unit_conversion_to_output_unit=1.0 #default velocity unit,
+Mass_input_unit_conversion_to_output_unit=1.0 #default mass unit,
+#assumes input is in 1e10 msun, Mpc and km/s and output units are the same
+Gravity=43.0211349 #for 1e10 Msun, km/s and Mpc
+Hubble_unit=100.0 # assuming units are km/s and Mpc, then value of Hubble in km/s/Mpc
+#converting hydro quantities
+Stellar_age_input_is_cosmological_scalefactor=1
+Metallicity_input_unit_conversion_to_output_unit=1.0
+Stellar_age_input_unit_conversion_to_output_unit=1.0
+Star_formation_rate_input_unit_conversion_to_output_unit=1.0
+
+#set the units of the output by providing conversion to a defined unit
+#conversion of output length units to kpc
+Length_unit_to_kpc=1000.0
+#conversion of output velocity units to km/s
+Velocity_to_kms=1.0
+#conversion of output mass units to solar masses
+Mass_to_solarmass=1.0e10
+#1 / 0.012
+Metallicity_to_solarmetallicity=83.33
+Star_formation_rate_to_solarmassperyear=97.78
+Stellar_age_to_yr=1.0
+#ensures that output is physical and not comoving distances per little h
+Comoving_units=0
+
+#sets the total buffer size in bytes used to store temporary particle information
+#of mpi read threads before they are broadcast to the appropriate waiting non-read threads
+#if not set, default value is equivalent to 1e6 particles per mpi process, quite large
+#but significantly minimises the number of send/receives
+#in this example the buffer size is roughly that for a send/receive of 10000 particles
+#for 100 mpi processes
+MPI_particle_total_buf_size=100000000
+
+################################
+#search related options
+################################
+
+#how to search a simulation
+Particle_search_type=1 #search dark matter particles only
+#for baryon search
+Baryon_searchflag=2 #if 1 search for baryons separately using phase-space search when identifying substructures, 2 allows special treatment in field FOF linking and phase-space substructure search, 0 treat the same as dark matter particles
+#for search for substruture
+Search_for_substructure=1 #if 0, end search once field objects are found
+#also useful for zoom simulations or simulations of individual objects, setting this flag means no field structure search is run
+Singlehalo_search=0 #if file is single halo in which one wishes to search for substructure. Here disabled.
+#additional option for field haloes
+Keep_FOF=0 #if field 6DFOF search is done, allows to keep structures found in 3DFOF (can be interpreted as the inter halo stellar mass when only stellar search is used).\n
+
+#minimum size for structures
+Minimum_size=20 #min 20 particles
+Minimum_halo_size=32 #if field halos have different minimum sizes, otherwise set to -1.
+
+#for field fof halo search
+FoF_Field_search_type=5 #5 3DFOF search for field halos, 4 for 6DFOF clean up of field halos, 3 for 6DFOF with velocity scale distinct for each initial 3D FOF candidate
+Halo_3D_linking_length=0.20
+
+#for mean field estimates and local velocity density distribution funciton estimator related quantiites, rarely need to change this
+Local_velocity_density_approximate_calculation=1 #calculates velocity density using approximative (and quicker) near neighbour search
+Cell_fraction = 0.01 #fraction of field fof halo used to determine mean velocity distribution function. Typical values are ~0.005-0.02
+Grid_type=1 #normal entropy based grid, shouldn't have to change
+Nsearch_velocity=32 #number of velocity neighbours used to calculate local velocity distribution function. Typial values are ~32
+Nsearch_physical=256 #numerof physical neighbours from which the nearest velocity neighbour set is based. Typical values are 128-512
+
+#for substructure search, rarely ever need to change this
+FoF_search_type=1 #default phase-space FOF search. Don't really need to change
+Iterative_searchflag=1 #iterative substructure search, for substructure find initial candidate substructures with smaller linking lengths then expand search region
+Outlier_threshold=2.5 #outlier threshold for a particle to be considered residing in substructure, that is how dynamically distinct a particle is. Typical values are >2
+Substructure_physical_linking_length=0.10
+Velocity_ratio=2.0 #ratio of speeds used in phase-space FOF
+Velocity_opening_angle=0.10 #angle between velocities. 18 degrees here, typical values are ~10-30
+Velocity_linking_length=0.20 #where scaled by structure dispersion
+Significance_level=1.0 #how significant a substructure is relative to Poisson noise. Values >= 1 are fine.
+
+#for iterative substructure search, rarely ever need to change this
+Iterative_threshold_factor=1.0 #change in threshold value when using iterative search. Here no increase in threshold if iterative or not
+Iterative_linking_length_factor=2.0 #increase in final linking final iterative substructure search
+Iterative_Vratio_factor=1.0 #change in Vratio when using iterative search. no change in vratio
+Iterative_ThetaOp_factor=1.0 #change in velocity opening angle. no change in velocity opening angle
+
+#for checking for halo merger remnants, which are defined as large, well separated phase-space density maxima
+Halo_core_search=2 # searches for separate 6dfof cores in field haloes, and then more than just flags halo as merging, assigns particles to each merging "halo". 2 is full separation, 1 is flagging, 0 is off
+#if searching for cores, linking lengths. likely does not need to change much
+Use_adaptive_core_search=0 #calculate dispersions in configuration & vel space to determine linking lengths
+Use_phase_tensor_core_growth=2 #use full stepped phase-space tensor assignment
+Halo_core_ellx_fac=0.7 #how linking lengths are changed when searching for local 6DFOF cores,
+Halo_core_ellv_fac=2.0 #how velocity lengths based on dispersions are changed when searching for local 6DFOF cores
+Halo_core_ncellfac=0.005 #fraction of total halo particle number setting min size of a local 6DFOF core
+Halo_core_num_loops=8 #number of loops to iteratively search for cores
+Halo_core_loop_ellx_fac=0.75 #how much to change the configuration space linking per iteration
+Halo_core_loop_ellv_fac=1.0 #how much to change the velocity space linking per iteration
+Halo_core_loop_elln_fac=1.2 #how much to change the min number of particles per iteration
+Halo_core_phase_significance=2.0 #how significant a core must be in terms of dispersions (sigma) significance
+
+################################
+#Unbinding options (VELOCIraptor is able to accurately identify tidal debris so particles need not be bound to a structure)
+################################
+
+#unbinding related items
+Unbind_flag=1 #run unbinding
+#objects must have particles that meet the allowed kinetic to potential ratio AND also have some total fraction that are completely bound.
+Unbinding_type=0
+#alpha factor used to determine whether particle is "bound" alaph*T+W<0. For standard subhalo catalogues use >0.9 but if interested in tidal debris 0.2-0.5
+Allowed_kinetic_potential_ratio=0.95
+Min_bound_mass_frac=0.65 #minimum bound mass fraction
+#run unbinding of field structures, aka halos. This is useful for sams and 6DFOF halos but may not be useful if interested in 3DFOF mass functions.
+Bound_halos=0
+#don't keep background potential when unbinding
+Keep_background_potential=1
+#use all particles to determine velocity frame for unbinding
+Frac_pot_ref=1.0
+Min_npot_ref=20
+#reference frame only meaningful if calculating velocity frame using subset of particles in object. Can use radially sorted fraction of particles about minimum potential or centre of mass
+Kinetic_reference_frame_type=0
+Unbinding_max_unbound_removal_fraction_per_iteration=0.5
+Unbinding_max_unbound_fraction=0.95
+Unbinding_max_unbound_fraction_allowed=0.005
+
+################################
+#Calculation of properties related options
+################################
+Virial_density=500 #user defined virial overdensity. Note that 200 rho_c, 200 rho_m and BN98 are already calculated.
+#when calculating properties, for field objects calculate inclusive masses
+Inclusive_halo_masses=3 #calculate inclusive masses for halos using full Spherical overdensity apertures
+#ensures that output is physical and not comoving distances per little h
+Comoving_units=0
+#calculate more (sub)halo properties (like angular momentum in spherical overdensity apertures, both inclusive and exclusive)
+Extensive_halo_properties_output=1
+Extensive_gas_properties_output=1
+Extensive_star_properties_output=1
+#calculate aperture masses
+Calculate_aperture_quantities=1
+Number_of_apertures=5
+Aperture_values_in_kpc=5,10,30,50,100,
+Number_of_projected_apertures=5
+Projected_aperture_values_in_kpc=5,10,30,50,100,
+#calculate radial profiles
+Calculate_radial_profiles=1
+Number_of_radial_profile_bin_edges=20
+#default radial normalisation log rad bins, normed by R200crit, Integer flag of 0 is log bins and R200crit norm.
+Radial_profile_norm=0
+Radial_profile_bin_edges=-2.,-1.87379263,-1.74758526,-1.62137789,-1.49517052,-1.36896316,-1.24275579,-1.11654842,-0.99034105,-0.86413368,-0.73792631,-0.61171894,-0.48551157,-0.3593042,-0.23309684,-0.10688947,0.0193179,0.14552527,0.27173264,0.39794001,
+Iterate_cm_flag=0 #do not interate to determine centre-of-mass
+Sort_by_binding_energy=1 #sort particles by binding energy
+Reference_frame_for_properties=2 #use the minimum potential as reference frame about which to calculate properties
+
+################################
+#output related
+################################
+
+Write_group_array_file=0 #do not write a group array file
+Separate_output_files=0 #do not separate output into field and substructure files similar to subfind
+Binary_output=2 #Use HDF5 output (binary output 1, ascii 0, and HDF 2)
+#output particles residing in the spherical overdensity apertures of halos, only the particles exclusively belonging to halos
+Spherical_overdensity_halo_particle_list_output=1
+
+#halo ids are adjusted by this value * 1000000000000 (or 1000000 if code compiled with the LONGINTS option turned off)
+#to ensure that halo ids are temporally unique. So if you had 100 snapshots, for snap 100 set this to 100 and 100*1000000000000 will
+#be added to the halo id as set for this snapshot, so halo 1 becomes halo 100*1000000000000+1 and halo 1 of snap 0 would just have ID=1
+
+#ALTER THIS as part of a script to get temporally unique ids
+Snapshot_value=SNAP
+
+################################
+#other options
+################################
+Verbose=0 #how talkative do you want the code to be, 0 not much, 1 a lot, 2 chatterbox
+
+
diff --git a/examples/EAGLE_ICs/README b/examples/EAGLE_ICs/README
index c6ed9b35f6..fa0deffdf8 100644
--- a/examples/EAGLE_ICs/README
+++ b/examples/EAGLE_ICs/README
@@ -4,7 +4,7 @@ and phases are the same as used in the original suite. The only
difference is the file format, adapted for SWIFT.
Compared to the original EAGLE runs of Schaye et al. 2015),
-the following changes have been made:
+the following changes have been made (at standard resolution):
- The dark matter softening lengths have been increased to 1.3
pkpc and 3.32 ckpc. The comoving baryon softening lengths have
@@ -23,6 +23,9 @@ the following changes have been made:
distribution and not using a fixed delay of 30 Myr any more.
- The black hole accretion rate is now limited to 100% of the
Eddington rate (from 100/h = 150 %).
+ - The circular velocity threshold for BH mergers is measured
+ at the actual distance linking the BHs not at the kernel support
+ length any more.
The scripts in this directory download the tables required to
run the EAGLE model. Plotting scripts are also provided
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index f77d2d177d..8c2b1b36f0 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -178,15 +178,23 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 32a28d2697..622a60dd1c 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -182,15 +182,23 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index b5649c9dfe..005fcf5660 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -173,15 +173,23 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 9e5f359366..b5c049fcc7 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -188,15 +188,23 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/examples/QuickLymanAlpha/L050N0752/qla_50.yml b/examples/QuickLymanAlpha/L050N0752/qla_50.yml
index 32a5fbf66d..9eeac77e73 100644
--- a/examples/QuickLymanAlpha/L050N0752/qla_50.yml
+++ b/examples/QuickLymanAlpha/L050N0752/qla_50.yml
@@ -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/
diff --git a/examples/main.c b/examples/main.c
index 5f7ac3ff9e..bf23cff630 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -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);
diff --git a/examples/main_fof.c b/examples/main_fof.c
index 6204d2d9ca..6b5841eed7 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -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);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 4f10f6bc0c..759349e1a8 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -216,12 +216,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:
@@ -482,15 +499,27 @@ GEARFeedback:
# EAGLE AGN model
EAGLEAGN:
- subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
- eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
- viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
- radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
- coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
- AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
- AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
- max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
- threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
- threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
+ use_subgrid_mass_from_ics: 1 # (Optional) Use subgrid masses specified in ICs [1, default], or initialise them to particle masses [0]?
+ with_subgrid_mass_check: 1 # (Optional) Verify that initial black hole subgrid masses are positive [1, default]. Only used if use_subgrid_mass_from_ics is 1.
+ multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
+ with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
+ viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
+ radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
+ coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
+ AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+ AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
+ max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
+ with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
+ max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
+ min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
+ set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
+ reposition_coefficient_upsilon: 0.0001 # Repositioning speed normalisation [km/s/M_sun]. Only meaningful if set_reposition_speed is 1.
+ reposition_exponent_xi: 1.0 # (Optional) Scaling of repositioning velocity with BH subgrid mass (default: 1.0, linear). Only meaningful if set_reposition_speed is 1.
+ threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
+ threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
+ merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
+ merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
diff --git a/src/Makefile.am b/src/Makefile.am
index facf21d15c..457c903b4d 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -52,14 +52,14 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \
logger_io.h tracers_io.h tracers.h tracers_struct.h star_formation_io.h fof.h fof_struct.h fof_io.h \
- multipole.h multipole_accept.h multipole_struct.h binomial.h integer_power.h \
+ multipole.h multipole_accept.h multipole_struct.h binomial.h integer_power.h sincos.h \
star_formation_struct.h star_formation.h star_formation_iact.h \
star_formation_logger.h star_formation_logger_struct.h \
pressure_floor.h pressure_floor_struct.h pressure_floor_iact.h \
velociraptor_struct.h velociraptor_io.h random.h memuse.h mpiuse.h memuse_rnodes.h \
black_holes.h black_holes_io.h black_holes_properties.h black_holes_struct.h \
feedback.h feedback_struct.h feedback_properties.h task_order.h \
- space_unique_id.h
+ space_unique_id.h line_of_sight.h
# source files for EAGLE cooling
QLA_COOLING_SOURCES =
@@ -109,7 +109,7 @@ AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c
collectgroup.c hydro_space.c equation_of_state.c \
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
outputlist.c velociraptor_dummy.c logger_io.c memuse.c mpiuse.c memuse_rnodes.c fof.c \
- hashmap.c pressure_floor.c space_unique_id.c \
+ hashmap.c pressure_floor.c space_unique_id.c line_of_sight.c \
$(QLA_COOLING_SOURCES) \
$(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES) \
$(GRACKLE_COOLING_SOURCES) $(GEAR_FEEDBACK_SOURCES)
diff --git a/src/black_holes/Default/black_holes.h b/src/black_holes/Default/black_holes.h
index 929524ebff..e794228621 100644
--- a/src/black_holes/Default/black_holes.h
+++ b/src/black_holes/Default/black_holes.h
@@ -206,7 +206,8 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
*/
__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
struct bpart* restrict bp, const struct black_holes_props* props,
- const struct phys_const* constants, const struct cosmology* cosmo) {}
+ const struct phys_const* constants, const struct cosmology* cosmo,
+ const double dt) {}
/**
* @brief Reset acceleration fields of a particle
diff --git a/src/black_holes/Default/black_holes_iact.h b/src/black_holes/Default/black_holes_iact.h
index bd9542f00c..67f36ede20 100644
--- a/src/black_holes/Default/black_holes_iact.h
+++ b/src/black_holes/Default/black_holes_iact.h
@@ -40,7 +40,8 @@ runner_iact_nonsym_bh_gas_density(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *bi, const struct part *pj, const struct xpart *xpj,
const int with_cosmology, const struct cosmology *cosmo,
- const struct gravity_props *grav_props, const integertime_t ti_current,
+ const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props, const integertime_t ti_current,
const double time) {
float wi, wi_dx;
@@ -94,6 +95,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
struct xpart *xpj, const int with_cosmology,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
const integertime_t ti_current,
const double time) {}
@@ -119,6 +121,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const struct bpart *bi, struct bpart *bj,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
const integertime_t ti_current) {}
/**
@@ -144,6 +147,7 @@ runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
struct xpart *xpj, const int with_cosmology,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
const integertime_t ti_current,
const double time) {
#ifdef DEBUG_INTERACTIONS_BH
diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index a731361211..e84bbd2991 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -31,6 +31,7 @@
/* Standard includes */
#include <float.h>
+#include <math.h>
/**
* @brief Computes the gravity time-step of a given black hole particle.
@@ -56,20 +57,35 @@ __attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
struct bpart* bp, const struct black_holes_props* props) {
bp->time_bin = 0;
- bp->subgrid_mass = bp->mass;
+ if (props->use_subgrid_mass_from_ics == 0)
+ bp->subgrid_mass = bp->mass;
+ else if (props->with_subgrid_mass_check && bp->subgrid_mass <= 0)
+ error(
+ "Black hole %lld has a subgrid mass of %f (internal units).\n"
+ "If this is because the ICs do not contain a 'SubgridMass' data "
+ "set, you should set the parameter "
+ "'EAGLEAGN:use_subgrid_mass_from_ics' to 0 to initialize the "
+ "black hole subgrid masses to the corresponding dynamical masses.\n"
+ "If the subgrid mass is intentionally set to this value, you can "
+ "disable this error by setting 'EAGLEAGN:with_subgrid_mass_check' "
+ "to 0.",
+ bp->id, bp->subgrid_mass);
bp->total_accreted_mass = 0.f;
bp->accretion_rate = 0.f;
bp->formation_time = -1.f;
bp->cumulative_number_seeds = 1;
bp->number_of_mergers = 0;
bp->number_of_gas_swallows = 0;
+ bp->number_of_direct_gas_swallows = 0;
+ bp->number_of_repositions = 0;
+ bp->number_of_reposition_attempts = 0;
+ bp->number_of_time_steps = 0;
bp->last_high_Eddington_fraction_scale_factor = -1.f;
bp->last_minor_merger_time = -1.;
bp->last_major_merger_time = -1.;
bp->swallowed_angular_momentum[0] = 0.f;
bp->swallowed_angular_momentum[1] = 0.f;
bp->swallowed_angular_momentum[2] = 0.f;
- bp->number_of_repositions = 0;
black_holes_mark_bpart_as_not_swallowed(&bp->merger_data);
}
@@ -105,6 +121,8 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
bp->reposition.delta_x[2] = -FLT_MAX;
bp->reposition.min_potential = FLT_MAX;
bp->reposition.potential = FLT_MAX;
+ bp->accretion_rate = 0.f; /* Optionally accumulated ngb-by-ngb */
+ bp->f_visc = FLT_MAX;
}
/**
@@ -196,28 +214,24 @@ __attribute__((always_inline)) INLINE static void black_holes_end_density(
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
- /* Finish the calculation by inserting the missing h-factors */
+ /* --- Finish the calculation by inserting the missing h factors --- */
bp->density.wcount *= h_inv_dim;
bp->density.wcount_dh *= h_inv_dim_plus_one;
bp->rho_gas *= h_inv_dim;
- bp->sound_speed_gas *= h_inv_dim;
- bp->velocity_gas[0] *= h_inv_dim;
- bp->velocity_gas[1] *= h_inv_dim;
- bp->velocity_gas[2] *= h_inv_dim;
- bp->circular_velocity_gas[0] *= h_inv_dim;
- bp->circular_velocity_gas[1] *= h_inv_dim;
- bp->circular_velocity_gas[2] *= h_inv_dim;
-
const float rho_inv = 1.f / bp->rho_gas;
- /* Finish the calculation by undoing the mass smoothing */
- bp->sound_speed_gas *= rho_inv;
- bp->velocity_gas[0] *= rho_inv;
- bp->velocity_gas[1] *= rho_inv;
- bp->velocity_gas[2] *= rho_inv;
- bp->circular_velocity_gas[0] *= rho_inv * h_inv;
- bp->circular_velocity_gas[1] *= rho_inv * h_inv;
- bp->circular_velocity_gas[2] *= rho_inv * h_inv;
+ /* For the following, we also have to undo the mass smoothing
+ * (N.B.: bp->velocity_gas is in BH frame, in internal units). */
+ bp->sound_speed_gas *= h_inv_dim * rho_inv;
+ bp->velocity_gas[0] *= h_inv_dim * rho_inv;
+ bp->velocity_gas[1] *= h_inv_dim * rho_inv;
+ bp->velocity_gas[2] *= h_inv_dim * rho_inv;
+
+ /* ... and for the circular velocity, convert from specifc angular
+ * momentum to circular velocity at the smoothing radius (extra h_inv) */
+ bp->circular_velocity_gas[0] *= h_inv_dim_plus_one * rho_inv;
+ bp->circular_velocity_gas[1] *= h_inv_dim_plus_one * rho_inv;
+ bp->circular_velocity_gas[2] *= h_inv_dim_plus_one * rho_inv;
}
/**
@@ -299,8 +313,10 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_part(
const float dr = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
message(
"BH %lld swallowing gas particle %lld "
- "(Delta_v = [%f, %f, %f] U_V, Delta_v_rad = %f)",
- bp->id, p->id, dv[0], dv[1], dv[2],
+ "(Delta_v = [%f, %f, %f] U_V, "
+ "Delta_x = [%f, %f, %f] U_L, "
+ "Delta_v_rad = %f)",
+ bp->id, p->id, -dv[0], -dv[1], -dv[2], -dx[0], -dx[1], -dx[2],
(dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]) / dr);
/* Update the BH metal masses */
@@ -310,6 +326,7 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_part(
/* This BH swallowed a gas particle */
bp->number_of_gas_swallows++;
+ bp->number_of_direct_gas_swallows++;
/* This BH lost a neighbour */
bp->num_ngbs--;
@@ -409,6 +426,9 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const struct phys_const* constants, const struct cosmology* cosmo,
const double time, const int with_cosmology, const double dt) {
+ /* Record that the black hole has another active time step */
+ bp->number_of_time_steps++;
+
if (dt == 0. || bp->rho_gas == 0.) return;
/* Gather some physical constants (all in internal units) */
@@ -426,37 +446,20 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const double delta_T = props->AGN_delta_T_desired;
const double delta_u = delta_T * props->temp_to_u_factor;
const double alpha_visc = props->alpha_visc;
+ const int with_angmom_limiter = props->with_angmom_limiter;
/* (Subgrid) mass of the BH (internal units) */
const double BH_mass = bp->subgrid_mass;
- /* Convert the quantities we gathered to physical frame (all internal units)
+ /* Convert the quantities we gathered to physical frame (all internal units).
* Note: for the velocities this means peculiar velocities */
- const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
const double gas_c_phys = bp->sound_speed_gas * cosmo->a_factor_sound_speed;
- const double gas_v_peculiar[3] = {bp->velocity_gas[0] * cosmo->a_inv,
- bp->velocity_gas[1] * cosmo->a_inv,
- bp->velocity_gas[2] * cosmo->a_inv};
-
- const double bh_v_peculiar[3] = {bp->v[0] * cosmo->a_inv,
- bp->v[1] * cosmo->a_inv,
- bp->v[2] * cosmo->a_inv};
-
+ const double gas_c_phys2 = gas_c_phys * gas_c_phys;
const double gas_v_circular[3] = {
bp->circular_velocity_gas[0] * cosmo->a_inv,
bp->circular_velocity_gas[1] * cosmo->a_inv,
bp->circular_velocity_gas[2] * cosmo->a_inv};
- /* Difference in peculiar velocity between the gas and the BH
- * Note that there is no need for a Hubble flow term here. We are
- * computing the gas velocity at the position of the black hole. */
- const double v_diff_peculiar[3] = {gas_v_peculiar[0] - bh_v_peculiar[0],
- gas_v_peculiar[1] - bh_v_peculiar[1],
- gas_v_peculiar[2] - bh_v_peculiar[2]};
- const double v_diff_norm2 = v_diff_peculiar[0] * v_diff_peculiar[0] +
- v_diff_peculiar[1] * v_diff_peculiar[1] +
- v_diff_peculiar[2] * v_diff_peculiar[2];
-
/* Norm of the circular velocity of the gas around the BH */
const double tangential_velocity2 = gas_v_circular[0] * gas_v_circular[0] +
gas_v_circular[1] * gas_v_circular[1] +
@@ -464,14 +467,48 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const double tangential_velocity = sqrt(tangential_velocity2);
/* We can now compute the Bondi accretion rate (internal units) */
- const double gas_c_phys2 = gas_c_phys * gas_c_phys;
- const double denominator2 = v_diff_norm2 + gas_c_phys2;
- const double denominator_inv = 1. / sqrt(denominator2);
- double Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass * gas_rho_phys *
- denominator_inv * denominator_inv * denominator_inv;
+ double Bondi_rate;
+
+ if (props->multi_phase_bondi) {
+
+ /* In this case, we are in 'multi-phase-Bondi' mode -- otherwise,
+ * the accretion_rate is still zero (was initialised to this) */
+ const float hi_inv = 1.f / bp->h;
+ const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
+
+ Bondi_rate = bp->accretion_rate *
+ (4. * M_PI * G * G * BH_mass * BH_mass * hi_inv_dim);
+ } else {
+
+ /* Standard approach: compute accretion rate for all gas simultaneously.
+ *
+ * Convert the quantities we gathered to physical frame (all internal
+ * units). Note: velocities are already in black hole frame. */
+ const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
+ const double gas_v_phys[3] = {bp->velocity_gas[0] * cosmo->a_inv,
+ bp->velocity_gas[1] * cosmo->a_inv,
+ bp->velocity_gas[2] * cosmo->a_inv};
+
+ const double gas_v_norm2 = gas_v_phys[0] * gas_v_phys[0] +
+ gas_v_phys[1] * gas_v_phys[1] +
+ gas_v_phys[2] * gas_v_phys[2];
+
+ const double denominator2 = gas_v_norm2 + gas_c_phys2;
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure that the denominator is strictly positive */
+ if (denominator2 <= 0)
+ error(
+ "Invalid denominator for black hole particle %lld in Bondi rate "
+ "calculation.",
+ bp->id);
+#endif
+ const double denominator_inv = 1. / sqrt(denominator2);
+ Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass * gas_rho_phys *
+ denominator_inv * denominator_inv * denominator_inv;
+ }
/* Compute the reduction factor from Rosas-Guevara et al. (2015) */
- if (props->with_angmom_limiter) {
+ if (with_angmom_limiter) {
const double Bondi_radius = G * BH_mass / gas_c_phys2;
const double Bondi_time = Bondi_radius / gas_c_phys;
const double r_times_v_tang = Bondi_radius * tangential_velocity;
@@ -479,17 +516,21 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
r_times_v_tang * r_times_v_tang * r_times_v_tang;
const double viscous_time = 2. * M_PI * r_times_v_tang_3 /
(1e-6 * alpha_visc * G * G * BH_mass * BH_mass);
+
const double f_visc = min(Bondi_time / viscous_time, 1.);
+ bp->f_visc = f_visc;
- /* Limit the Bondi rate by the Bondi viscuous time ratio */
+ /* Limit the accretion rate by the Bondi-to-viscous time ratio */
Bondi_rate *= f_visc;
+ } else {
+ bp->f_visc = 1.0;
}
/* Compute the Eddington rate (internal units) */
const double Eddington_rate =
4. * M_PI * G * BH_mass * proton_mass / (epsilon_r * c * sigma_Thomson);
- /* Shall we record this high rate? */
+ /* Should we record this time as the most recent high accretion rate? */
if (Bondi_rate > f_Edd_recording * Eddington_rate) {
if (with_cosmology) {
bp->last_high_Eddington_fraction_scale_factor = cosmo->a;
@@ -498,7 +539,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
}
}
- /* Limit the accretion rate to the Eddington fraction */
+ /* Limit the accretion rate to a fraction of the Eddington rate */
const double accr_rate = min(Bondi_rate, f_Edd * Eddington_rate);
bp->accretion_rate = accr_rate;
@@ -567,24 +608,68 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
+ * @param dt The black hole particle's time step.
*/
__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
struct bpart* restrict bp, const struct black_holes_props* props,
- const struct phys_const* constants, const struct cosmology* cosmo) {
-
- const float potential = gravity_get_comoving_potential(bp->gpart);
+ const struct phys_const* constants, const struct cosmology* cosmo,
+ const double dt) {
- /* Is the potential lower (i.e. the BH is at the bottom already)
- * OR is the BH massive enough that we don't reposition? */
- if (potential < bp->reposition.min_potential ||
- bp->subgrid_mass > props->max_reposition_mass) {
+ /* First check: did we find any eligible neighbour particle to jump to? */
+ if (bp->reposition.min_potential != FLT_MAX) {
- /* No need to reposition */
- bp->reposition.min_potential = FLT_MAX;
- bp->reposition.delta_x[0] = -FLT_MAX;
- bp->reposition.delta_x[1] = -FLT_MAX;
- bp->reposition.delta_x[2] = -FLT_MAX;
- }
+ /* Record that we have a (possible) repositioning situation */
+ bp->number_of_reposition_attempts++;
+
+ /* Is the potential lower (i.e. the BH is at the bottom already)
+ * OR is the BH massive enough that we don't reposition? */
+ const float potential = gravity_get_comoving_potential(bp->gpart);
+ if (potential < bp->reposition.min_potential ||
+ bp->subgrid_mass > props->max_reposition_mass) {
+
+ /* No need to reposition */
+ bp->reposition.min_potential = FLT_MAX;
+ bp->reposition.delta_x[0] = -FLT_MAX;
+ bp->reposition.delta_x[1] = -FLT_MAX;
+ bp->reposition.delta_x[2] = -FLT_MAX;
+
+ } else if (props->set_reposition_speed) {
+
+ /* If we are re-positioning, move the BH a fraction of delta_x, so
+ * that we have a well-defined re-positioning velocity. We have
+ * checked already that reposition_coefficient_upsilon is positive. */
+ const float repos_vel =
+ props->reposition_coefficient_upsilon *
+ pow(bp->subgrid_mass / constants->const_solar_mass,
+ props->reposition_exponent_xi);
+
+ const double dx = bp->reposition.delta_x[0];
+ const double dy = bp->reposition.delta_x[1];
+ const double dz = bp->reposition.delta_x[2];
+ const double d = sqrt(dx * dx + dy * dy + dz * dz);
+
+ /* Convert target reposition velocity to a fractional reposition
+ * along reposition.delta_x */
+
+ /* Exclude the pathological case of repositioning by zero distance */
+ if (d > 0) {
+ double repos_frac = repos_vel * dt / d;
+
+ /* We should never get negative repositioning fractions... */
+ if (repos_frac < 0)
+ error("Wanting to reposition by negative fraction (%g)?", repos_frac);
+
+ /* ... but fractions > 1 can occur if the target velocity is high.
+ * We do not want this, because it could lead to overshooting the
+ * actual potential minimum. */
+ if (repos_frac > 1) repos_frac = 1.;
+
+ bp->reposition.delta_x[0] *= repos_frac;
+ bp->reposition.delta_x[1] *= repos_frac;
+ bp->reposition.delta_x[2] *= repos_frac;
+ }
+ } /* ends section for fractional repositioning */
+ } /* ends section if we found eligible repositioning target(s) */
}
/**
@@ -668,9 +753,12 @@ INLINE static void black_holes_create_from_gas(
bp->cumulative_number_seeds = 1;
bp->number_of_mergers = 0;
bp->number_of_gas_swallows = 0;
+ bp->number_of_direct_gas_swallows = 0;
+ bp->number_of_time_steps = 0;
- /* We haven't repositioned yet */
+ /* We haven't repositioned yet, nor attempted it */
bp->number_of_repositions = 0;
+ bp->number_of_reposition_attempts = 0;
/* Initial metal masses */
const float gas_mass = hydro_get_mass(p);
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index b0626bb646..30186f892e 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -38,16 +38,20 @@
* @param bi First particle (black hole).
* @param pj Second particle (gas, not updated).
* @param xpj The extended data of the second particle (not updated).
+ * @param with_cosmology Are we doing a cosmological run?
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
+ * @param bh_props The properties of the BH scheme
* @param ti_current Current integer time value (for random numbers).
+ * @param time Current physical time in the simulation.
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_density(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *bi, const struct part *pj, const struct xpart *xpj,
const int with_cosmology, const struct cosmology *cosmo,
- const struct gravity_props *grav_props, const integertime_t ti_current,
+ const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props, const integertime_t ti_current,
const double time) {
float wi, wi_dx;
@@ -83,20 +87,56 @@ runner_iact_nonsym_bh_gas_density(
/* Contribution to the smoothed sound speed */
bi->sound_speed_gas += mj * cj * wi;
- /* Neighbour peculiar drifted velocity */
- const float vj[3] = {pj->v[0], pj->v[1], pj->v[2]};
+ /* Neighbour's (drifted) velocity in the frame of the black hole
+ * (we don't include a Hubble term since we are interested in the
+ * velocity contribution at the location of the black hole) */
+ const float dv[3] = {pj->v[0] - bi->v[0], pj->v[1] - bi->v[1],
+ pj->v[2] - bi->v[2]};
- /* Contribution to the smoothed velocity */
- bi->velocity_gas[0] += mj * vj[0] * wi;
- bi->velocity_gas[1] += mj * vj[1] * wi;
- bi->velocity_gas[2] += mj * vj[2] * wi;
+ /* Contribution to the smoothed velocity (gas w.r.t. black hole) */
+ bi->velocity_gas[0] += mj * dv[0] * wi;
+ bi->velocity_gas[1] += mj * dv[1] * wi;
+ bi->velocity_gas[2] += mj * dv[2] * wi;
- /* Contribution to the circular valocity */
- const float dv[3] = {bi->v[0] - vj[0], bi->v[1] - vj[1], bi->v[2] - vj[2]};
+ /* Contribution to the specific angular momentum of gas, which is later
+ * converted to the circular velocity at the smoothing length */
bi->circular_velocity_gas[0] += mj * wi * (dx[1] * dv[2] - dx[2] * dv[1]);
bi->circular_velocity_gas[1] += mj * wi * (dx[2] * dv[0] - dx[0] * dv[2]);
bi->circular_velocity_gas[2] += mj * wi * (dx[0] * dv[1] - dx[1] * dv[0]);
+ if (bh_props->multi_phase_bondi) {
+ /* Contribution to BH accretion rate
+ *
+ * i) Calculate denominator in Bondi formula */
+ const double gas_v_phys[3] = {dv[0] * cosmo->a_inv, dv[1] * cosmo->a_inv,
+ dv[2] * cosmo->a_inv};
+ const double gas_v_norm2 = gas_v_phys[0] * gas_v_phys[0] +
+ gas_v_phys[1] * gas_v_phys[1] +
+ gas_v_phys[2] * gas_v_phys[2];
+
+ const double gas_c_phys = cj * cosmo->a_factor_sound_speed;
+ const double gas_c_phys2 = gas_c_phys * gas_c_phys;
+ const double denominator2 = gas_v_norm2 + gas_c_phys2;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure that the denominator is strictly positive */
+ if (denominator2 <= 0)
+ error(
+ "Invalid denominator for BH particle %lld and gas particle "
+ "%lld in Bondi rate calculation.",
+ bi->id, pj->id);
+#endif
+ const double denominator_inv = 1. / sqrt(denominator2);
+
+ /* ii) Contribution of gas particle to the BH accretion rate
+ * (without constant pre-factor)
+ * N.B.: rhoj is the weighted contribution to BH gas density. */
+ const float rhoj = mj * wi * cosmo->a3_inv;
+ bi->accretion_rate +=
+ rhoj * denominator_inv * denominator_inv * denominator_inv;
+
+ } /* End of accretion contribution calculation */
+
#ifdef DEBUG_INTERACTIONS_BH
/* Update ngb counters */
if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_BH)
@@ -120,9 +160,12 @@ runner_iact_nonsym_bh_gas_density(
* @param bi First particle (black hole).
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
+ * @param with_cosmology Are we doing a cosmological run?
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
+ * @param bh_props The properties of the BH scheme
* @param ti_current Current integer time value (for random numbers).
+ * @param time Current physical time in the simulation.
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
@@ -131,6 +174,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
struct xpart *xpj, const int with_cosmology,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
const integertime_t ti_current,
const double time) {
@@ -152,26 +196,46 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
const float max_dist_repos2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv *
- const_max_repositioning_distance_ratio *
- const_max_repositioning_distance_ratio * grav_props->epsilon_baryon_cur *
+ bh_props->max_reposition_distance_ratio *
+ bh_props->max_reposition_distance_ratio * grav_props->epsilon_baryon_cur *
grav_props->epsilon_baryon_cur;
- /* This gas neighbour is close enough that we can consider it's potential
- for repositioning */
+ /* Is this gas neighbour close enough that we can consider its potential
+ for repositioning? */
if (r2 < max_dist_repos2) {
- /* Compute relative peculiar velocity between the two BHs
- * Recall that in SWIFT v is (v_pec * a) */
- const float delta_v[3] = {bi->v[0] - pj->v[0], bi->v[1] - pj->v[1],
- bi->v[2] - pj->v[2]};
- const float v2 = delta_v[0] * delta_v[0] + delta_v[1] * delta_v[1] +
- delta_v[2] * delta_v[2];
-
- const float v2_pec = v2 * cosmo->a2_inv;
+ /* Flag to check whether neighbour is slow enough to be considered
+ * as repositioning target. Always true if velocity cut is switched off. */
+ int neighbour_is_slow_enough = 1;
+ if (bh_props->with_reposition_velocity_threshold) {
+
+ /* Compute relative peculiar velocity between the two BHs
+ * Recall that in SWIFT v is (v_pec * a) */
+ const float delta_v[3] = {bi->v[0] - pj->v[0], bi->v[1] - pj->v[1],
+ bi->v[2] - pj->v[2]};
+ const float v2 = delta_v[0] * delta_v[0] + delta_v[1] * delta_v[1] +
+ delta_v[2] * delta_v[2];
+ const float v2_pec = v2 * cosmo->a2_inv;
+
+ /* Compute the maximum allowed velocity */
+ float v2_max = bh_props->max_reposition_velocity_ratio *
+ bh_props->max_reposition_velocity_ratio *
+ bi->sound_speed_gas * bi->sound_speed_gas;
+
+ /* If desired, limit the value of the threshold (v2_max) to be no
+ * smaller than a user-defined value */
+ if (bh_props->min_reposition_velocity_threshold > 0) {
+ const float v2_min_thresh =
+ bh_props->min_reposition_velocity_threshold *
+ bh_props->min_reposition_velocity_threshold;
+ v2_max = max(v2_max, v2_min_thresh);
+ }
- /* Check the velocity criterion */
- if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
+ /* Is the neighbour too fast to jump to? */
+ if (v2_pec >= v2_max) neighbour_is_slow_enough = 0;
+ }
+ if (neighbour_is_slow_enough) {
const float potential = pj->black_holes_data.potential;
/* Is the potential lower? */
@@ -235,6 +299,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
* @param bj Second particle (black hole)
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
+ * @param bh_props The properties of the BH scheme
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
@@ -243,6 +308,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
struct bpart *bi, struct bpart *bj,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
const integertime_t ti_current) {
/* Compute relative peculiar velocity between the two BHs
@@ -258,17 +324,38 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float max_dist_repos2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv *
- const_max_repositioning_distance_ratio *
- const_max_repositioning_distance_ratio * grav_props->epsilon_baryon_cur *
+ bh_props->max_reposition_distance_ratio *
+ bh_props->max_reposition_distance_ratio * grav_props->epsilon_baryon_cur *
grav_props->epsilon_baryon_cur;
- /* This BH neighbour is close enough that we can consider it's potential
- for repositioning */
+ /* Is this BH neighbour close enough that we can consider its potential
+ for repositioning? */
if (r2 < max_dist_repos2) {
- /* Check the velocity criterion */
- if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
+ /* Flag to check whether neighbour is slow enough to be considered
+ * as repositioning target. Always true if velocity cut switched off */
+ int neighbour_is_slow_enough = 1;
+ if (bh_props->with_reposition_velocity_threshold) {
+
+ /* Compute the maximum allowed velocity */
+ float v2_max = bh_props->max_reposition_velocity_ratio *
+ bh_props->max_reposition_velocity_ratio *
+ bi->sound_speed_gas * bi->sound_speed_gas;
+
+ /* If desired, limit the value of the threshold (v2_max) to be no
+ * smaller than a user-defined value */
+ if (bh_props->min_reposition_velocity_threshold > 0) {
+ const float v2_min_thresh =
+ bh_props->min_reposition_velocity_threshold *
+ bh_props->min_reposition_velocity_threshold;
+ v2_max = max(v2_max, v2_min_thresh);
+ }
+
+ /* Is the neighbour too fast to jump to? */
+ if (v2_pec >= v2_max) neighbour_is_slow_enough = 0;
+ }
+ if (neighbour_is_slow_enough) {
const float potential = bj->reposition.potential;
/* Is the potential lower? */
@@ -285,16 +372,19 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
/* Find the most massive of the two BHs */
float M = bi->subgrid_mass;
+ float h = hi;
if (bj->subgrid_mass > M) {
M = bj->subgrid_mass;
+ h = hj;
}
/* (Square of) max swallowing distance allowed based on the softening */
const float max_dist_merge2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv *
- const_max_merging_distance_ratio * const_max_merging_distance_ratio *
- grav_props->epsilon_baryon_cur * grav_props->epsilon_baryon_cur;
+ bh_props->max_merging_distance_ratio *
+ bh_props->max_merging_distance_ratio * grav_props->epsilon_baryon_cur *
+ grav_props->epsilon_baryon_cur;
const float G_Newton = grav_props->G_Newton;
@@ -305,10 +395,40 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
if ((bj->subgrid_mass < bi->subgrid_mass) ||
(bj->subgrid_mass == bi->subgrid_mass && bj->id < bi->id)) {
+ /* Merge if gravitationally bound AND if within max distance
+ * Note that we use the kernel support here as the size and not just the
+ * smoothing length */
+
/* Maximum velocity difference between BHs allowed to merge */
- const float v2_threshold = 2.f * G_Newton * M / sqrt(r2);
+ float v2_threshold;
+
+ if (bh_props->merger_threshold_type == 0) {
+
+ /* 'Old-style' merger threshold using circular velocity at the
+ * edge of the more massive BH's kernel */
+ v2_threshold = G_Newton * M / (kernel_gamma * h);
+ } else {
+
+ /* Arguably better merger threshold using the escape velocity at
+ * the distance of the lower-mass BH */
+ const float r_12 = sqrt(r2);
+
+ if ((bh_props->merger_threshold_type == 1) &&
+ (r_12 < grav_props->epsilon_baryon_cur)) {
+
+ /* If BHs are within softening range, take this into account */
+ float w_grav;
+ kernel_grav_pot_eval(r_12 / grav_props->epsilon_baryon_cur, &w_grav);
+ const float r_mod = w_grav / grav_props->epsilon_baryon_cur;
+ v2_threshold = 2.f * G_Newton * M / (r_mod);
+
+ } else {
+
+ /* Standard formula if BH interactions are not softened */
+ v2_threshold = 2.f * G_Newton * M / (r_12);
+ }
+ } /* Ends sections for different merger thresholds */
- /* Merge if gravitationally bound AND if within max distance */
if ((v2_pec < v2_threshold) && (r2 < max_dist_merge2)) {
/* This particle is swallowed by the BH with the largest ID of all the
@@ -346,6 +466,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
* @param with_cosmology Are we doing a cosmological run?
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
+ * @param bh_props The properties of the BH scheme
* @param ti_current Current integer time value (for random numbers).
* @param time current physical time in the simulation
*/
@@ -356,6 +477,7 @@ runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
struct xpart *xpj, const int with_cosmology,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
const integertime_t ti_current,
const double time) {
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index 9bd5ac199e..31fed833aa 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -35,7 +35,7 @@ INLINE static void black_holes_read_particles(struct bpart* bparts,
int* num_fields) {
/* Say how much we want to read */
- *num_fields = 6;
+ *num_fields = 7;
/* List what we want to read */
list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
@@ -50,6 +50,8 @@ INLINE static void black_holes_read_particles(struct bpart* bparts,
UNIT_CONV_LENGTH, bparts, h);
list[5] = io_make_input_field("EnergyReservoir", FLOAT, 1, OPTIONAL,
UNIT_CONV_ENERGY, bparts, energy_reservoir);
+ list[6] = io_make_input_field("SubgridMasses", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_MASS, bparts, subgrid_mass);
}
INLINE static void convert_bpart_pos(const struct engine* e,
@@ -94,7 +96,7 @@ INLINE static void convert_bpart_vel(const struct engine* e,
ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
- /* Conversion from internal units to peculiar velocities */
+ /* Conversion from internal to physical units */
ret[0] *= cosmo->a_inv;
ret[1] *= cosmo->a_inv;
ret[2] *= cosmo->a_inv;
@@ -105,18 +107,10 @@ INLINE static void convert_bpart_gas_vel(const struct engine* e,
const struct cosmology* cosmo = e->cosmology;
- /* Convert velocities to peculiar velocities */
- const double gas_v_peculiar[3] = {bp->velocity_gas[0] * cosmo->a_inv,
- bp->velocity_gas[1] * cosmo->a_inv,
- bp->velocity_gas[2] * cosmo->a_inv};
-
- const double bh_v_peculiar[3] = {bp->v[0] * cosmo->a_inv,
- bp->v[1] * cosmo->a_inv,
- bp->v[2] * cosmo->a_inv};
-
- ret[0] = gas_v_peculiar[0] - bh_v_peculiar[0];
- ret[1] = gas_v_peculiar[1] - bh_v_peculiar[1];
- ret[2] = gas_v_peculiar[2] - bh_v_peculiar[2];
+ /* Convert relative velocities to physical units */
+ ret[0] = bp->velocity_gas[0] * cosmo->a_inv;
+ ret[1] = bp->velocity_gas[1] * cosmo->a_inv;
+ ret[2] = bp->velocity_gas[2] * cosmo->a_inv;
}
INLINE static void convert_bpart_gas_circular_vel(const struct engine* e,
@@ -125,7 +119,7 @@ INLINE static void convert_bpart_gas_circular_vel(const struct engine* e,
const struct cosmology* cosmo = e->cosmology;
- /* Conversion from internal units to peculiar velocities */
+ /* Conversion from internal to physical units */
ret[0] = bp->circular_velocity_gas[0] * cosmo->a_inv;
ret[1] = bp->circular_velocity_gas[1] * cosmo->a_inv;
ret[2] = bp->circular_velocity_gas[2] * cosmo->a_inv;
@@ -145,7 +139,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
int with_cosmology) {
/* Say how much we want to write */
- *num_fields = 22;
+ *num_fields = 27;
/* List what we want to write */
list[0] = io_make_output_field_convert_bpart(
@@ -188,8 +182,8 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
"Co-moving densities of the gas around the particles");
list[8] = io_make_output_field(
- "GasSoundSpeeds", FLOAT, 1, UNIT_CONV_SPEED, 1.5f * hydro_gamma_minus_one,
- bparts, sound_speed_gas,
+ "GasSoundSpeeds", FLOAT, 1, UNIT_CONV_SPEED,
+ -1.5f * hydro_gamma_minus_one, bparts, sound_speed_gas,
"Co-moving sound-speeds of the gas around the particles");
list[9] = io_make_output_field(
@@ -208,7 +202,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
"Total mass accreted onto the particles since its birth");
list[12] = io_make_output_field(
- "CumulativeNumberSeeds", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ "CumulativeNumberOfSeeds", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
cumulative_number_seeds,
"Total number of BH seeds that have merged into this black hole");
@@ -273,23 +267,56 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
list[19] = io_make_output_field_convert_bpart(
"GasCircularVelocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, bparts,
convert_bpart_gas_circular_vel,
- "Peculiar circular velocities of the gas particles around the black "
- "holes. This is the curl of a * dx/dt where x is the co-moving position "
- "of the particles.");
+ "Circular velocities of the gas around the black hole at the "
+ "smoothing radius. This is j / h_BH, where j is the smoothed, peculiar "
+ "specific angular momentum of gas around the black holes, and h_BH is "
+ "the smoothing length of each black hole.");
- list[20] = io_make_output_field(
+ list[20] =
+ io_make_output_field("TimeBins", CHAR, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ time_bin, "Time-bins of the particles");
+
+ list[21] = io_make_output_field(
"NumberOfSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
number_of_gas_swallows,
"Number of gas particles the black holes have swallowed. "
- "This includes the particles swallowed by any of the BHs that "
+ "This includes the particles swallowed by any of the black holes that "
"merged into this one.");
- list[21] = io_make_output_field(
+ list[22] = io_make_output_field(
+ "NumberOfDirectSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ number_of_direct_gas_swallows,
+ "Number of gas particles the black holes have swallowed. "
+ "This does not include any particles swallowed by any of the black holes "
+ "that merged into this one.");
+
+ list[23] = io_make_output_field(
"NumberOfRepositions", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
number_of_repositions,
"Number of repositioning events the black holes went through. This does "
"not include the number of reposition events accumulated by any merged "
- "black hole.");
+ "black holes.");
+
+ list[24] = io_make_output_field(
+ "NumberOfRepositionAttempts", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ number_of_reposition_attempts,
+ "Number of time steps in which the black holes had an eligible particle "
+ "to reposition to. They may or may not have ended up moving there, "
+ "depending on their subgrid mass and on whether these particles were at "
+ "a lower or higher potential than the black holes themselves. It does "
+ "not include attempted repositioning events accumulated by any merged "
+ "black holes.");
+
+ list[25] = io_make_output_field(
+ "NumberOfTimeSteps", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ number_of_time_steps,
+ "Total number of time steps at which the black holes were active.");
+
+ list[26] = io_make_output_field(
+ "ViscosityFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts, f_visc,
+ "Multiplicative factors by which the Bondi-Hoyle-Lyttleton accretion "
+ "rates have been suppressed by the Rosas-Guevara et al. (2015) "
+ "accretion disc model.");
#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
diff --git a/src/black_holes/EAGLE/black_holes_parameters.h b/src/black_holes/EAGLE/black_holes_parameters.h
index cb6438bc85..c34ce6eb9e 100644
--- a/src/black_holes/EAGLE/black_holes_parameters.h
+++ b/src/black_holes/EAGLE/black_holes_parameters.h
@@ -26,6 +26,8 @@
* @file EAGLE/black_holes_parameters.h
* @brief Parameters of the EAGLE black holes
* model that need to be defined at compile time.
+ *
+ * @note In this branch, these properties are not used anywhere!
*/
/*! Maximal distance for merging particles in units of the (spline not Plummer)
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index 0d19657ebc..ed2f3bca72 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -93,12 +93,18 @@ struct bpart {
/*! Smoothed sound speed of the gas surrounding the black hole. */
float sound_speed_gas;
- /*! Smoothed velocity (peculiar) of the gas surrounding the black hole */
+ /*! Smoothed velocity of the gas surrounding the black hole,
+ * in the frame of the black hole (internal units) */
float velocity_gas[3];
- /*! Curl of the velocity field around the black hole */
+ /*! Circular velocity of the gas around the black hole at the smoothing
+ * radius (calculated as j_gas / h_BH, where j is specific ang. mom.) */
float circular_velocity_gas[3];
+ /*! Multiplicative factor for accretion rates, from Rosas-Guevara et al.
+ * (2015) angular momentum based accretion disc model */
+ float f_visc;
+
/*! Total mass of the gas neighbours. */
float ngb_mass;
@@ -111,16 +117,29 @@ struct bpart {
/*! Total number of BH merger events (i.e. not including all progenies) */
int number_of_mergers;
- /*! Total number of gas particles that have been swallowed */
+ /*! Total number of gas particles swallowed (including particles swallowed
+ * by merged-in black holes) */
int number_of_gas_swallows;
- /*! Total (physical) angular momentum accumulated by swallowing particles */
- float swallowed_angular_momentum[3];
+ /*! Total number of gas particles swallowed (excluding particles swallowed
+ * by merged-in black holes) */
+ int number_of_direct_gas_swallows;
- /*! Total number of times this BH was repositioned (not including progenies)
- */
+ /*! Total number of times the black hole has been repositioned (excluding
+ * repositionings of merged-in black holes) */
int number_of_repositions;
+ /*! Total number of times a black hole attempted repositioning (including
+ * cases where it was aborted because the black hole was already at a
+ * lower potential than all eligible neighbours) */
+ int number_of_reposition_attempts;
+
+ /*! Total number of time steps in which the black hole was active. */
+ int number_of_time_steps;
+
+ /*! Total (physical) angular momentum accumulated by swallowing particles */
+ float swallowed_angular_momentum[3];
+
/*! Union for the last high Eddington ratio point in time */
union {
diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h
index 54068b3484..fba2df4ba6 100644
--- a/src/black_holes/EAGLE/black_holes_properties.h
+++ b/src/black_holes/EAGLE/black_holes_properties.h
@@ -52,28 +52,38 @@ struct black_holes_props {
/*! Mass of a BH seed at creation time */
float subgrid_seed_mass;
- /* ----- Properties of the accretion model ------ */
+ /*! Should we use the subgrid mass specified in ICs? */
+ int use_subgrid_mass_from_ics;
- /*! Maximal fraction of the Eddington rate allowed. */
- float f_Edd;
+ /*! Should we enforce positive subgrid masses initially? */
+ int with_subgrid_mass_check;
- /*! Radiative efficiency of the black holes. */
- float epsilon_r;
+ /* ----- Properties of the accretion model ------ */
- /*! Feedback coupling efficiency of the black holes. */
- float epsilon_f;
+ /*! Calculate Bondi accretion rate for individual neighbours? */
+ int multi_phase_bondi;
- /*! Are we using the Rosas-Guevara et al. (2015) term? */
+ /*! Are we applying the angular-momentum-based multiplicative term from
+ * Rosas-Guevara et al. (2015)? */
int with_angmom_limiter;
/*! Normalisation of the viscuous angular momentum accretion reduction */
float alpha_visc;
+ /*! Radiative efficiency of the black holes. */
+ float epsilon_r;
+
+ /*! Maximal fraction of the Eddington rate allowed. */
+ float f_Edd;
+
/*! Eddington fraction threshold for recording */
float f_Edd_recording;
/* ---- Properties of the feedback model ------- */
+ /*! Feedback coupling efficiency of the black holes. */
+ float epsilon_f;
+
/*! Temperature increase induced by AGN feedback (Kelvin) */
float AGN_delta_T_desired;
@@ -85,6 +95,29 @@ struct black_holes_props {
/*! Maximal mass of BH to reposition */
float max_reposition_mass;
+ /*! Maximal distance to reposition, in units of softening length */
+ float max_reposition_distance_ratio;
+
+ /*! Switch to enable a relative velocity limit for particles to which the
+ * black holes can reposition */
+ int with_reposition_velocity_threshold;
+
+ /*! Maximal velocity offset of particles to which the black hole can
+ * reposition, in units of the ambient sound speed of the black hole */
+ float max_reposition_velocity_ratio;
+
+ /*! Minimum value of the velocity repositioning threshold */
+ float min_reposition_velocity_threshold;
+
+ /*! Switch to enable repositioning at fixed (maximum) speed */
+ int set_reposition_speed;
+
+ /*! Normalisation factor for repositioning velocity */
+ float reposition_coefficient_upsilon;
+
+ /*! Repositioning velocity scaling with black hole mass */
+ float reposition_exponent_xi;
+
/* ---- Properties of the merger model ---------- */
/*! Mass ratio above which a merger is considered 'minor' */
@@ -93,6 +126,12 @@ struct black_holes_props {
/*! Mass ratio above which a merger is considered 'major' */
float major_merger_threshold;
+ /*! Type of merger threshold (0: standard, 1: improved) */
+ int merger_threshold_type;
+
+ /*! Maximal distance over which BHs merge, in units of softening length */
+ float max_merging_distance_ratio;
+
/* ---- Common conversion factors --------------- */
/*! Conversion factor from temperature to internal energy */
@@ -159,15 +198,16 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
/* Convert to internal units */
bp->subgrid_seed_mass *= phys_const->const_solar_mass;
+ bp->use_subgrid_mass_from_ics =
+ parser_get_opt_param_int(params, "EAGLEAGN:use_subgrid_mass_from_ics", 1);
+ if (bp->use_subgrid_mass_from_ics)
+ bp->with_subgrid_mass_check =
+ parser_get_opt_param_int(params, "EAGLEAGN:with_subgrid_mass_check", 1);
+
/* Accretion parameters ---------------------------------- */
- bp->f_Edd = parser_get_param_float(params, "EAGLEAGN:max_eddington_fraction");
- bp->f_Edd_recording = parser_get_param_float(
- params, "EAGLEAGN:eddington_fraction_for_recording");
- bp->epsilon_r =
- parser_get_param_float(params, "EAGLEAGN:radiative_efficiency");
- bp->epsilon_f =
- parser_get_param_float(params, "EAGLEAGN:coupling_efficiency");
+ bp->multi_phase_bondi =
+ parser_get_param_int(params, "EAGLEAGN:multi_phase_bondi");
/* Rosas-Guevara et al. (2015) model */
bp->with_angmom_limiter =
@@ -175,8 +215,18 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
if (bp->with_angmom_limiter)
bp->alpha_visc = parser_get_param_float(params, "EAGLEAGN:viscous_alpha");
+ bp->epsilon_r =
+ parser_get_param_float(params, "EAGLEAGN:radiative_efficiency");
+
+ bp->f_Edd = parser_get_param_float(params, "EAGLEAGN:max_eddington_fraction");
+ bp->f_Edd_recording = parser_get_param_float(
+ params, "EAGLEAGN:eddington_fraction_for_recording");
+
/* Feedback parameters ---------------------------------- */
+ bp->epsilon_f =
+ parser_get_param_float(params, "EAGLEAGN:coupling_efficiency");
+
bp->AGN_delta_T_desired =
parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_K");
@@ -187,9 +237,51 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->max_reposition_mass =
parser_get_param_float(params, "EAGLEAGN:max_reposition_mass");
-
/* Convert to internal units */
bp->max_reposition_mass *= phys_const->const_solar_mass;
+ bp->max_reposition_distance_ratio =
+ parser_get_param_float(params, "EAGLEAGN:max_reposition_distance_ratio");
+
+ bp->with_reposition_velocity_threshold = parser_get_param_int(
+ params, "EAGLEAGN:with_reposition_velocity_threshold");
+
+ if (bp->with_reposition_velocity_threshold) {
+ bp->max_reposition_velocity_ratio = parser_get_param_float(
+ params, "EAGLEAGN:max_reposition_velocity_ratio");
+
+ /* Prevent nonsensical input */
+ if (bp->max_reposition_velocity_ratio <= 0)
+ error("max_reposition_velocity_ratio must be positive, not %f.",
+ bp->max_reposition_velocity_ratio);
+
+ bp->min_reposition_velocity_threshold = parser_get_param_float(
+ params, "EAGLEAGN:min_reposition_velocity_threshold");
+ /* Convert from km/s to internal units */
+ bp->min_reposition_velocity_threshold *=
+ (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs));
+ }
+
+ bp->set_reposition_speed =
+ parser_get_param_int(params, "EAGLEAGN:set_reposition_speed");
+
+ if (bp->set_reposition_speed) {
+ bp->reposition_coefficient_upsilon = parser_get_param_float(
+ params, "EAGLEAGN:reposition_coefficient_upsilon");
+
+ /* Prevent the user from making silly wishes */
+ if (bp->reposition_coefficient_upsilon <= 0)
+ error(
+ "reposition_coefficient_upsilon must be positive, not %f "
+ "km/s/M_sun.",
+ bp->reposition_coefficient_upsilon);
+
+ /* Convert from km/s to internal units */
+ bp->reposition_coefficient_upsilon *=
+ (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs));
+
+ bp->reposition_exponent_xi = parser_get_opt_param_float(
+ params, "EAGLEAGN:reposition_exponent_xi", 1.0);
+ }
/* Merger parameters ------------------------------------- */
@@ -199,6 +291,12 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->major_merger_threshold =
parser_get_param_float(params, "EAGLEAGN:threshold_major_merger");
+ bp->merger_threshold_type =
+ parser_get_param_int(params, "EAGLEAGN:merger_threshold_type");
+
+ bp->max_merging_distance_ratio =
+ parser_get_param_float(params, "EAGLEAGN:merger_max_distance_ratio");
+
/* Common conversion factors ----------------------------- */
/* Calculate temperature to internal energy conversion factor (all internal
diff --git a/src/cooling/EAGLE/cooling_io.h b/src/cooling/EAGLE/cooling_io.h
index e642185882..d6ab4f2254 100644
--- a/src/cooling/EAGLE/cooling_io.h
+++ b/src/cooling/EAGLE/cooling_io.h
@@ -24,6 +24,7 @@
/* Local includes */
#include "cooling.h"
+#include "engine.h"
#include "io_properties.h"
#ifdef HAVE_HDF5
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index 69424fa331..92efcafa4e 100644
--- a/src/cooling/const_du/cooling.h
+++ b/src/cooling/const_du/cooling.h
@@ -41,6 +41,7 @@
#include <math.h>
/* Local includes. */
+#include "cosmology.h"
#include "entropy_floor.h"
#include "hydro.h"
#include "parser.h"
diff --git a/src/cooling/const_du/cooling_io.h b/src/cooling/const_du/cooling_io.h
index 0ecb2ebe51..4fd7bc06a5 100644
--- a/src/cooling/const_du/cooling_io.h
+++ b/src/cooling/const_du/cooling_io.h
@@ -35,6 +35,7 @@
/* Local includes */
#include "cooling.h"
+#include "engine.h"
#include "io_properties.h"
#ifdef HAVE_HDF5
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 1597687bf6..b8f28f5999 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -37,6 +37,7 @@
/* Local includes. */
#include "const.h"
+#include "cosmology.h"
#include "entropy_floor.h"
#include "error.h"
#include "hydro.h"
diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h
index 091d567066..d68543ae6c 100644
--- a/src/cooling/const_lambda/cooling_io.h
+++ b/src/cooling/const_lambda/cooling_io.h
@@ -33,6 +33,7 @@
/* Local includes */
#include "cooling.h"
+#include "engine.h"
#include "io_properties.h"
#ifdef HAVE_HDF5
diff --git a/src/cooling/none/cooling_io.h b/src/cooling/none/cooling_io.h
index 848a90dd5e..d24db8f99d 100644
--- a/src/cooling/none/cooling_io.h
+++ b/src/cooling/none/cooling_io.h
@@ -24,6 +24,7 @@
/* Local includes */
#include "cooling.h"
+#include "engine.h"
#include "io_properties.h"
#ifdef HAVE_HDF5
diff --git a/src/distributed_io.c b/src/distributed_io.c
index 8a5a21f4fc..355fdb05ac 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -391,6 +391,7 @@ void write_output_distributed(struct engine* e,
swift_type_count);
io_write_attribute_i(h_grp, "NumFilesPerSnapshot", numFiles);
io_write_attribute_i(h_grp, "ThisFile", mpi_rank);
+ io_write_attribute_s(h_grp, "OutputType", "Snapshot");
/* Close header */
H5Gclose(h_grp);
diff --git a/src/engine.c b/src/engine.c
index 3ac99c047c..4fc82e59c9 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -70,6 +70,7 @@
#include "gravity.h"
#include "gravity_cache.h"
#include "hydro.h"
+#include "line_of_sight.h"
#include "logger.h"
#include "logger_io.h"
#include "map.h"
@@ -125,7 +126,8 @@ const char *engine_policy_names[] = {"none",
"fof search",
"time-step limiter",
"time-step sync",
- "logger"};
+ "logger",
+ "line of sight"};
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
@@ -2761,9 +2763,9 @@ void engine_step(struct engine *e) {
* @param e The #engine.
*/
void engine_check_for_dumps(struct engine *e) {
-
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int with_stf = (e->policy & engine_policy_structure_finding);
+ const int with_los = (e->policy & engine_policy_line_of_sight);
/* What kind of output are we getting? */
enum output_type {
@@ -2771,6 +2773,7 @@ void engine_check_for_dumps(struct engine *e) {
output_snapshot,
output_statistics,
output_stf,
+ output_los,
};
/* What kind of output do we want? And at which time ?
@@ -2806,6 +2809,16 @@ void engine_check_for_dumps(struct engine *e) {
}
}
+ /* Do we want to write a line of sight file? */
+ if (with_los) {
+ if (e->ti_end_min > e->ti_next_los && e->ti_next_los > 0) {
+ if (e->ti_next_los < ti_output) {
+ ti_output = e->ti_next_los;
+ type = output_los;
+ }
+ }
+ }
+
/* Store information before attempting extra dump-related drifts */
const integertime_t ti_current = e->ti_current;
const timebin_t max_active_bin = e->max_active_bin;
@@ -2892,6 +2905,16 @@ void engine_check_for_dumps(struct engine *e) {
#endif
break;
+ case output_los:
+
+ /* Compute the LoS */
+ do_line_of_sight(e);
+
+ /* Move on */
+ engine_compute_next_los_time(e);
+
+ break;
+
default:
error("Invalid dump type");
}
@@ -2928,6 +2951,16 @@ void engine_check_for_dumps(struct engine *e) {
}
}
+ /* Do line of sight ? */
+ if (with_los) {
+ if (e->ti_end_min > e->ti_next_los && e->ti_next_los > 0) {
+ if (e->ti_next_los < ti_output) {
+ ti_output = e->ti_next_los;
+ type = output_los;
+ }
+ }
+ }
+
} /* While loop over output types */
/* Restore the information we stored */
@@ -3819,7 +3852,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
const struct chemistry_global_data *chemistry,
- struct fof_props *fof_properties) {
+ struct fof_props *fof_properties,
+ struct los_props *los_properties) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
@@ -3868,6 +3902,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
units_init_default(e->snapshot_units, params, "Snapshots", internal_units);
e->snapshot_output_count = 0;
e->stf_output_count = 0;
+ e->los_output_count = 0;
e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
e->dt_max_RMS_displacement = FLT_MAX;
@@ -3900,6 +3935,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->fof_properties = fof_properties;
e->parameter_file = params;
e->stf_this_timestep = 0;
+ e->los_properties = los_properties;
#ifdef WITH_MPI
e->usertime_last_step = 0.0;
e->systime_last_step = 0.0;
@@ -3975,6 +4011,16 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
parser_get_opt_param_double(params, "StructureFinding:delta_time", -1.);
}
+ /* Initialise line of sight output. */
+ if (e->policy & engine_policy_line_of_sight) {
+ e->time_first_los =
+ parser_get_opt_param_double(params, "LineOfSight:time_first", 0.);
+ e->a_first_los = parser_get_opt_param_double(
+ params, "LineOfSight:scale_factor_first", 0.1);
+ e->delta_time_los =
+ parser_get_opt_param_double(params, "LineOfSight:delta_time", -1.);
+ }
+
/* Initialise FoF calls frequency. */
if (e->policy & engine_policy_fof) {
@@ -4456,6 +4502,11 @@ void engine_config(int restart, int fof, struct engine *e,
/* Find the time of the first statistics output */
engine_compute_next_statistics_time(e);
+ /* Find the time of the first line of sight output */
+ if (e->policy & engine_policy_line_of_sight) {
+ engine_compute_next_los_time(e);
+ }
+
/* Find the time of the first stf output */
if (e->policy & engine_policy_structure_finding) {
engine_compute_next_stf_time(e);
@@ -4868,6 +4919,77 @@ void engine_compute_next_statistics_time(struct engine *e) {
}
}
+/**
+ * @brief Computes the next time (on the time line) for a line of sight dump
+ *
+ * @param e The #engine.
+ */
+void engine_compute_next_los_time(struct engine *e) {
+ /* Do output_list file case */
+ if (e->output_list_los) {
+ output_list_read_next_time(e->output_list_los, e, "line of sights",
+ &e->ti_next_los);
+ return;
+ }
+
+ /* Find upper-bound on last output */
+ double time_end;
+ if (e->policy & engine_policy_cosmology)
+ time_end = e->cosmology->a_end * e->delta_time_los;
+ else
+ time_end = e->time_end + e->delta_time_los;
+
+ /* Find next los above current time */
+ double time;
+ if (e->policy & engine_policy_cosmology)
+ time = e->a_first_los;
+ else
+ time = e->time_first_los;
+
+ int found_los_time = 0;
+ while (time < time_end) {
+
+ /* Output time on the integer timeline */
+ if (e->policy & engine_policy_cosmology)
+ e->ti_next_los = log(time / e->cosmology->a_begin) / e->time_base;
+ else
+ e->ti_next_los = (time - e->time_begin) / e->time_base;
+
+ /* Found it? */
+ if (e->ti_next_los > e->ti_current) {
+ found_los_time = 1;
+ break;
+ }
+
+ if (e->policy & engine_policy_cosmology)
+ time *= e->delta_time_los;
+ else
+ time += e->delta_time_los;
+ }
+
+ /* Deal with last line of sight */
+ if (!found_los_time) {
+ e->ti_next_los = -1;
+ if (e->verbose) message("No further LOS output time.");
+ } else {
+
+ /* Be nice, talk... */
+ if (e->policy & engine_policy_cosmology) {
+ const double next_los_time =
+ exp(e->ti_next_los * e->time_base) * e->cosmology->a_begin;
+ if (e->verbose)
+ message("Next output time for line of sight set to a=%e.",
+ next_los_time);
+ } else {
+ const double next_los_time =
+ e->ti_next_los * e->time_base + e->time_begin;
+ if (e->verbose)
+ message("Next output time for line of sight set to t=%e.",
+ next_los_time);
+ }
+ }
+}
+
/**
* @brief Computes the next time (on the time line) for structure finding
*
@@ -5041,6 +5163,19 @@ void engine_init_output_lists(struct engine *e, struct swift_params *params) {
else
e->time_first_stf_output = stf_time_first;
}
+
+ /* Deal with line of sight */
+ double los_time_first;
+ e->output_list_los = NULL;
+ output_list_init(&e->output_list_los, e, "LineOfSight", &e->delta_time_los,
+ &los_time_first);
+
+ if (e->output_list_los) {
+ if (e->policy & engine_policy_cosmology)
+ e->a_first_los = los_time_first;
+ else
+ e->time_first_los = los_time_first;
+ }
}
/**
@@ -5256,12 +5391,14 @@ void engine_clean(struct engine *e, const int fof, const int restart) {
#ifdef WITH_FOF
free((void *)e->fof_properties);
#endif
+ free((void *)e->los_properties);
#ifdef WITH_MPI
free((void *)e->reparttype);
#endif
if (e->output_list_snapshots) free((void *)e->output_list_snapshots);
if (e->output_list_stats) free((void *)e->output_list_stats);
if (e->output_list_stf) free((void *)e->output_list_stf);
+ if (e->output_list_los) free((void *)e->output_list_los);
#ifdef WITH_LOGGER
if (e->policy & engine_policy_logger) free((void *)e->logger);
#endif
@@ -5310,12 +5447,14 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
#ifdef WITH_FOF
fof_struct_dump(e->fof_properties, stream);
#endif
+ los_struct_dump(e->los_properties, stream);
parser_struct_dump(e->parameter_file, stream);
if (e->output_list_snapshots)
output_list_struct_dump(e->output_list_snapshots, stream);
if (e->output_list_stats)
output_list_struct_dump(e->output_list_stats, stream);
if (e->output_list_stf) output_list_struct_dump(e->output_list_stf, stream);
+ if (e->output_list_los) output_list_struct_dump(e->output_list_los, stream);
#ifdef WITH_LOGGER
if (e->policy & engine_policy_logger) {
@@ -5441,6 +5580,11 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
e->fof_properties = fof_props;
#endif
+ struct los_props *los_properties =
+ (struct los_props *)malloc(sizeof(struct los_props));
+ los_struct_restore(los_properties, stream);
+ e->los_properties = los_properties;
+
struct swift_params *parameter_file =
(struct swift_params *)malloc(sizeof(struct swift_params));
parser_struct_restore(parameter_file, stream);
@@ -5467,6 +5611,13 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
e->output_list_stf = output_list_stf;
}
+ if (e->output_list_los) {
+ struct output_list *output_list_los =
+ (struct output_list *)malloc(sizeof(struct output_list));
+ output_list_struct_restore(output_list_los, stream);
+ e->output_list_los = output_list_los;
+ }
+
#ifdef WITH_LOGGER
if (e->policy & engine_policy_logger) {
struct logger_writer *log =
diff --git a/src/engine.h b/src/engine.h
index 3543667a56..e1ace452ec 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -79,8 +79,9 @@ enum engine_policy {
engine_policy_timestep_limiter = (1 << 21),
engine_policy_timestep_sync = (1 << 22),
engine_policy_logger = (1 << 23),
+ engine_policy_line_of_sight = (1 << 24),
};
-#define engine_maxpolicy 24
+#define engine_maxpolicy 25
extern const char *engine_policy_names[engine_maxpolicy + 1];
/**
@@ -488,6 +489,17 @@ struct engine {
/* Has there been an stf this timestep? */
char stf_this_timestep;
+ /* Line of sight properties. */
+ struct los_props *los_properties;
+
+ /* Line of sight outputs information. */
+ struct output_list *output_list_los;
+ double a_first_los;
+ double time_first_los;
+ double delta_time_los;
+ integertime_t ti_next_los;
+ int los_output_count;
+
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Run brute force checks only on steps when all gparts active? */
int force_checks_only_all_active;
@@ -510,6 +522,7 @@ void engine_compute_next_snapshot_time(struct engine *e);
void engine_compute_next_stf_time(struct engine *e);
void engine_compute_next_fof_time(struct engine *e);
void engine_compute_next_statistics_time(struct engine *e);
+void engine_compute_next_los_time(struct engine *e);
void engine_recompute_displacement_constraint(struct engine *e);
void engine_unskip(struct engine *e);
void engine_unskip_timestep_communications(struct engine *e);
@@ -538,7 +551,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
const struct chemistry_global_data *chemistry,
- struct fof_props *fof_properties);
+ struct fof_props *fof_properties,
+ struct los_props *los_properties);
void engine_config(int restart, int fof, struct engine *e,
struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int verbose,
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
index 329563d4e6..281ad02715 100644
--- a/src/hydro/AnarchyPU/hydro.h
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -887,23 +887,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
- internal_energy_ratio *= p->u;
-
- /* Now we can use this to 'update' the value of the smoothed pressure. To
- * truly update this variable, we would need another loop over neighbours
- * using the new internal energies of everyone, but that's not feasible. */
- p->pressure_bar *= internal_energy_ratio;
-
- /* Check against entropy floor */
- const float floor_A = entropy_floor(p, cosmo, floor_props);
- const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
-
- /* Check against absolute minimum */
- const float min_u =
- hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
-
- p->u = max(p->u, floor_u);
- p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
@@ -927,6 +910,26 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->pressure_bar *= expf_exact;
}
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, floor_u);
+ p->u = max(p->u, min_u);
+
+ /* Now that p->u has been properly bounded, we can use it to apply the
+ * drift for the pressure */
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
+
/* Compute the new sound speed */
const float soundspeed =
gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index c64ffada10..065a38e196 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -825,25 +825,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
- /* Check against entropy floor */
- const float floor_A = entropy_floor(p, cosmo, floor_props);
- const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
-
- /* Check against absolute minimum */
- const float min_u =
- hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
-
- p->u = max(p->u, floor_u);
- p->u = max(p->u, min_u);
-
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt_drift;
- if (fabsf(w1) < 0.2f)
+ if (fabsf(w1) < 0.2f) {
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
- else
+ } else {
p->h *= expf(w1);
+ }
/* Predict density and weighted pressure */
const float w2 = -hydro_dimension * w1;
@@ -856,6 +846,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf_exact;
}
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, floor_u);
+ p->u = max(p->u, min_u);
+
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 2172dae966..d195c68b15 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -718,6 +718,24 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the entropy */
p->entropy += p->entropy_dt * dt_therm;
+ const float h_inv = 1.f / p->h;
+
+ /* Predict smoothing length */
+ const float w1 = p->force.h_dt * h_inv * dt_drift;
+ if (fabsf(w1) < 0.2f) {
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ } else {
+ p->h *= expf(w1);
+ }
+
+ /* Predict density */
+ const float w2 = -hydro_dimension * w1;
+ if (fabsf(w2) < 0.2f) {
+ p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+ } else {
+ p->rho *= expf(w2);
+ }
+
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
@@ -732,22 +750,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->entropy = max(p->entropy, floor_A);
p->entropy = max(p->entropy, min_A);
- const float h_inv = 1.f / p->h;
-
- /* Predict smoothing length */
- const float w1 = p->force.h_dt * h_inv * dt_drift;
- if (fabsf(w1) < 0.2f)
- p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
- else
- p->h *= expf(w1);
-
- /* Predict density */
- const float w2 = -hydro_dimension * w1;
- if (fabsf(w2) < 0.2f)
- p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
- else
- p->rho *= expf(w2);
-
/* Re-compute the pressure */
float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
comoving_pressure =
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index e49575e769..36eb28ccf4 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -705,32 +705,34 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
- /* Check against entropy floor */
- const float floor_A = entropy_floor(p, cosmo, floor_props);
- const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
-
- /* Check against absolute minimum */
- const float min_u =
- hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
-
- p->u = max(p->u, floor_u);
- p->u = max(p->u, min_u);
-
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt_drift;
- if (fabsf(w1) < 0.2f)
+ if (fabsf(w1) < 0.2f) {
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
- else
+ } else {
p->h *= expf(w1);
+ }
/* Predict density */
const float w2 = -hydro_dimension * w1;
- if (fabsf(w2) < 0.2f)
+ if (fabsf(w2) < 0.2f) {
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
- else
+ } else {
p->rho *= expf(w2);
+ }
+
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, floor_u);
+ p->u = max(p->u, min_u);
/* Compute the new pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index fddafa52a9..6f68f92d1a 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -750,23 +750,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
- internal_energy_ratio *= p->u;
-
- /* Now we can use this to 'update' the value of the smoothed pressure. To
- * truly update this variable, we would need another loop over neighbours
- * using the new internal energies of everyone, but that's not feasible. */
- p->pressure_bar *= internal_energy_ratio;
-
- /* Check against entropy floor */
- const float floor_A = entropy_floor(p, cosmo, floor_props);
- const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
-
- /* Check against absolute minimum */
- const float min_u =
- hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
-
- p->u = max(p->u, floor_u);
- p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
@@ -790,6 +773,26 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->pressure_bar *= expf_exact;
}
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, floor_u);
+ p->u = max(p->u, min_u);
+
+ /* Now that p->u has been properly bounded, we can use it to apply the
+ * drift for the pressure */
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
+
/* Compute the new sound speed */
hydro_update_soundspeed(p, cosmo);
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index 27508dd7ea..11b4659086 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -747,23 +747,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
- internal_energy_ratio *= p->u;
-
- /* Now we can use this to 'update' the value of the smoothed pressure. To
- * truly update this variable, we would need another loop over neighbours
- * using the new internal energies of everyone, but that's not feasible. */
- p->pressure_bar *= internal_energy_ratio;
-
- /* Check against entropy floor */
- const float floor_A = entropy_floor(p, cosmo, floor_props);
- const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
-
- /* Check against absolute minimum */
- const float min_u =
- hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
-
- p->u = max(p->u, floor_u);
- p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
@@ -787,6 +770,26 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->pressure_bar *= expf_exact;
}
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, floor_u);
+ p->u = max(p->u, min_u);
+
+ /* Now that p->u has been properly bounded, we can use it to apply the
+ * drift for the pressure */
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
+
/* Compute the new sound speed */
const float soundspeed =
gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index fd19256ebf..f0504b382e 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -683,17 +683,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the entropy */
p->entropy += p->entropy_dt * dt_therm;
- /* Check against entropy floor */
- const float floor_A = entropy_floor(p, cosmo, floor_props);
-
- /* Check against absolute minimum */
- const float min_u = hydro_props->minimal_internal_energy;
- const float min_A =
- gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u);
-
- p->entropy = max(p->entropy, floor_A);
- p->entropy = max(p->entropy, min_A);
-
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
@@ -713,6 +702,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho_bar *= expf(w2);
}
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+
+ /* Check against absolute minimum */
+ const float min_u = hydro_props->minimal_internal_energy;
+ const float min_A =
+ gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u);
+
+ p->entropy = max(p->entropy, floor_A);
+ p->entropy = max(p->entropy, min_A);
+
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h
index b0da3ac888..ed158d7d10 100644
--- a/src/hydro/SPHENIX/hydro.h
+++ b/src/hydro/SPHENIX/hydro.h
@@ -746,26 +746,33 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Construct the source term for the AV; if shock detected this is _positive_
* as div_v_dt should be _negative_ before the shock hits */
- const float S = kernel_support_physical * kernel_support_physical *
- max(0.f, -1.f * div_v_dt);
- /* 0.25 factor comes from our definition of v_sig (sum of soundspeeds rather
- * than mean). */
- /* Note this is v_sig_physical squared, not comoving */
- const float v_sig_square = 0.25 * v_sig_physical * v_sig_physical;
+ /* Source term is only activated if flow is converging (i.e. in the pre-
+ * shock region) */
+ const float S = p->viscosity.div_v < 0.f
+ ? kernel_support_physical * kernel_support_physical *
+ max(0.f, -1.f * div_v_dt)
+ : 0.f;
+
+ /* We want the decay to occur based on the thermodynamic properties
+ * of the gas - hence use the soundspeed instead of the signal velocity */
+ const float soundspeed_square = soundspeed_physical * soundspeed_physical;
/* Calculate the current appropriate value of the AV based on the above */
const float alpha_loc =
- hydro_props->viscosity.alpha_max * S / (v_sig_square + S);
+ hydro_props->viscosity.alpha_max * S / (soundspeed_square + S);
if (alpha_loc > p->viscosity.alpha) {
/* Reset the value of alpha to the appropriate value */
p->viscosity.alpha = alpha_loc;
} else {
/* Integrate the alpha forward in time to decay back to alpha = alpha_loc */
- p->viscosity.alpha =
- alpha_loc + (p->viscosity.alpha - alpha_loc) *
- expf(-dt_alpha * sound_crossing_time_inverse *
- hydro_props->viscosity.length);
+ /* This type of integration is stable w.r.t. different time-step lengths
+ * (Price 2018) */
+ const float timescale_ratio =
+ dt_alpha * sound_crossing_time_inverse * hydro_props->viscosity.length;
+
+ p->viscosity.alpha += alpha_loc * timescale_ratio;
+ p->viscosity.alpha /= (1.f + timescale_ratio);
}
/* Check that we did not hit the minimum */
@@ -901,17 +908,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
- /* Check against entropy floor */
- const float floor_A = entropy_floor(p, cosmo, floor_props);
- const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
-
- /* Check against absolute minimum */
- const float min_u =
- hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
-
- p->u = max(p->u, floor_u);
- p->u = max(p->u, min_u);
-
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
@@ -932,6 +928,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf_exact;
}
+ /* Check against entropy floor - explicitly do this after drifting the
+ * density as this has a density dependence. */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, floor_u);
+ p->u = max(p->u, min_u);
+
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
const float pressure_including_floor =
diff --git a/src/hydro/SPHENIX/hydro_parameters.h b/src/hydro/SPHENIX/hydro_parameters.h
index d77693a50e..06b77e9eec 100644
--- a/src/hydro/SPHENIX/hydro_parameters.h
+++ b/src/hydro/SPHENIX/hydro_parameters.h
@@ -71,7 +71,7 @@
#define hydro_props_default_viscosity_alpha_max 2.0f
/*! Decay length for the viscosity scheme. This is scheme dependent. */
-#define hydro_props_default_viscosity_length 0.25f
+#define hydro_props_default_viscosity_length 0.05f
/* Diffusion parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
diff --git a/src/line_of_sight.c b/src/line_of_sight.c
new file mode 100644
index 0000000000..cbd5426fb0
--- /dev/null
+++ b/src/line_of_sight.c
@@ -0,0 +1,1006 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Stuart McAlpine (stuart.mcalpine@helsinki.fi)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* MPI headers. */
+#ifdef WITH_MPI
+#include <mpi.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "atomic.h"
+#include "chemistry_io.h"
+#include "cooling_io.h"
+#include "engine.h"
+#include "fof_io.h"
+#include "hydro_io.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+#include "line_of_sight.h"
+#include "periodic.h"
+#include "star_formation_io.h"
+#include "tracers_io.h"
+#include "velociraptor_io.h"
+
+/**
+ * @brief Will the line of sight intersect a given cell?
+ *
+ * Also return 0 if the cell is empty.
+ *
+ * @param c The top level cell.
+ * @param los The line of sight structure.
+ */
+static INLINE int does_los_intersect(const struct cell *c,
+ const struct line_of_sight *los) {
+
+ /* Empty cell? */
+ if (c->hydro.count == 0) return 0;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->hydro.h_max <= 0.) error("Invalid h_max for does_los_intersect");
+#endif
+
+ /* Distance from LOS to left and bottom cell edge. */
+ const double cx_min = c->loc[los->xaxis];
+ const double cy_min = c->loc[los->yaxis];
+
+ /* Distance from LOS to right and top cell edge. */
+ const double cx_max = c->loc[los->xaxis] + c->width[los->xaxis];
+ const double cy_max = c->loc[los->yaxis] + c->width[los->yaxis];
+
+ /* Maximum smoothing length of a part in this top level cell. */
+ const double hsml = c->hydro.h_max * kernel_gamma;
+ const double hsml2 = hsml * hsml;
+
+ double dx, dy;
+
+ if (los->periodic) {
+ dx = min(fabs(nearest(cx_min - los->Xpos, los->dim[los->xaxis])),
+ fabs(nearest(cx_max - los->Xpos, los->dim[los->xaxis])));
+ dy = min(fabs(nearest(cy_min - los->Ypos, los->dim[los->yaxis])),
+ fabs(nearest(cy_max - los->Ypos, los->dim[los->yaxis])));
+ } else {
+ dx = fabs(cx_max - los->Xpos);
+ dy = fabs(cy_max - los->Ypos);
+ }
+
+ /* Is sightline directly within this top level cell? */
+ if (dx < (1.01 * c->width[los->xaxis]) / 2. &&
+ dy < (1.01 * c->width[los->yaxis]) / 2.) {
+ return 1;
+ /* Could a part from this top level cell smooth into the sightline? */
+ } else if (dx * dx + dy * dy < hsml2) {
+ return 1;
+ /* Don't need to work with this top level cell. */
+ } else {
+ return 0;
+ }
+}
+
+/**
+ * @brief Reads the LOS properties from the param file.
+ *
+ * @param dim Space dimensions.
+ * @param los_params Sightline parameters to save into.
+ * @param params Swift params to read from.
+ */
+void los_init(const double dim[3], struct los_props *los_params,
+ struct swift_params *params) {
+
+ /* How many line of sights in each plane. */
+ los_params->num_along_x =
+ parser_get_param_int(params, "LineOfSight:num_along_x");
+ los_params->num_along_y =
+ parser_get_param_int(params, "LineOfSight:num_along_y");
+ los_params->num_along_z =
+ parser_get_param_int(params, "LineOfSight:num_along_z");
+
+ /* Min/max range across x,y and z (simulation axes) where random
+ * LOS's are allowed. */
+ los_params->allowed_losrange_x[0] = 0.;
+ los_params->allowed_losrange_x[1] = dim[0];
+ parser_get_opt_param_double_array(params, "LineOfSight:allowed_los_range_x",
+ 2, los_params->allowed_losrange_x);
+ los_params->allowed_losrange_y[0] = 0.;
+ los_params->allowed_losrange_y[1] = dim[1];
+ parser_get_opt_param_double_array(params, "LineOfSight:allowed_los_range_y",
+ 2, los_params->allowed_losrange_y);
+ los_params->allowed_losrange_z[0] = 0.;
+ los_params->allowed_losrange_z[1] = dim[2];
+ parser_get_opt_param_double_array(params, "LineOfSight:allowed_los_range_z",
+ 2, los_params->allowed_losrange_z);
+
+ /* Compute total number of sightlines. */
+ los_params->num_tot = los_params->num_along_z + los_params->num_along_x +
+ los_params->num_along_y;
+
+ /* Where are we saving them? */
+ parser_get_param_string(params, "LineOfSight:basename", los_params->basename);
+
+ /* Min/max range allowed when sightline is shooting down x,y and z
+ * (simulation axes). */
+ los_params->range_when_shooting_down_axis[0][0] = 0.;
+ los_params->range_when_shooting_down_axis[0][1] = dim[0];
+ parser_get_opt_param_double_array(
+ params, "LineOfSight:range_when_shooting_down_x", 2,
+ los_params->range_when_shooting_down_axis[0]);
+ los_params->range_when_shooting_down_axis[1][0] = 0.;
+ los_params->range_when_shooting_down_axis[1][1] = dim[1];
+ parser_get_opt_param_double_array(
+ params, "LineOfSight:range_when_shooting_down_y", 2,
+ los_params->range_when_shooting_down_axis[1]);
+ los_params->range_when_shooting_down_axis[2][0] = 0.;
+ los_params->range_when_shooting_down_axis[2][1] = dim[2];
+ parser_get_opt_param_double_array(
+ params, "LineOfSight:range_when_shooting_down_z", 2,
+ los_params->range_when_shooting_down_axis[2]);
+}
+
+/**
+ * @brief Create a #line_of_sight object from its attributes
+ */
+void create_sightline(const double Xpos, const double Ypos,
+ enum los_direction xaxis, enum los_direction yaxis,
+ enum los_direction zaxis, const int periodic,
+ const double dim[3], struct line_of_sight *los,
+ const double range_when_shooting_down_axis[2]) {
+ los->Xpos = Xpos;
+ los->Ypos = Ypos;
+ los->particles_in_los_local = 0;
+ los->particles_in_los_total = 0;
+ los->xaxis = xaxis;
+ los->yaxis = yaxis;
+ los->zaxis = zaxis;
+ los->periodic = periodic;
+ los->dim[0] = dim[0];
+ los->dim[1] = dim[1];
+ los->dim[2] = dim[2];
+ los->num_intersecting_top_level_cells = 0;
+ los->range_when_shooting_down_axis[0] = range_when_shooting_down_axis[0];
+ los->range_when_shooting_down_axis[1] = range_when_shooting_down_axis[1];
+}
+
+/**
+ * @brief Generates random sightline positions.
+ *
+ * Independent sightlines are made for the XY, YZ and XZ planes.
+ *
+ * @param LOS Structure to store sightlines.
+ * @param params Sightline parameters.
+ */
+void generate_sightlines(struct line_of_sight *Los,
+ const struct los_props *params, const int periodic,
+ const double dim[3]) {
+
+ /* Keep track of number of sightlines. */
+ int count = 0;
+
+ /* Sightlines in XY plane, shoots down Z. */
+ for (int i = 0; i < params->num_along_z; i++) {
+ double Xpos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_x[1] - params->allowed_losrange_x[0])) +
+ params->allowed_losrange_x[0];
+ double Ypos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_y[1] - params->allowed_losrange_y[0])) +
+ params->allowed_losrange_y[0];
+
+ create_sightline(Xpos, Ypos, simulation_x_axis, simulation_y_axis,
+ simulation_z_axis, periodic, dim, &Los[count],
+ params->range_when_shooting_down_axis[simulation_z_axis]);
+ count++;
+ }
+
+ /* Sightlines in YZ plane, shoots down X. */
+ for (int i = 0; i < params->num_along_x; i++) {
+ double Xpos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_y[1] - params->allowed_losrange_y[0])) +
+ params->allowed_losrange_y[0];
+ double Ypos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_z[1] - params->allowed_losrange_z[0])) +
+ params->allowed_losrange_z[0];
+
+ create_sightline(Xpos, Ypos, simulation_y_axis, simulation_z_axis,
+ simulation_x_axis, periodic, dim, &Los[count],
+ params->range_when_shooting_down_axis[simulation_x_axis]);
+ count++;
+ }
+
+ /* Sightlines in XZ plane, shoots down Y. */
+ for (int i = 0; i < params->num_along_y; i++) {
+ double Xpos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_x[1] - params->allowed_losrange_x[0])) +
+ params->allowed_losrange_x[0];
+ double Ypos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_z[1] - params->allowed_losrange_z[0])) +
+ params->allowed_losrange_z[0];
+
+ create_sightline(Xpos, Ypos, simulation_x_axis, simulation_z_axis,
+ simulation_y_axis, periodic, dim, &Los[count],
+ params->range_when_shooting_down_axis[simulation_y_axis]);
+ count++;
+ }
+
+ /* Make sure we made the correct ammount */
+ if (count != params->num_tot)
+ error("Could not make the right number of sightlines");
+}
+
+/**
+ * @brief Print line_of_sight information.
+ *
+ * @param Los Structure to print.
+ */
+void print_los_info(const struct line_of_sight *Los, const int i) {
+
+ message(
+ "[LOS %i] Xpos:%g Ypos:%g parts_in_los:%i "
+ "num_intersecting_top_level_cells:%i",
+ i, Los[i].Xpos, Los[i].Ypos, Los[i].particles_in_los_total,
+ Los[i].num_intersecting_top_level_cells);
+}
+
+/**
+ * @brief Writes dataset for a given part attribute.
+ *
+ * @param p io_props dataset for this attribute.
+ * @param N number of parts in this line of sight.
+ * @param j Line of sight ID.
+ * @param e The engine.
+ * @param grp HDF5 group to write to.
+ */
+void write_los_hdf5_dataset(const struct io_props props, const size_t N,
+ const int j, const struct engine *e, hid_t grp) {
+
+ /* Create data space */
+ const hid_t h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0)
+ error("Error while creating data space for field '%s'.", props.name);
+
+ int rank = 0;
+ hsize_t shape[2];
+ hsize_t chunk_shape[2];
+ if (props.dimension > 1) {
+ rank = 2;
+ shape[0] = N;
+ shape[1] = props.dimension;
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
+ chunk_shape[1] = props.dimension;
+ } else {
+ rank = 1;
+ shape[0] = N;
+ shape[1] = 0;
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
+ chunk_shape[1] = 0;
+ }
+
+ /* Make sure the chunks are not larger than the dataset */
+ if (chunk_shape[0] > N) chunk_shape[0] = N;
+
+ /* Change shape of data space */
+ hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, shape);
+ if (h_err < 0)
+ error("Error while changing data space shape for field '%s'.", props.name);
+
+ /* Dataset properties */
+ const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
+
+ /* Set chunk size */
+ h_err = H5Pset_chunk(h_prop, rank, chunk_shape);
+ if (h_err < 0)
+ error("Error while setting chunk size (%llu, %llu) for field '%s'.",
+ chunk_shape[0], chunk_shape[1], props.name);
+
+ /* Impose check-sum to verify data corruption */
+ h_err = H5Pset_fletcher32(h_prop);
+ if (h_err < 0)
+ error("Error while setting checksum options for field '%s'.", props.name);
+
+ /* Impose data compression */
+ if (e->snapshot_compression > 0) {
+ h_err = H5Pset_shuffle(h_prop);
+ if (h_err < 0)
+ error("Error while setting shuffling options for field '%s'.",
+ props.name);
+
+ h_err = H5Pset_deflate(h_prop, e->snapshot_compression);
+ if (h_err < 0)
+ error("Error while setting compression options for field '%s'.",
+ props.name);
+ }
+
+ /* Allocate temporary buffer */
+ const size_t num_elements = N * props.dimension;
+ const size_t typeSize = io_sizeof_type(props.type);
+ void *temp = NULL;
+ if (swift_memalign("writebuff", (void **)&temp, IO_BUFFER_ALIGNMENT,
+ num_elements * typeSize) != 0)
+ error("Unable to allocate temporary i/o buffer");
+
+ /* Copy particle data to temp buffer */
+ io_copy_temp_buffer(temp, e, props, N, e->internal_units, e->snapshot_units);
+
+ /* Create dataset */
+ char att_name[200];
+ sprintf(att_name, "/LOS_%04i/%s", j, props.name);
+ const hid_t h_data = H5Dcreate(grp, att_name, io_hdf5_type(props.type),
+ h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace '%s'.", props.name);
+
+ /* Write dataset */
+ herr_t status = H5Dwrite(h_data, io_hdf5_type(props.type), H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, temp);
+ if (status < 0) error("Error while writing data array '%s'.", props.name);
+
+ /* Write unit conversion factors for this data set */
+ char buffer[FIELD_BUFFER_SIZE] = {0};
+ units_cgs_conversion_string(buffer, e->snapshot_units, props.units,
+ props.scale_factor_exponent);
+ float baseUnitsExp[5];
+ units_get_base_unit_exponents_array(baseUnitsExp, props.units);
+ io_write_attribute_f(h_data, "U_M exponent", baseUnitsExp[UNIT_MASS]);
+ io_write_attribute_f(h_data, "U_L exponent", baseUnitsExp[UNIT_LENGTH]);
+ io_write_attribute_f(h_data, "U_t exponent", baseUnitsExp[UNIT_TIME]);
+ io_write_attribute_f(h_data, "U_I exponent", baseUnitsExp[UNIT_CURRENT]);
+ io_write_attribute_f(h_data, "U_T exponent", baseUnitsExp[UNIT_TEMPERATURE]);
+ io_write_attribute_f(h_data, "h-scale exponent", 0.f);
+ io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
+ io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
+
+ /* Write the actual number this conversion factor corresponds to */
+ const double factor =
+ units_cgs_conversion_factor(e->snapshot_units, props.units);
+ io_write_attribute_d(
+ h_data,
+ "Conversion factor to CGS (not including cosmological corrections)",
+ factor);
+ io_write_attribute_d(
+ h_data,
+ "Conversion factor to physical CGS (including cosmological corrections)",
+ factor * pow(e->cosmology->a, props.scale_factor_exponent));
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (strlen(props.description) == 0)
+ error("Invalid (empty) description of the field '%s'", props.name);
+#endif
+
+ /* Write the full description */
+ io_write_attribute_s(h_data, "Description", props.description);
+
+ /* Free and close everything */
+ swift_free("writebuff", temp);
+ H5Pclose(h_prop);
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+}
+
+/**
+ * @brief Write parts in LOS to HDF5 file.
+ *
+ * @param grp HDF5 group of this LOS.
+ * @param j LOS ID.
+ * @param N number of parts in this line of sight.
+ * @param parts the list of parts in this LOS.
+ * @param e The engine.
+ * @param xparts the list of xparts in this LOS.
+ */
+void write_los_hdf5_datasets(hid_t grp, const int j, const size_t N,
+ const struct part *parts, const struct engine *e,
+ const struct xpart *xparts) {
+
+ /* What kind of run are we working with? */
+ struct swift_params *params = e->parameter_file;
+ const int with_cosmology = e->policy & engine_policy_cosmology;
+ const int with_cooling = e->policy & engine_policy_cooling;
+ const int with_temperature = e->policy & engine_policy_temperature;
+ const int with_fof = e->policy & engine_policy_fof;
+#ifdef HAVE_VELOCIRAPTOR
+ const int with_stf = (e->policy & engine_policy_structure_finding) &&
+ (e->s->gpart_group_data != NULL);
+#else
+ const int with_stf = 0;
+#endif
+
+ int num_fields = 0;
+ struct io_props list[100];
+
+ /* Find all the gas output fields */
+ hydro_write_particles(parts, xparts, list, &num_fields);
+ num_fields += chemistry_write_particles(parts, list + num_fields);
+ if (with_cooling || with_temperature) {
+ num_fields += cooling_write_particles(parts, xparts, list + num_fields,
+ e->cooling_func);
+ }
+ if (with_fof) {
+ num_fields += fof_write_parts(parts, xparts, list + num_fields);
+ }
+ if (with_stf) {
+ num_fields += velociraptor_write_parts(parts, xparts, list + num_fields);
+ }
+ num_fields +=
+ tracers_write_particles(parts, xparts, list + num_fields, with_cosmology);
+ num_fields +=
+ star_formation_write_particles(parts, xparts, list + num_fields);
+
+ /* Loop over each output field */
+ for (int i = 0; i < num_fields; i++) {
+
+ /* Did the user cancel this field? */
+ char field[PARSER_MAX_LINE_SIZE];
+ sprintf(field, "SelectOutputLOS:%.*s", FIELD_BUFFER_SIZE, list[i].name);
+ int should_write = parser_get_opt_param_int(params, field, 1);
+
+ /* Write (if selected) */
+ if (should_write) write_los_hdf5_dataset(list[i], N, j, e, grp);
+ }
+}
+
+/**
+ * @brief Writes HDF5 headers and information groups for this line of sight.
+ *
+ * @param h_file HDF5 file reference.
+ * @param e The engine.
+ * @param LOS_params The line of sight params.
+ */
+void write_hdf5_header(hid_t h_file, const struct engine *e,
+ const struct los_props *LOS_params,
+ const size_t total_num_parts_in_los) {
+
+ /* Open header to write simulation properties */
+ hid_t h_grp =
+ H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating file header\n");
+
+ /* Convert basic output information to snapshot units */
+ const double factor_time = units_conversion_factor(
+ e->internal_units, e->snapshot_units, UNIT_CONV_TIME);
+ const double factor_length = units_conversion_factor(
+ e->internal_units, e->snapshot_units, UNIT_CONV_LENGTH);
+ const double dblTime = e->time * factor_time;
+ const double dim[3] = {e->s->dim[0] * factor_length,
+ e->s->dim[1] * factor_length,
+ e->s->dim[2] * factor_length};
+
+ io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
+ io_write_attribute_d(h_grp, "Time", dblTime);
+ io_write_attribute_d(h_grp, "Dimension", (int)hydro_dimension);
+ io_write_attribute_d(h_grp, "Redshift", e->cosmology->z);
+ io_write_attribute_d(h_grp, "Scale-factor", e->cosmology->a);
+ io_write_attribute_s(h_grp, "Code", "SWIFT");
+ io_write_attribute_s(h_grp, "RunName", e->run_name);
+
+ /* Store the time at which the snapshot was written */
+ time_t tm = time(NULL);
+ struct tm *timeinfo = localtime(&tm);
+ char snapshot_date[64];
+ strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
+ io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+
+ /* GADGET-2 legacy values */
+ /* Number of particles of each type */
+ long long N_total[swift_type_count] = {0};
+ N_total[0] = total_num_parts_in_los;
+ unsigned int numParticles[swift_type_count] = {0};
+ unsigned int numParticlesHighWord[swift_type_count] = {0};
+ for (int ptype = 0; ptype < swift_type_count; ++ptype) {
+ numParticles[ptype] = (unsigned int)N_total[ptype];
+ numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
+ }
+ io_write_attribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total,
+ swift_type_count);
+ io_write_attribute(h_grp, "NumPart_Total", UINT, numParticles,
+ swift_type_count);
+ io_write_attribute(h_grp, "NumPart_Total_HighWord", UINT,
+ numParticlesHighWord, swift_type_count);
+ io_write_attribute_i(h_grp, "NumFilesPerSnapshot", 1);
+ io_write_attribute_i(h_grp, "ThisFile", 0);
+ io_write_attribute_s(h_grp, "OutputType", "LineOfSight");
+
+ /* Close group */
+ H5Gclose(h_grp);
+
+ io_write_meta_data(h_file, e, e->internal_units, e->snapshot_units);
+
+ /* Print the LOS properties */
+ h_grp = H5Gcreate(h_file, "/LineOfSightParameters", H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating LOS group");
+
+ /* Record this LOS attributes */
+ const int num_los_per_axis[3] = {LOS_params->num_along_x,
+ LOS_params->num_along_y,
+ LOS_params->num_along_z};
+ io_write_attribute(h_grp, "NumLineOfSight_PerAxis", INT, num_los_per_axis, 3);
+ io_write_attribute(h_grp, "NumLineOfSight_Total", INT, &LOS_params->num_tot,
+ 1);
+ io_write_attribute(h_grp, "AllowedLOSRangeX", DOUBLE,
+ LOS_params->allowed_losrange_x, 2);
+ io_write_attribute(h_grp, "AllowedLOSRangeY", DOUBLE,
+ LOS_params->allowed_losrange_y, 2);
+ io_write_attribute(h_grp, "AllowedLOSRangeZ", DOUBLE,
+ LOS_params->allowed_losrange_z, 2);
+ H5Gclose(h_grp);
+}
+
+/**
+ * @brief Loop over each part to see which ones intersect the LOS.
+ *
+ * @param map_data The parts.
+ * @param count The number of parts.
+ * @param extra_data The line_of_sight structure for this LOS.
+ */
+void los_first_loop_mapper(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct line_of_sight *LOS_list = (struct line_of_sight *)extra_data;
+ const struct part *parts = (struct part *)map_data;
+
+ size_t los_particle_count = 0;
+
+ /* Loop over each part to find those in LOS. */
+ for (int i = 0; i < count; i++) {
+
+ /* Don't consider inhibited parts. */
+ if (parts[i].time_bin == time_bin_inhibited) continue;
+
+ /* Don't consider part if outwith allowed z-range. */
+ if (parts[i].x[LOS_list->zaxis] <
+ LOS_list->range_when_shooting_down_axis[0] ||
+ parts[i].x[LOS_list->zaxis] >
+ LOS_list->range_when_shooting_down_axis[1])
+ continue;
+
+ /* Distance from this part to LOS along x dim. */
+ double dx = parts[i].x[LOS_list->xaxis] - LOS_list->Xpos;
+
+ /* Periodic wrap. */
+ if (LOS_list->periodic) dx = nearest(dx, LOS_list->dim[LOS_list->xaxis]);
+
+ /* Square. */
+ const double dx2 = dx * dx;
+
+ /* Smoothing length of this part. */
+ const double hsml = parts[i].h * kernel_gamma;
+ const double hsml2 = hsml * hsml;
+
+ /* Does this particle fall into our LOS? */
+ if (dx2 < hsml2) {
+
+ /* Distance from this part to LOS along y dim. */
+ double dy = parts[i].x[LOS_list->yaxis] - LOS_list->Ypos;
+
+ /* Periodic wrap. */
+ if (LOS_list->periodic) dy = nearest(dy, LOS_list->dim[LOS_list->yaxis]);
+
+ /* Square. */
+ const double dy2 = dy * dy;
+
+ /* Does this part still fall into our LOS? */
+ if (dy2 < hsml2) {
+
+ /* 2D distance to LOS. */
+ if (dx2 + dy2 <= hsml2) {
+
+ /* We've found one. */
+ los_particle_count++;
+ }
+ }
+ }
+ } /* End of loop over all parts */
+
+ atomic_add(&LOS_list->particles_in_los_local, los_particle_count);
+}
+
+/**
+ * @brief Find all top level cells that a LOS will intersect.
+ *
+ * This includes the top level cells the LOS directly passes through
+ * and the neighbouring top level cells whose parts could smooth into the LOS.
+ *
+ * @param e The engine.
+ * @param los The line_of_sight structure.
+ */
+void find_intersecting_top_level_cells(
+ const struct engine *e, struct line_of_sight *los, int *los_cells_top,
+ const struct cell *cells, const int *local_cells_with_particles,
+ const int nr_local_cells_with_particles) {
+
+ /* Keep track of how many top level cells we intersect. */
+ int num_intersecting_top_level_cells = 0;
+
+ /* Loop over each top level cell */
+ for (int n = 0; n < nr_local_cells_with_particles; n++) {
+
+ /* This top level cell. */
+ const struct cell *c = &cells[local_cells_with_particles[n]];
+
+ if (does_los_intersect(c, los)) {
+ num_intersecting_top_level_cells++;
+ los_cells_top[n] = 1;
+ }
+ }
+
+#ifdef WITH_MPI
+ if (MPI_Allreduce(MPI_IN_PLACE, &num_intersecting_top_level_cells, 1, MPI_INT,
+ MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
+ error("Failed to allreduce num_intersecting_top_level_cells.");
+#endif
+
+ /* Store how many top level cells this LOS intersects. */
+ los->num_intersecting_top_level_cells = num_intersecting_top_level_cells;
+}
+
+/**
+ * @brief Main work function for computing line of sights.
+ *
+ * 1) Construct N random line of sight positions.
+ * 2) Loop over each line of sight.
+ * - 2.1) Find which top level cells sightline intersects.
+ * - 2.2) Loop over each part in these top level cells to see which intersect
+ * sightline.
+ * - 2.3) Use this count to construct a LOS parts/xparts array.
+ * - 2.4) Loop over each part and extract those in sightline to new array.
+ * - 2.5) Save sightline parts to HDF5 file.
+ *
+ * @param e The engine.
+ */
+void do_line_of_sight(struct engine *e) {
+
+ /* Start counting. */
+ const ticks tic = getticks();
+
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+ const struct los_props *LOS_params = e->los_properties;
+ const int verbose = e->verbose;
+
+ /* Start by generating the random sightline positions. */
+ struct line_of_sight *LOS_list = (struct line_of_sight *)malloc(
+ LOS_params->num_tot * sizeof(struct line_of_sight));
+
+ if (e->nodeID == 0) {
+ generate_sightlines(LOS_list, LOS_params, periodic, dim);
+ if (verbose)
+ message("Generated %i random sightlines.", LOS_params->num_tot);
+ }
+
+#ifdef WITH_MPI
+ /* Share the list of LoS with all the MPI ranks */
+ MPI_Bcast(LOS_list, LOS_params->num_tot * sizeof(struct line_of_sight),
+ MPI_BYTE, 0, MPI_COMM_WORLD);
+#endif
+
+ /* Node 0 creates the HDF5 file. */
+ hid_t h_file = -1, h_grp = -1;
+ char fileName[256], groupName[200];
+
+ if (e->nodeID == 0) {
+ sprintf(fileName, "%s_%04i.hdf5", LOS_params->basename,
+ e->los_output_count);
+ if (verbose) message("Creating LOS file: %s", fileName);
+ h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_file < 0) error("Error while opening file '%s'.", fileName);
+ }
+
+#ifdef WITH_MPI
+ MPI_Barrier(MPI_COMM_WORLD);
+#endif
+
+ /* Keep track of the total number of parts in all sightlines. */
+ size_t total_num_parts_in_los = 0;
+
+ /* ------------------------------- */
+ /* Main loop over each random LOS. */
+ /* ------------------------------- */
+
+ /* Get list of local non-empty top level cells. */
+ const struct cell *cells = e->s->cells_top;
+ const int *local_cells_with_particles = e->s->local_cells_with_particles_top;
+ const int nr_local_cells_with_particles = s->nr_local_cells_with_particles;
+
+ /* Loop over each random LOS. */
+ for (int j = 0; j < LOS_params->num_tot; j++) {
+
+ /* Create empty top level cell list for this LOS */
+ int *los_cells_top = (int *)swift_malloc(
+ "tl_cells_los", nr_local_cells_with_particles * sizeof(int));
+ bzero(los_cells_top, nr_local_cells_with_particles * sizeof(int));
+
+ /* Find all top level cells this LOS will intersect. */
+ find_intersecting_top_level_cells(e, &LOS_list[j], los_cells_top, cells,
+ local_cells_with_particles,
+ nr_local_cells_with_particles);
+
+ /* Next count all the parts that intersect with this line of sight */
+ for (int n = 0; n < nr_local_cells_with_particles; n++) {
+ if (los_cells_top[n] == 1) {
+ const struct cell *c = &cells[local_cells_with_particles[n]];
+ threadpool_map(&s->e->threadpool, los_first_loop_mapper, c->hydro.parts,
+ c->hydro.count, sizeof(struct part),
+ threadpool_auto_chunk_size, &LOS_list[j]);
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Confirm we are capturing all the parts that intersect the LOS by redoing
+ * the count looping over all parts in the space (not just those in the
+ * subset of top level cells). */
+
+ struct part *parts = s->parts;
+ const size_t nr_parts = s->nr_parts;
+ const int old_particles_in_los_local = LOS_list[j].particles_in_los_local;
+ LOS_list[j].particles_in_los_local = 0;
+
+ /* Count all parts that intersect with this line of sight. */
+ threadpool_map(&s->e->threadpool, los_first_loop_mapper, parts, nr_parts,
+ sizeof(struct part), threadpool_auto_chunk_size,
+ &LOS_list[j]);
+
+ /* Make sure we get the same answer as above. */
+ if (old_particles_in_los_local != LOS_list[j].particles_in_los_local)
+ error("Space vs cells don't match s:%d != c:%d",
+ LOS_list[j].particles_in_los_local, old_particles_in_los_local);
+#endif
+
+#ifdef WITH_MPI
+ /* Make sure all nodes know how many parts are in this LOS */
+ int *counts = (int *)malloc(sizeof(int) * e->nr_nodes);
+ int *offsets = (int *)malloc(sizeof(int) * e->nr_nodes);
+
+ /* How many parts does each rank have for this LOS? */
+ MPI_Allgather(&LOS_list[j].particles_in_los_local, 1, MPI_INT, counts, 1,
+ MPI_INT, MPI_COMM_WORLD);
+
+ int offset_count = 0;
+ for (int k = 0; k < e->nr_nodes; k++) {
+
+ /* Total parts in this LOS. */
+ LOS_list[j].particles_in_los_total += counts[k];
+
+ /* Counts and offsets for Gatherv. */
+ offsets[k] = offset_count;
+ offset_count += counts[k];
+ }
+#else
+ LOS_list[j].particles_in_los_total = LOS_list[j].particles_in_los_local;
+#endif
+ total_num_parts_in_los += LOS_list[j].particles_in_los_total;
+
+ /* Print information about this LOS */
+ if (e->nodeID == 0) print_los_info(LOS_list, j);
+
+ /* Don't work with empty LOS */
+ if (LOS_list[j].particles_in_los_total == 0) {
+ if (e->nodeID == 0) {
+ message("*WARNING* LOS %i is empty", j);
+ print_los_info(LOS_list, j);
+ }
+#ifdef WITH_MPI
+ free(counts);
+ free(offsets);
+#endif
+ swift_free("tl_cells_los", los_cells_top);
+ continue;
+ }
+
+ /* Setup LOS part and xpart structures. */
+ struct part *LOS_parts = NULL;
+ struct xpart *LOS_xparts = NULL;
+
+ /* Rank 0 allocates more space as it will gather all the data for writing */
+ if (e->nodeID == 0) {
+ if ((LOS_parts = (struct part *)swift_malloc(
+ "los_parts_array",
+ sizeof(struct part) * LOS_list[j].particles_in_los_total)) ==
+ NULL)
+ error("Failed to allocate LOS part memory.");
+ if ((LOS_xparts = (struct xpart *)swift_malloc(
+ "los_xparts_array",
+ sizeof(struct xpart) * LOS_list[j].particles_in_los_total)) ==
+ NULL)
+ error("Failed to allocate LOS xpart memory.");
+ } else {
+ if ((LOS_parts = (struct part *)swift_malloc(
+ "los_parts_array",
+ sizeof(struct part) * LOS_list[j].particles_in_los_local)) ==
+ NULL)
+ error("Failed to allocate LOS part memory.");
+ if ((LOS_xparts = (struct xpart *)swift_malloc(
+ "los_xparts_array",
+ sizeof(struct xpart) * LOS_list[j].particles_in_los_local)) ==
+ NULL)
+ error("Failed to allocate LOS xpart memory.");
+ }
+
+ /* Loop over each part again, pulling out those in LOS. */
+ int count = 0;
+
+ for (int n = 0; n < e->s->nr_local_cells_with_particles; ++n) {
+
+ if (los_cells_top[n] == 0) continue;
+
+ const struct cell *c = &cells[local_cells_with_particles[n]];
+ const struct part *cell_parts = c->hydro.parts;
+ const struct xpart *cell_xparts = c->hydro.xparts;
+ const size_t num_parts_in_cell = c->hydro.count;
+
+ for (size_t i = 0; i < num_parts_in_cell; i++) {
+
+ /* Don't consider inhibited parts. */
+ if (cell_parts[i].time_bin == time_bin_inhibited) continue;
+
+ /* Don't consider part if outwith allowed z-range. */
+ if (cell_parts[i].x[LOS_list[j].zaxis] <
+ LOS_list[j].range_when_shooting_down_axis[0] ||
+ cell_parts[i].x[LOS_list[j].zaxis] >
+ LOS_list[j].range_when_shooting_down_axis[1])
+ continue;
+
+ /* Distance from this part to LOS along x dim. */
+ double dx = cell_parts[i].x[LOS_list[j].xaxis] - LOS_list[j].Xpos;
+
+ /* Periodic wrap. */
+ if (LOS_list[j].periodic)
+ dx = nearest(dx, LOS_list[j].dim[LOS_list[j].xaxis]);
+
+ /* Square */
+ const double dx2 = dx * dx;
+
+ /* Smoothing length of this part. */
+ const double hsml = cell_parts[i].h * kernel_gamma;
+ const double hsml2 = hsml * hsml;
+
+ /* Does this part fall into our LOS? */
+ if (dx2 < hsml2) {
+
+ /* Distance from this part to LOS along y dim. */
+ double dy = cell_parts[i].x[LOS_list[j].yaxis] - LOS_list[j].Ypos;
+
+ /* Periodic wrap. */
+ if (LOS_list[j].periodic)
+ dy = nearest(dy, LOS_list[j].dim[LOS_list[j].yaxis]);
+
+ /* Square */
+ const double dy2 = dy * dy;
+
+ /* Does this part still fall into our LOS? */
+ if (dy2 < hsml2) {
+ /* 2D distance to LOS. */
+
+ if (dx2 + dy2 <= hsml2) {
+
+ /* Store part and xpart properties. */
+ memcpy(&LOS_parts[count], &cell_parts[i], sizeof(struct part));
+ memcpy(&LOS_xparts[count], &cell_xparts[i], sizeof(struct xpart));
+
+ count++;
+ }
+ }
+ }
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (count != LOS_list[j].particles_in_los_local)
+ error("LOS counts don't add up");
+#endif
+
+#ifdef WITH_MPI
+ /* Collect all parts in this LOS to rank 0. */
+ if (e->nodeID == 0) {
+ MPI_Gatherv(MPI_IN_PLACE, 0, part_mpi_type, LOS_parts, counts, offsets,
+ part_mpi_type, 0, MPI_COMM_WORLD);
+ MPI_Gatherv(MPI_IN_PLACE, 0, xpart_mpi_type, LOS_xparts, counts, offsets,
+ xpart_mpi_type, 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Gatherv(LOS_parts, LOS_list[j].particles_in_los_local, part_mpi_type,
+ LOS_parts, counts, offsets, part_mpi_type, 0, MPI_COMM_WORLD);
+ MPI_Gatherv(LOS_xparts, LOS_list[j].particles_in_los_local,
+ xpart_mpi_type, LOS_xparts, counts, offsets, xpart_mpi_type,
+ 0, MPI_COMM_WORLD);
+ }
+#endif
+
+ /* Rank 0 writes particles to file. */
+ if (e->nodeID == 0) {
+
+ /* Create HDF5 group for this LOS */
+ sprintf(groupName, "/LOS_%04i", j);
+ h_grp =
+ H5Gcreate(h_file, groupName, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating LOS HDF5 group\n");
+
+ /* Record this LOS attributes */
+ io_write_attribute(h_grp, "NumParts", INT,
+ &LOS_list[j].particles_in_los_total, 1);
+ io_write_attribute(h_grp, "Xaxis", INT, &LOS_list[j].xaxis, 1);
+ io_write_attribute(h_grp, "Yaxis", INT, &LOS_list[j].yaxis, 1);
+ io_write_attribute(h_grp, "Zaxis", INT, &LOS_list[j].zaxis, 1);
+ io_write_attribute(h_grp, "Xpos", DOUBLE, &LOS_list[j].Xpos, 1);
+ io_write_attribute(h_grp, "Ypos", DOUBLE, &LOS_list[j].Ypos, 1);
+
+ /* Write the data for this LOS */
+ write_los_hdf5_datasets(h_grp, j, LOS_list[j].particles_in_los_total,
+ LOS_parts, e, LOS_xparts);
+
+ /* Close HDF5 group */
+ H5Gclose(h_grp);
+ }
+
+ /* Free up some memory */
+#ifdef WITH_MPI
+ free(counts);
+ free(offsets);
+#endif
+ swift_free("tl_cells_los", los_cells_top);
+ swift_free("los_parts_array", LOS_parts);
+ swift_free("los_xparts_array", LOS_xparts);
+
+ } /* End of loop over each LOS */
+
+ if (e->nodeID == 0) {
+ /* Write header */
+ write_hdf5_header(h_file, e, LOS_params, total_num_parts_in_los);
+
+ /* Close HDF5 file */
+ H5Fclose(h_file);
+ }
+
+ /* Up the LOS counter. */
+ e->los_output_count++;
+
+ /* How long did we take? */
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
+/**
+ * @brief Write a los_props struct to the given FILE as a stream of bytes.
+ *
+ * @param internal_los the struct
+ * @param stream the file stream
+ */
+void los_struct_dump(const struct los_props *internal_los, FILE *stream) {
+ restart_write_blocks((void *)internal_los, sizeof(struct los_props), 1,
+ stream, "losparams", "los params");
+}
+
+/**
+ * @brief Restore a los_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param internal_los the struct
+ * @param stream the file stream
+ */
+void los_struct_restore(const struct los_props *internal_los, FILE *stream) {
+ restart_read_blocks((void *)internal_los, sizeof(struct los_props), 1, stream,
+ NULL, "los params");
+}
diff --git a/src/line_of_sight.h b/src/line_of_sight.h
new file mode 100644
index 0000000000..0e79a9a501
--- /dev/null
+++ b/src/line_of_sight.h
@@ -0,0 +1,122 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Stuart McAlpine (stuart.mcalpine@helsinki.fi)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_LOS_H
+#define SWIFT_LOS_H
+
+/* Config parameters. */
+#include "../config.h"
+#include "engine.h"
+#include "io_properties.h"
+
+/**
+ * @brief Maps the LOS axis geometry to the simulation axis geometry.
+ *
+ * Sigtlines will always shoot down the los_direction_z,
+ * which can map to x,y or z of the simulation geometry.
+ *
+ * The remainng two axes, los_direction_x/y, then create
+ * the plane orthogonal to the LOS direction. The random
+ * sightline positions are created on this plane.
+ */
+enum los_direction { simulation_x_axis, simulation_y_axis, simulation_z_axis };
+
+/**
+ * @brief Properties of a single line-of-sight
+ */
+struct line_of_sight {
+
+ /*! Simulation axis the LOS shoots down. */
+ enum los_direction zaxis;
+
+ /*! The two remaining axes defining the plane orthogonal to the sightline. */
+ enum los_direction xaxis, yaxis;
+
+ /*! Sightline position along los_direction_x. */
+ double Xpos;
+
+ /*! Sightline position along los_direction_y. */
+ double Ypos;
+
+ /*! Number of parts in LOS. */
+ int particles_in_los_total;
+
+ /*! Number of parts in LOS on this node. */
+ int particles_in_los_local;
+
+ /*! Is the simulation periodic? */
+ int periodic;
+
+ /*! Dimensions of the space. */
+ double dim[3];
+
+ /*! How many top level cells does ths LOS intersect? */
+ int num_intersecting_top_level_cells;
+
+ /*! The min--max range to consider for parts in LOS. */
+ double range_when_shooting_down_axis[2];
+};
+
+/**
+ * @brief Properties of the line-of-sight computation
+ */
+struct los_props {
+
+ /*! Number of sightlines shooting down simulation z axis. */
+ int num_along_z;
+
+ /*! Number of sightlines shooting down simulation x axis. */
+ int num_along_x;
+
+ /*! Number of sightlines shooting down simulation y axis. */
+ int num_along_y;
+
+ /*! Total number of sightlines. */
+ int num_tot;
+
+ /*! The min--max range along the simulation x axis random sightlines are
+ * allowed. */
+ double allowed_losrange_x[2];
+
+ /*! The min--max range along the simulation y axis random sightlines are
+ * allowed. */
+ double allowed_losrange_y[2];
+
+ /*! The min--max range along the simulation z axis random sightlines are
+ * allowed. */
+ double allowed_losrange_z[2];
+
+ /*! The min--max range to consider when LOS is shooting down each simulation
+ * axis. */
+ double range_when_shooting_down_axis[3][2];
+
+ /*! Base name for line of sight HDF5 files. */
+ char basename[200];
+};
+
+void print_los_info(const struct line_of_sight *Los, const int i);
+void do_line_of_sight(struct engine *e);
+void los_init(const double dim[3], struct los_props *los_params,
+ struct swift_params *params);
+
+void los_struct_dump(const struct los_props *internal_los, FILE *stream);
+void los_struct_restore(const struct los_props *internal_los, FILE *stream);
+
+#endif /* SWIFT_LOS_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 5f26423800..be820e2ecd 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1139,6 +1139,7 @@ void prepare_file(struct engine* e, const char* fileName,
swift_type_count);
io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
io_write_attribute_i(h_grp, "ThisFile", 0);
+ io_write_attribute_s(h_grp, "OutputType", "Snapshot");
/* Close header */
H5Gclose(h_grp);
diff --git a/src/parser.h b/src/parser.h
index 87e618503d..d5445f374c 100644
--- a/src/parser.h
+++ b/src/parser.h
@@ -32,7 +32,7 @@
/* Some constants. */
#define PARSER_MAX_LINE_SIZE 256
-#define PARSER_MAX_NO_OF_PARAMS 512
+#define PARSER_MAX_NO_OF_PARAMS 600
#define PARSER_MAX_NO_OF_SECTIONS 64
/* A parameter in the input file */
diff --git a/src/part.c b/src/part.c
index 4e811de4ca..be685f1d86 100644
--- a/src/part.c
+++ b/src/part.c
@@ -435,6 +435,7 @@ MPI_Datatype xpart_mpi_type;
MPI_Datatype gpart_mpi_type;
MPI_Datatype spart_mpi_type;
MPI_Datatype bpart_mpi_type;
+MPI_Datatype lospart_mpi_type;
/**
* @brief Registers MPI particle types.
@@ -481,5 +482,6 @@ void part_free_mpi_types(void) {
MPI_Type_free(&gpart_mpi_type);
MPI_Type_free(&spart_mpi_type);
MPI_Type_free(&bpart_mpi_type);
+ MPI_Type_free(&lospart_mpi_type);
}
#endif
diff --git a/src/part.h b/src/part.h
index 3f2d637115..b01618dae8 100644
--- a/src/part.h
+++ b/src/part.h
@@ -148,6 +148,7 @@ extern MPI_Datatype xpart_mpi_type;
extern MPI_Datatype gpart_mpi_type;
extern MPI_Datatype spart_mpi_type;
extern MPI_Datatype bpart_mpi_type;
+extern MPI_Datatype lospart_mpi_type;
void part_create_mpi_types(void);
void part_free_mpi_types(void);
diff --git a/src/runner_black_holes.c b/src/runner_black_holes.c
index 2414c468f0..422b9f1e6c 100644
--- a/src/runner_black_holes.c
+++ b/src/runner_black_holes.c
@@ -320,9 +320,6 @@ void runner_do_bh_swallow(struct runner *r, struct cell *c, int timer) {
const long long swallow_id =
black_holes_get_bpart_swallow_id(&cell_bp->merger_data);
- /* message("OO id=%lld swallow_id = %lld", cell_bp->id, */
- /* swallow_id); */
-
/* Has this particle been flagged for swallowing? */
if (swallow_id >= 0) {
diff --git a/src/runner_doiact_functions_black_holes.h b/src/runner_doiact_functions_black_holes.h
index c190603c8a..41b573eeec 100644
--- a/src/runner_doiact_functions_black_holes.h
+++ b/src/runner_doiact_functions_black_holes.h
@@ -102,7 +102,8 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, with_cosmology, cosmo,
- e->gravity_properties, ti_current, e->time);
+ e->gravity_properties, e->black_holes_properties,
+ ti_current, e->time);
}
} /* loop over the parts in ci. */
} /* loop over the bparts in ci. */
@@ -157,7 +158,7 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
if (r2 < hig2) {
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, e->gravity_properties,
- ti_current);
+ e->black_holes_properties, ti_current);
}
} /* loop over the bparts in ci. */
} /* loop over the bparts in ci. */
@@ -255,7 +256,8 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, with_cosmology, cosmo,
- e->gravity_properties, ti_current, e->time);
+ e->gravity_properties, e->black_holes_properties,
+ ti_current, e->time);
}
} /* loop over the parts in cj. */
} /* loop over the bparts in ci. */
@@ -310,7 +312,7 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
if (r2 < hig2) {
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, e->gravity_properties,
- ti_current);
+ e->black_holes_properties, ti_current);
}
} /* loop over the bparts in cj. */
} /* loop over the bparts in ci. */
@@ -422,7 +424,8 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, with_cosmology, cosmo,
- e->gravity_properties, ti_current, e->time);
+ e->gravity_properties, e->black_holes_properties,
+ ti_current, e->time);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -499,7 +502,8 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, pj->h, bi, pj, xpj, with_cosmology, cosmo,
- e->gravity_properties, ti_current, e->time);
+ e->gravity_properties, e->black_holes_properties,
+ ti_current, e->time);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index a1f38c9936..db655da20a 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -827,10 +827,6 @@ void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
if (bpart_is_active(bp, e)) {
- /* Compute the final operations for repositioning of this BH */
- black_holes_end_reposition(bp, e->black_holes_properties,
- e->physical_constants, e->cosmology);
-
/* Get particle time-step */
double dt;
if (with_cosmology) {
@@ -844,6 +840,10 @@ void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
dt = get_timestep(bp->time_bin, e->time_base);
}
+ /* Compute the final operations for repositioning of this BH */
+ black_holes_end_reposition(bp, e->black_holes_properties,
+ e->physical_constants, e->cosmology, dt);
+
/* Compute variables required for the feedback loop */
black_holes_prepare_feedback(bp, e->black_holes_properties,
e->physical_constants, e->cosmology,
diff --git a/src/serial_io.c b/src/serial_io.c
index 58b8f28451..478058299a 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -1006,6 +1006,7 @@ void write_output_serial(struct engine* e,
swift_type_count);
io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
io_write_attribute_i(h_grp, "ThisFile", 0);
+ io_write_attribute_s(h_grp, "OutputType", "Snapshot");
/* Close header */
H5Gclose(h_grp);
diff --git a/src/sincos.h b/src/sincos.h
new file mode 100644
index 0000000000..fbbee10fe5
--- /dev/null
+++ b/src/sincos.h
@@ -0,0 +1,84 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_SINCOS_H
+#define SWIFT_SINCOS_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "inline.h"
+
+#if !defined(HAVE_SINCOS) && !defined(HAVE___SINCOS)
+
+/**
+ * @brief Compute both the sin() and cos() of a number.
+ *
+ * This function is only used as a replacement for compilers that do
+ * not implement GNU extensions to the C language.
+ *
+ * @param x The input value.
+ * @param sin (return) The sine of x.
+ * @param cos (return) The cosine of x.
+ */
+__attribute__((always_inline)) INLINE static void sincos(const double x,
+ double *sin,
+ double *cos) {
+
+ *sin = sin(x);
+ *cos = cos(x);
+}
+
+#endif
+
+#if !defined(HAVE_SINCOSF) && !defined(HAVE___SINCOSF)
+
+/**
+ * @brief Compute both the sin() and cos() of a number.
+ *
+ * This function is only used as a replacement for compilers that do
+ * not implement GNU extensions to the C language.
+ *
+ * @param x The input value.
+ * @param sin (return) The sine of x.
+ * @param cos (return) The cosine of x.
+ */
+__attribute__((always_inline)) INLINE static void sincosf(const float x,
+ float *sin,
+ float *cos) {
+
+ *sin = sinf(x);
+ *cos = cosf(x);
+}
+
+#endif
+
+/* Use the __sincos and __sincosf versions if needed. */
+#if !defined(HAVE_SINCOS) && defined(HAVE___SINCOS)
+#define sincos(x, s, c) __sincos(x, s, c)
+#endif
+
+#if !defined(HAVE_SINCOSF) && defined(HAVE___SINCOSF)
+#define sincosf(x, s, c) __sincosf(x, s, c)
+#endif
+
+#endif /* SWIFT_SINCOS_H */
diff --git a/src/single_io.c b/src/single_io.c
index 2cc44cc59a..d25da91279 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -851,6 +851,7 @@ void write_output_single(struct engine* e,
swift_type_count);
io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
io_write_attribute_i(h_grp, "ThisFile", 0);
+ io_write_attribute_s(h_grp, "OutputType", "Snapshot");
/* Close header */
H5Gclose(h_grp);
diff --git a/src/swift.h b/src/swift.h
index c7a5f8d5ac..3c659aa7dc 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -50,6 +50,7 @@
#include "hashmap.h"
#include "hydro.h"
#include "hydro_properties.h"
+#include "line_of_sight.h"
#include "lock.h"
#include "logger.h"
#include "logger_io.h"
diff --git a/src/task.c b/src/task.c
index 9440bf8887..70741af01b 100644
--- a/src/task.c
+++ b/src/task.c
@@ -1373,6 +1373,10 @@ enum task_categories task_get_category(const struct task *t) {
case task_type_end_grav_force:
return task_category_gravity;
+ case task_type_fof_self:
+ case task_type_fof_pair:
+ return task_category_fof;
+
case task_type_self:
case task_type_pair:
case task_type_sub_self:
diff --git a/tools/analyse_runtime.py b/tools/analyse_runtime.py
index 363b142a12..bc2addeece 100755
--- a/tools/analyse_runtime.py
+++ b/tools/analyse_runtime.py
@@ -118,6 +118,7 @@ labels = [
["fof_search_tree:", 2],
["engine_launch: \(fof\)", 2],
["engine_launch: \(fof comms\)", 2],
+ ["do_line_of_sight:", 0],
]
times = np.zeros(len(labels))
counts = np.zeros(len(labels))
--
GitLab