diff --git a/README b/README
index 97e3d4827671c7f1aa1ea8685bed7979c7cd5ce5..8587a13b3583b9e1ccb1612b0a1abf905b7a2218 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 99af7ba4f745197dcf440c39f835ff8688d62d26..7223f4719ab99fc01c103c6813ca319f89424a49 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 2181f2167b9372b6471981c490f71b9a77131880..46d9fb2a8a09fa3c010899b3e9496095c00a7a9a 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 a5f9e7eff5c72a0c968691d14244fec018217f97..2de6c6bb254aadc940db35de6806fc7785152f3e 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 dda984293e7cce7aeffd9886be8c4a633ee03c14..f56126f1ae9142a9b67c00885a76fa9e873ab2fc 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 0000000000000000000000000000000000000000..42558cee4e2daf5bfd3155aa936374a4a9606ed7
--- /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 ee854dcc89f42aa441dd235a68668d60d8df468a..61cd9c68bca0026292e398bf31f25fbe067eb085 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 ae75332e7a6caa5d1116a53ce0af642d3a22f01b..381f857728399024afe6a9b97c6e50935ea6d80d 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 4ca7e0b5bf6b50fdd6eabee36622a874efb01dd2..0000000000000000000000000000000000000000
--- 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 7945009d7dfac7dca8025c3e374926455cc9986b..ca539e3dcf262fc025d8e182e45c541c22a8c7d7 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 e908b552de266a29b34bee18e462987c2f2bfffe..78dc4692fa887ede0a0dd40d7dbf0156699bc84f 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 0ca0f37bcfb8a5290790bc604ea9eb1bf11dbad2..e598e96436f0fc3391e46e6bd0fcd0226fb2838e 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 17dcb98069ced54348b9c4050fa68462c1c41c18..5b5ca5a4e917b6d164fb84b6e00351f9ffeb7cf4 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 b4aa82a79cf77c5d19c5f2f475ba5c968216758b..825669a8e2d6b065dc8ae3da5869be577f62bf2c 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 18e8f5cfce0600fad349a66f650d9b1557e9d5ff..2afd0de07b42c2072932ff091ac0ea58a820eccf 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 0000000000000000000000000000000000000000..3b68fbf9cafed56e52f46a0194900ebdbcbcfe8f
--- /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 0000000000000000000000000000000000000000..8e0f6b6114de1d2526ca0cc6f4e7d4dd3c8923fa
--- /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 0000000000000000000000000000000000000000..fd1f7120ba2739a6d2955f31839ce7ad42979e19
--- /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 0000000000000000000000000000000000000000..592ab8483d015fe1bfafe5cc603fabc230b25589
--- /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 0000000000000000000000000000000000000000..3022d6f44f3964c02fa56d24550e7b0baa80c8ea
--- /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 0000000000000000000000000000000000000000..8590cbf5bc77e8d7a956d210339cced4bbdc692c
--- /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 c6ed9b35f687dcc116f8eba73a3cc5090ae45d6e..fa0deffdf894bce6534845df872776d90cb20dd8 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 f77d2d177df2c4a4dec77d8a37b8dc9657008bfe..8c2b1b36f0843a9bdadd5a01186d0cbcd7c46621 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 32a28d2697b0e779ceb43173f2c49b198e598c31..622a60dd1cf381360140059b19e359fec651638d 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 b5649c9dfe17a6bb44d86329d9beecfaba28133d..005fcf566036edc6d1f0b53e9130e5ba95a31355 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 d4c9890c4e9fefd86c3363745225a6d67c5c3a65..dba44b7e1422c5a79b0ebfff2e1d0cbbabd15427 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 32a5fbf66d96e68203b99d6cd920f98d6f75e090..9eeac77e732eeeb32c5440dd5586199c401fedf1 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 5f7ac3ff9e7442348185f6ba81036613eae1c81b..bf23cff630db22814b1bb4deeb54a7d17d46d6fd 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 6204d2d9caa9284b7c5e762f305deb7d989bcb7e..6b5841eed72b2be8b842bba5487455aaaa5f4ed7 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 4f10f6bc0c7dcd561bf1b3151333dbc2799794d7..759349e1a88a42ea78d3c6482e094b73cdf7669a 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 facf21d15cd1768f6af54789c58a4b718fc41d45..457c903b4d00ee6341f49dc2e2eb2226d7da70cc 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 929524ebffb418136c7a753b831fd1c1d6e12f63..e7942286214da45e821863db9874390a6d23a08b 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 bd9542f00c94ddd413d0d3e08703f1d00e45f280..67f36ede20a0dbbe48c92a67e02a2b97fc80c86a 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 a731361211f402c7d9f733f053c481743ee50620..e84bbd2991e438637a0d3f597562f5c2e959e75d 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 b0626bb6469bd626f9c4570bc113bc9af94a92eb..30186f892e1eaf9b34e07ddcc5a419312696a78b 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 9bd5ac199e45d0267e89a54009f0eae71c5d26c7..31fed833aa94298b6c25943a4212f03042fdefbb 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 cb6438bc85b4f113e5f1b1f6122ebf9fcc7088c6..c34ce6eb9e874aa8b330da050399ff9f6d0df9ea 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 0d19657ebc52b6498c886c0f55d21ed53db6f2bb..ed2f3bca72d2c548aed20915fdb7295e2b9e5a4d 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 54068b34840da1e4c829474ab1ca2614017c9f8f..fba2df4ba6576dd5fc736530c1ddf0485044702c 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 e64218588293720a5caa164a2f9b456b8ab02a8b..d6ab4f225423284bfedcfe2fde90c285f764a5ab 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 69424fa331b60964fd7448b08ecc59d7751a492d..92efcafa4ee25006d0651d9f97db82e85c88b3a1 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 0ecb2ebe517f322f32ee4701d63d2153cfd2900d..4fd7bc06a5cbd79f84a9fb04ac894b8a59ff806a 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 1597687bf6feec49f88f8a2f2b4bba05ecbf5939..b8f28f5999d303d8970293026ade508ed9d0d597 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 091d567066a9f490050d77ce27f96d4140ecb39d..d68543ae6caeb59b348a16d3b5d8e206849cd1b7 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 848a90dd5edf36c845cde02b14b1cc6996c40405..d24db8f99df4324be88bcab0d98e2ed42f75f74b 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 8a5a21f4fc555d4d27b33b34853669f2010b5e22..355fdb05ac33b40ffdb0c2b9630b37a2e4b136c5 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 3ac99c047ce3b99bb8e1c9bfc1389ac0af93aef7..4fc82e59c954b91467f0b2a55aa0eacb9f6b97e2 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 3543667a563587e126cf191968d93c7c5e7b4e0a..e1ace452ecec54ff28d1fa6b474ff064183c2897 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 329563d4e6c4b71e01e5271137da37a5a8b99db5..281ad02715b06dead893f7cbc679782803b383ca 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 c64ffada10891f98a0c6bcb9ab3b412e276c3b3b..065a38e1969c48682b7c4619cc317e2dd576e7e2 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 2172dae966036c3bd62e9481068011caa5e57689..d195c68b15278cf2758ec307541d0422606a5bd9 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 e49575e769461de259285d92433f6dc1f6ac266c..36eb28ccf4eb804130e7930bb5a99e8826424c02 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 fddafa52a93fe5823c7808660858b83768efec85..6f68f92d1a11403e3398ba9b3ab19e87601024c4 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 27508dd7ea3db86cf4c65e3762694488a53616bf..11b4659086d3c2b0791789a411934e1ff8b22c32 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 fd19256ebf2c6fc746dfeb367ae242794bcb69cb..f0504b382e911ca85c2c6e68d54f6dfc50841a94 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 b0da3ac888ec1a7421ee19471bb5736a01130682..ed158d7d10fb720d8c638b1bd76c759741ce11e4 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 d77693a50e2b33c4e5917badbd937af50da4892e..06b77e9eecaebccea5fc444e21772025d961af87 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 0000000000000000000000000000000000000000..cbd5426fb0facc37e528a8b47c1a8dc0ace1a785
--- /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 0000000000000000000000000000000000000000..0e79a9a50126a76cee152d9fc9612b7b33961407
--- /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 5f26423800357ecd8d5a8c81a94dfb047ff43e13..be820e2ecd3434f7d41635f8b0397ff5a1ef368c 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 87e618503dca567f7cf381475fcc91f1575b3913..d5445f374cd4b1e0808c188e6bcae6a159728695 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 4e811de4ca9812358fa2182ec8178be04ac810dc..be685f1d862e5fe69fd45c73998fbfc0a0b3b2ff 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 3f2d6371156e5e259f095f88ac744be50be845eb..b01618dae8f659af93eec6c71bf2a7305e945f91 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 2414c468f0dd1f031cbaedb5612ab90684a27e59..422b9f1e6cbabbc263c5bb927097dc29ed0ea9e2 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 c190603c8aac663314b505af0d8c69ba1896541f..41b573eeecdfa4332d266161c5d6cecfebbdae83 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 a1f38c9936578e63a31097cdd1b577c8fc6f762b..db655da20ad5475da888334fdeed2f0665acac30 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 58b8f284515086866c9dd4da4f0879800d5ee254..478058299a7b777b1351adcf64ddbf96b7bfb2bc 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 0000000000000000000000000000000000000000..fbbee10fe54ce5ac48d1eaaca083cdc0eb51f52c
--- /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 2cc44cc59ad376b483bae87ebe922989ac1b9c8d..d25da91279e5f89ed456c9519d913a39ff05596a 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 c7a5f8d5acbcc17dd467cbcfb914eac025894dae..3c659aa7dce1d78d9887187da31e4ee500ad6848 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 9440bf88876198487be31ef44cd8f7308c6af9e0..70741af01bbf869ec22b1f5fd8ac5c076784db2a 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 363b142a12caac910519129dab6258bac0022add..bc2addeece6947a9da4ef318db1964c32d285441 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))