diff --git a/configure.ac b/configure.ac
index 1130f5d875ccbaa6450371640f0156cefd9a56ca..71f3e0d5ec7c07fd18219d84e61114359939fdc5 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1340,7 +1340,7 @@ case "$with_subgrid" in
with_subgrid_stars=EAGLE
with_subgrid_star_formation=EAGLE
with_subgrid_feedback=EAGLE
- with_subgrid_black_holes=none
+ with_subgrid_black_holes=EAGLE
;;
*)
AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
@@ -1789,8 +1789,11 @@ case "$with_black_holes" in
none)
AC_DEFINE([BLACK_HOLES_NONE], [1], [No black hole model])
;;
+ EAGLE)
+ AC_DEFINE([BLACK_HOLES_EAGLE], [1], [EAGLE black hole model])
+ ;;
*)
- AC_MSG_ERROR([Unknown stellar model: $with_black_holes])
+ AC_MSG_ERROR([Unknown black-hole model: $with_black_holes])
;;
esac
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index a3da1b4f8edc8c1af225ac941441a92b14a13f99..c6b3046d2d3591c937dfd98cf75fb7697b90110f 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -776,7 +776,7 @@ INPUT += @top_srcdir@/src/star_formation/EAGLE
INPUT += @top_srcdir@/src/tracers/EAGLE
INPUT += @top_srcdir@/src/stars/EAGLE
INPUT += @top_srcdir@/src/feedback/EAGLE
-INPUT += @top_srcdir@/src/black_holes/Default
+INPUT += @top_srcdir@/src/black_holes/EAGLE
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 246d0d4834ff8ce4d08421d6fa754dd81f0b1383..4dcc1e0381adcb0405dd5b502c3d64896269a10f 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -147,3 +147,11 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index f6693c89b0d281e99adb36d905e2897e49495a1e..e1d9100a8f87bb8670fcb82ce4451084ecb1375a 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -148,3 +148,11 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index 5490b3a263badb6fca4b3bf134375a72c773e367..3efe5ba972903dbc3644581d214285603e75ca78 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -148,3 +148,13 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
+
+
\ No newline at end of file
diff --git a/examples/EAGLE_low_z/EAGLE_12/README b/examples/EAGLE_low_z/EAGLE_12/README
index 9fb2930f625fd2e85293b9ebcc91f46e59f02745..0530021202191bac74690e47106ca3488011d1b5 100644
--- a/examples/EAGLE_low_z/EAGLE_12/README
+++ b/examples/EAGLE_low_z/EAGLE_12/README
@@ -13,4 +13,4 @@ The particle load of the main EAGLE simulation can be reproduced by
running these ICs on 8 cores.
MD5 checksum of the ICs:
-d1ebf1c7ffb0488c3628e08ed6d2dbf4 EAGLE_ICs_12.hdf5
+6f96624175df28f2c4a02265b871b3d5 EAGLE_ICs_12.hdf5
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 873b8ec1084a78f830f6a927f173d42be9c5a210..1b81cbc213a0e12111d35e5efde31606d60fcd59 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -140,3 +140,11 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
diff --git a/examples/EAGLE_low_z/EAGLE_25/README b/examples/EAGLE_low_z/EAGLE_25/README
index 0fe37964d3734b63acd10195822e11094219586a..80907ce9ffff04382afceb7c906dca089a05b8bb 100644
--- a/examples/EAGLE_low_z/EAGLE_25/README
+++ b/examples/EAGLE_low_z/EAGLE_25/README
@@ -13,4 +13,4 @@ The particle load of the main EAGLE simulation can be reproduced by
running these ICs on 64 cores.
MD5 checksum of the ICs:
-592faecfc51239a00396f5a4acc2fa4d EAGLE_ICs_25.hdf5
+e6a5de0e962d8ffb7589671b9613daa0 EAGLE_ICs_25.hdf5
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 51bfe6af3d904c822bbc31720ad64fb95652eede..569dbdd1487faea311318fa508d5a14832f074d4 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -147,3 +147,11 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
diff --git a/examples/EAGLE_low_z/EAGLE_50/README b/examples/EAGLE_low_z/EAGLE_50/README
index 211360be35f83423160ceb9912f7f0ce9ee444c8..f63cd92438fb94dc9d34b1c048b0203fb903e020 100644
--- a/examples/EAGLE_low_z/EAGLE_50/README
+++ b/examples/EAGLE_low_z/EAGLE_50/README
@@ -13,4 +13,4 @@ The particle load of the main EAGLE simulation can be reproduced by
running these ICs on 512 cores.
MD5 checksum of the ICs:
-203985136c665e8677a4faa77f94110a EAGLE_ICs_50.hdf5
+fca1c2efea5c81f74a0d260fefb19595 EAGLE_ICs_50.hdf5
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index 07c527158d7b2dd38205a202d6a8ad5fae6a2b8f..4e1fbb8cdaf025d10ef07af09085290c91223dca 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -142,3 +142,11 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
diff --git a/examples/EAGLE_low_z/EAGLE_6/README b/examples/EAGLE_low_z/EAGLE_6/README
index 844c779e222b1cd7ebefaa61d637092a5a8c1e79..a236f338649e2199220c7f1c28e5531b1fd980c7 100644
--- a/examples/EAGLE_low_z/EAGLE_6/README
+++ b/examples/EAGLE_low_z/EAGLE_6/README
@@ -10,4 +10,4 @@ been removed.
Everything is ready to be run without cosmological integration.
MD5 checksum of the ICs:
-1923db3cf59b2e80b23e4e1aee13e012 EAGLE_ICs_6.hdf5
+ddd26b959815b131cd1af1479e752da7 EAGLE_ICs_6.hdf5
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 3632fcd82d50cf3f1d53b40bf329a1dbb35124fe..978f21ea5f7df1c1bb13329cbea796254d3d3b2a 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -68,7 +68,7 @@ InitialConditions:
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
-
+
EAGLEChemistry: # Solar abundances
init_abundance_metal: 0.014
init_abundance_Hydrogen: 0.70649785
@@ -151,3 +151,12 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
+
diff --git a/examples/main.c b/examples/main.c
index 08e8eacf5510bf7af77175d9a74460651f656800..4efd82d799f339a337df019c4f9ef12cbff07075 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -98,6 +98,7 @@ int main(int argc, char *argv[]) {
struct stars_props stars_properties;
struct feedback_props feedback_properties;
struct entropy_floor_properties entropy_floor;
+ struct black_holes_props black_holes_properties;
struct part *parts = NULL;
struct phys_const prog_const;
struct space s;
@@ -389,6 +390,16 @@ int main(int argc, char *argv[]) {
return 1;
}
+ if (!with_hydro && with_black_holes) {
+ if (myrank == 0) {
+ argparse_usage(&argparse);
+ printf(
+ "\nError: Cannot process black holes 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 &&
@@ -754,6 +765,16 @@ int main(int argc, char *argv[]) {
} else
bzero(&feedback_properties, sizeof(struct feedback_props));
+ /* Initialise the black holes properties */
+ if (with_black_holes) {
+#ifdef BLACK_HOLES_NONE
+ error("ERROR: Running with black_holes but compiled without it.");
+#endif
+ black_holes_props_init(&black_holes_properties, &prog_const, &us, params,
+ &hydro_properties, &cosmo);
+ } else
+ bzero(&black_holes_properties, sizeof(struct black_holes_props));
+
/* Initialise the gravity properties */
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology,
@@ -866,7 +887,7 @@ int main(int argc, char *argv[]) {
long long N_total[4] = {0, 0, 0, 0};
#if defined(WITH_MPI)
long long N_long[4] = {Ngas, Ngpart, Nspart, Nbpart};
- MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
+ MPI_Allreduce(&N_long, &N_total, 4, MPI_LONG_LONG_INT, MPI_SUM,
MPI_COMM_WORLD);
#else
N_total[0] = Ngas;
@@ -1004,8 +1025,9 @@ int main(int argc, char *argv[]) {
engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
&hydro_properties, &entropy_floor, &gravity_properties,
- &stars_properties, &feedback_properties, &mesh, &potential,
- &cooling_func, &starform, &chemistry);
+ &stars_properties, &black_holes_properties,
+ &feedback_properties, &mesh, &potential, &cooling_func,
+ &starform, &chemistry);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
talking, restart_file);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index f0a3c398e3afb895fb1932b654816437f2b0f429..375c5abc22d117b7322dd49e37109430b3ef0e32 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -374,3 +374,12 @@ EAGLEFeedback:
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:
+ max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ 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.
+
diff --git a/src/Makefile.am b/src/Makefile.am
index 3f7a61bc097802e3e28744234e0c96116918187c..adf9f30ae367854a6742adb333005b0d11b9d8b4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -53,7 +53,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
star_formation_struct.h star_formation.h star_formation_iact.h \
star_formation_logger.h star_formation_logger_struct.h \
velociraptor_struct.h velociraptor_io.h random.h memuse.h black_holes.h black_holes_io.h \
- feedback.h feedback_struct.h feedback_properties.h
+ black_holes_properties.h feedback.h feedback_struct.h feedback_properties.h
# source files for EAGLE cooling
EAGLE_COOLING_SOURCES =
@@ -84,9 +84,9 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
gravity_iact.h kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h \
- runner_doiact_nosort.h runner_doiact_stars.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h \
- adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
- logger_io.h timestep_limiter.h \
+ runner_doiact_nosort.h runner_doiact_stars.h runner_doiact_black_holes.h units.h intrinsics.h minmax.h \
+ kick.h timestep.h drift.h adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h \
+ dump.h logger.h sign.h logger_io.h timestep_limiter.h \
gravity.h gravity_io.h gravity_cache.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
@@ -196,7 +196,11 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
feedback/EAGLE/feedback_properties.h feedback/EAGLE/imf.h feedback/EAGLE/interpolate.h \
feedback/EAGLE/yield_tables.h \
black_holes/Default/black_holes.h black_holes/Default/black_holes_io.h \
- black_holes/Default/black_holes_part.h
+ black_holes/Default/black_holes_part.h black_holes/Default/black_holes_iact.h \
+ black_holes/Default/black_holes_properties.h \
+ black_holes/EAGLE/black_holes.h black_holes/EAGLE/black_holes_io.h \
+ black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h \
+ black_holes/EAGLE/black_holes_properties.h
# Sources and flags for regular library
diff --git a/src/active.h b/src/active.h
index a0f793b5531eaa19f6d1e92a5f243d753759d277..e7cd35fa570f15694ab8fa3011d6a02b65837886 100644
--- a/src/active.h
+++ b/src/active.h
@@ -372,7 +372,7 @@ __attribute__((always_inline)) INLINE static int bpart_is_active(
if (ti_end < ti_current)
error(
- "s-particle in an impossible time-zone! bp->ti_end=%lld "
+ "b-particle in an impossible time-zone! bp->ti_end=%lld "
"e->ti_current=%lld",
ti_end, ti_current);
#endif
diff --git a/src/black_holes.h b/src/black_holes.h
index 0decc3fbf5eddd30cc5b5a15b73d0a4c2447e5ee..8b658e45f37e4e49e1a14f1e3e853b8dbed476d8 100644
--- a/src/black_holes.h
+++ b/src/black_holes.h
@@ -25,6 +25,10 @@
/* Select the correct star model */
#if defined(BLACK_HOLES_NONE)
#include "./black_holes/Default/black_holes.h"
+#include "./black_holes/Default/black_holes_iact.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes.h"
+#include "./black_holes/EAGLE/black_holes_iact.h"
#else
#error "Invalid choice of star model"
#endif
diff --git a/src/black_holes/Default/black_holes.h b/src/black_holes/Default/black_holes.h
index f0f53bdf46cdd6afe062777c993e333ac5ef56c6..3c35d0862df99ff194fd609d6058c1be7065b0bf 100644
--- a/src/black_holes/Default/black_holes.h
+++ b/src/black_holes/Default/black_holes.h
@@ -20,6 +20,9 @@
#define SWIFT_DEFAULT_BLACK_HOLES_H
#include <float.h>
+
+/* Local includes */
+#include "black_holes_properties.h"
#include "dimension.h"
#include "kernel_hydro.h"
#include "minmax.h"
@@ -42,9 +45,10 @@ __attribute__((always_inline)) INLINE static float black_holes_compute_timestep(
* read in to do some conversions.
*
* @param bp The particle to act upon
+ * @param props The properties of the black holes model.
*/
__attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
- struct bpart* bp) {
+ struct bpart* bp, const struct black_holes_props* props) {
bp->time_bin = 0;
}
@@ -135,6 +139,23 @@ black_holes_bpart_has_no_neighbours(struct bpart* restrict bp,
bp->density.wcount_dh = 0.f;
}
+/**
+ * @brief Compute the accretion rate of the black hole and all the quantites
+ * required for the feedback loop.
+ *
+ * Nothing to do here.
+ *
+ * @param bp The black hole particle.
+ * @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 time-step size (in physical internal units).
+ */
+__attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
+ struct bpart* restrict bp, const struct black_holes_props* props,
+ 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
new file mode 100644
index 0000000000000000000000000000000000000000..4aeb27b4e13a614036962a564d99a80c6990ebec
--- /dev/null
+++ b/src/black_holes/Default/black_holes_iact.h
@@ -0,0 +1,96 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 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_DEFAULT_BH_IACT_H
+#define SWIFT_DEFAULT_BH_IACT_H
+
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @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 cosmo The cosmological model.
+ * @param ti_current Current integer time value (for random numbers).
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct bpart *restrict bi, const struct part *restrict pj,
+ const struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const integertime_t ti_current) {
+
+ float wi, wi_dx;
+
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Compute the kernel function */
+ const float hi_inv = 1.0f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute contribution to the number of neighbours */
+ bi->density.wcount += wi;
+ bi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+#ifdef DEBUG_INTERACTIONS_BH
+ /* Update ngb counters */
+ if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_BH)
+ bi->ids_ngbs_density[si->num_ngb_density] = pj->id;
+
+ /* Update ngb counters */
+ ++si->num_ngb_density;
+#endif
+}
+
+/**
+ * @brief Feedback interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param bi First particle (black hole).
+ * @param pj Second particle (gas)
+ * @param xpj The extended data of the second particle.
+ * @param cosmo The cosmological model.
+ * @param ti_current Current integer time value (for random numbers).
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_bh_feedback(const float r2, const float *dx, const float hi,
+ const float hj, struct bpart *restrict bi,
+ struct part *restrict pj,
+ const struct xpart *restrict xpj,
+ const struct cosmology *cosmo,
+ const integertime_t ti_current) {
+#ifdef DEBUG_INTERACTIONS_BH
+ /* Update ngb counters */
+ if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
+ bi->ids_ngbs_force[si->num_ngb_force] = pj->id;
+
+ /* Update ngb counters */
+ ++si->num_ngb_force;
+#endif
+}
+
+#endif
diff --git a/src/black_holes/Default/black_holes_properties.h b/src/black_holes/Default/black_holes_properties.h
new file mode 100644
index 0000000000000000000000000000000000000000..e6e66d5b1ce14e920a498a3a97013c0d33596b38
--- /dev/null
+++ b/src/black_holes/Default/black_holes_properties.h
@@ -0,0 +1,126 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 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_DEFAULT_BLACK_HOLES_PROPERTIES_H
+#define SWIFT_DEFAULT_BLACK_HOLES_PROPERTIES_H
+
+#include "chemistry.h"
+#include "hydro_properties.h"
+
+/**
+ * @brief Properties of the black hole scheme.
+ *
+ * In this default scheme, we only have the properties
+ * required by the neighbour search.
+ */
+struct black_holes_props {
+
+ /*! Resolution parameter */
+ float eta_neighbours;
+
+ /*! Target weightd number of neighbours (for info only)*/
+ float target_neighbours;
+
+ /*! Smoothing length tolerance */
+ float h_tolerance;
+
+ /*! Tolerance on neighbour number (for info only)*/
+ float delta_neighbours;
+
+ /*! Maximal number of iterations to converge h */
+ int max_smoothing_iterations;
+
+ /*! Maximal change of h over one time-step */
+ float log_max_h_change;
+};
+
+/**
+ * @brief Initialise the black hole properties from the parameter file.
+ *
+ * We read the defaults from the hydro properties.
+ *
+ * @param bp The #black_holes_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ * @param hydro_props The already read-in properties of the hydro scheme.
+ * @param cosmo The cosmological model.
+ */
+static INLINE void black_holes_props_init(struct black_holes_props *bp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *hydro_props,
+ const struct cosmology *cosmo) {
+
+ /* Kernel properties */
+ bp->eta_neighbours = parser_get_opt_param_float(
+ params, "BlackHoles:resolution_eta", hydro_props->eta_neighbours);
+
+ /* Tolerance for the smoothing length Newton-Raphson scheme */
+ bp->h_tolerance = parser_get_opt_param_float(params, "BlackHoles:h_tolerance",
+ hydro_props->h_tolerance);
+
+ /* Get derived properties */
+ bp->target_neighbours = pow_dimension(bp->eta_neighbours) * kernel_norm;
+ const float delta_eta = bp->eta_neighbours * (1.f + bp->h_tolerance);
+ bp->delta_neighbours =
+ (pow_dimension(delta_eta) - pow_dimension(bp->eta_neighbours)) *
+ kernel_norm;
+
+ /* Number of iterations to converge h */
+ bp->max_smoothing_iterations =
+ parser_get_opt_param_int(params, "BlackHoles:max_ghost_iterations",
+ hydro_props->max_smoothing_iterations);
+
+ /* Time integration properties */
+ const float max_volume_change =
+ parser_get_opt_param_float(params, "BlackHoles:max_volume_change", -1);
+ if (max_volume_change == -1)
+ bp->log_max_h_change = hydro_props->log_max_h_change;
+ else
+ bp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+}
+
+/**
+ * @brief Write a black_holes_props struct to the given FILE as a stream of
+ * bytes.
+ *
+ * @param props the black hole properties struct
+ * @param stream the file stream
+ */
+INLINE static void black_holes_struct_dump(
+ const struct black_holes_props *props, FILE *stream) {
+ restart_write_blocks((void *)props, sizeof(struct black_holes_props), 1,
+ stream, "black_holes props", "black holes props");
+}
+
+/**
+ * @brief Restore a black_holes_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param props the black hole properties struct
+ * @param stream the file stream
+ */
+INLINE static void black_holes_struct_restore(
+ const struct black_holes_props *props, FILE *stream) {
+ restart_read_blocks((void *)props, sizeof(struct black_holes_props), 1,
+ stream, NULL, "black holes props");
+}
+
+#endif /* SWIFT_DEFAULT_BLACK_HOLES_PROPERTIES_H */
diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
new file mode 100644
index 0000000000000000000000000000000000000000..909057bb962cd9a21af9a708d8796f5cd18d6a8d
--- /dev/null
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -0,0 +1,315 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 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_EAGLE_BLACK_HOLES_H
+#define SWIFT_EAGLE_BLACK_HOLES_H
+
+/* Local includes */
+#include "black_holes_properties.h"
+#include "cosmology.h"
+#include "dimension.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+#include "physical_constants.h"
+
+/* Standard includes */
+#include <float.h>
+
+/**
+ * @brief Computes the gravity time-step of a given black hole particle.
+ *
+ * @param bp Pointer to the s-particle data.
+ */
+__attribute__((always_inline)) INLINE static float black_holes_compute_timestep(
+ const struct bpart* const bp) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Initialises the b-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param bp The particle to act upon
+ * @param props The properties of the black holes model.
+ */
+__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;
+ bp->total_accreted_mass = 0.f;
+ bp->accretion_rate = 0.f;
+ bp->formation_time = -1.f;
+}
+
+/**
+ * @brief Prepares a b-particle for its interactions
+ *
+ * @param bp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_init_bpart(
+ struct bpart* bp) {
+
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ bp->ids_ngbs_density[i] = -1;
+ bp->num_ngb_density = 0;
+#endif
+
+ bp->density.wcount = 0.f;
+ bp->density.wcount_dh = 0.f;
+ bp->rho_gas = 0.f;
+ bp->sound_speed_gas = 0.f;
+ bp->velocity_gas[0] = 0.f;
+ bp->velocity_gas[1] = 0.f;
+ bp->velocity_gas[2] = 0.f;
+ bp->ngb_mass = 0.f;
+ bp->num_ngbs = 0;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param bp The particle
+ * @param dt_drift The drift time-step for positions.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_predict_extra(
+ struct bpart* restrict bp, float dt_drift) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param bp The particle.
+ */
+__attribute__((always_inline)) INLINE static void
+black_holes_reset_predicted_values(struct bpart* restrict bp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param bp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void black_holes_kick_extra(
+ struct bpart* bp, float dt) {}
+
+/**
+ * @brief Finishes the calculation of density on black holes
+ *
+ * @param bp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_end_density(
+ struct bpart* bp, const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = bp->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ 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 */
+ 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;
+
+ 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;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #spart has 0
+ * ngbs.
+ *
+ * @param bp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+black_holes_bpart_has_no_neighbours(struct bpart* restrict bp,
+ const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = bp->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ bp->density.wcount = kernel_root * h_inv_dim;
+ bp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Compute the accretion rate of the black hole and all the quantites
+ * required for the feedback loop.
+ *
+ * @param bp The black hole particle.
+ * @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 time-step size (in physical internal units).
+ */
+__attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
+ struct bpart* restrict bp, const struct black_holes_props* props,
+ const struct phys_const* constants, const struct cosmology* cosmo,
+ const double dt) {
+
+ /* Gather some physical constants (all in internal units) */
+ const double G = constants->const_newton_G;
+ const double c = constants->const_speed_light_c;
+ const double proton_mass = constants->const_proton_mass;
+ const double sigma_Thomson = constants->const_thomson_cross_section;
+
+ /* Gather the parameters of the model */
+ const double f_Edd = props->f_Edd;
+ const double epsilon_r = props->epsilon_r;
+ const double epsilon_f = props->epsilon_f;
+ const double num_ngbs_to_heat = props->num_ngbs_to_heat;
+ const double delta_T = props->AGN_delta_T_desired;
+ const double delta_u = delta_T * props->temp_to_u_factor;
+
+ /* (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)
+ */
+ 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};
+
+ /* 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];
+
+ /* 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);
+ const double Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass *
+ gas_rho_phys * denominator_inv * denominator_inv *
+ denominator_inv;
+
+ /* Compute the Eddington rate (internal units) */
+ const double Eddington_rate =
+ 4. * M_PI * G * BH_mass * proton_mass / (epsilon_r * c * sigma_Thomson);
+
+ /* Limit the accretion rate to the Eddington fraction */
+ const double accr_rate = max(Bondi_rate, f_Edd * Eddington_rate);
+ bp->accretion_rate = accr_rate;
+
+ /* Factor in the radiative efficiency */
+ const double mass_rate = (1. - epsilon_r) * accr_rate;
+ const double luminosity = epsilon_r * accr_rate * c * c;
+
+ /* Integrate forward in time */
+ bp->subgrid_mass += mass_rate * dt;
+ bp->total_accreted_mass += mass_rate * dt;
+ bp->energy_reservoir += luminosity * epsilon_f * dt;
+
+ /* Energy required to have a feedback event */
+ const double mean_ngb_mass = bp->ngb_mass / ((double)bp->num_ngbs);
+ const double E_feedback_event = num_ngbs_to_heat * delta_u * mean_ngb_mass;
+
+ /* Are we doing some feedback? */
+ if (bp->energy_reservoir > E_feedback_event) {
+
+ /* Default probability of heating */
+ double target_prob = bp->energy_reservoir / (delta_u * bp->ngb_mass);
+
+ /* Calculate the change in internal energy of the gas particles that get
+ * heated. Adjust the prbability if needed. */
+ double gas_delta_u;
+ double prob;
+ if (target_prob <= 1.) {
+
+ /* Normal case */
+ prob = target_prob;
+ gas_delta_u = delta_u;
+
+ } else {
+
+ /* Special case: we need to adjust the energy irrespective of the
+ * desired deltaT to ensure we inject all the available energy. */
+
+ prob = 1.;
+ gas_delta_u = bp->energy_reservoir / bp->ngb_mass;
+ }
+
+ /* Store all of this in the black hole for delivery onto the gas. */
+ bp->to_distribute.AGN_heating_probability = prob;
+ bp->to_distribute.AGN_delta_u = gas_delta_u;
+
+ /* Decrement the energy in the reservoir by the mean expected energy */
+ const double energy_used = bp->energy_reservoir / max(prob, 1.);
+ bp->energy_reservoir -= energy_used;
+
+ } else {
+
+ /* Flag that we don't want to heat anyone */
+ bp->to_distribute.AGN_heating_probability = 0.f;
+ bp->to_distribute.AGN_delta_u = 0.f;
+ }
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is the equivalent of hydro_reset_acceleration.
+ * We do not compute the acceleration on black hole, therefore no need to use
+ * it.
+ *
+ * @param bp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_reset_feedback(
+ struct bpart* restrict bp) {
+
+ bp->to_distribute.AGN_heating_probability = 0.f;
+ bp->to_distribute.AGN_delta_u = 0.f;
+
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ bp->ids_ngbs_force[i] = -1;
+ bp->num_ngb_force = 0;
+#endif
+}
+
+#endif /* SWIFT_EAGLE_BLACK_HOLES_H */
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..3174dea57ade8a9e3a1d9db3aa54cda428e0bc9f
--- /dev/null
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -0,0 +1,150 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 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_EAGLE_BH_IACT_H
+#define SWIFT_EAGLE_BH_IACT_H
+
+/* Local includes */
+#include "hydro.h"
+#include "random.h"
+
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @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 cosmo The cosmological model.
+ * @param ti_current Current integer time value (for random numbers).
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct bpart *restrict bi, const struct part *restrict pj,
+ const struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const integertime_t ti_current) {
+
+ float wi, wi_dx;
+
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Compute the kernel function */
+ const float hi_inv = 1.0f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute contribution to the number of neighbours */
+ bi->density.wcount += wi;
+ bi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ /* Contribution to the number of neighbours */
+ bi->num_ngbs += 1;
+
+ /* Neighbour gas mass */
+ const float mj = hydro_get_mass(pj);
+
+ /* Contribution to the BH gas density */
+ bi->rho_gas += mj * wi;
+
+ /* Contribution to the total neighbour mass */
+ bi->ngb_mass += mj;
+
+ /* Neighbour sounds speed */
+ const float cj = hydro_get_comoving_soundspeed(pj);
+
+ /* 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]};
+
+ /* 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;
+
+#ifdef DEBUG_INTERACTIONS_BH
+ /* Update ngb counters */
+ if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_BH)
+ bi->ids_ngbs_density[si->num_ngb_density] = pj->id;
+
+ /* Update ngb counters */
+ ++si->num_ngb_density;
+#endif
+}
+
+/**
+ * @brief Feedback interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param bi First particle (black hole).
+ * @param pj Second particle (gas)
+ * @param xpj The extended data of the second particle.
+ * @param cosmo The cosmological model.
+ * @param ti_current Current integer time value (for random numbers).
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_bh_feedback(const float r2, const float *dx, const float hi,
+ const float hj, struct bpart *restrict bi,
+ struct part *restrict pj,
+ struct xpart *restrict xpj,
+ const struct cosmology *cosmo,
+ const integertime_t ti_current) {
+
+ /* Get the heating probability */
+ const float prob = bi->to_distribute.AGN_heating_probability;
+
+ /* Are we doing some feedback? */
+ if (prob > 0.f) {
+
+ /* Draw a random number (Note mixing both IDs) */
+ const float rand = random_unit_interval(bi->id + pj->id, ti_current,
+ random_number_BH_feedback);
+
+ /* Are we lucky? */
+ if (rand < prob) {
+
+ /* Compute new energy of this particle */
+ const double u_init = hydro_get_physical_internal_energy(pj, xpj, cosmo);
+ const float delta_u = bi->to_distribute.AGN_delta_u;
+ const double u_new = u_init + delta_u;
+
+ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new);
+ hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new);
+ }
+ }
+
+#ifdef DEBUG_INTERACTIONS_BH
+ /* Update ngb counters */
+ if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
+ bi->ids_ngbs_force[si->num_ngb_force] = pj->id;
+
+ /* Update ngb counters */
+ ++si->num_ngb_force;
+#endif
+}
+
+#endif /* SWIFT_EAGLE_BH_IACT_H */
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc6cb962db61f9e45ded990b9dca7f00db075733
--- /dev/null
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -0,0 +1,114 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 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_EAGLE_BLACK_HOLES_IO_H
+#define SWIFT_EAGLE_BLACK_HOLES_IO_H
+
+#include "black_holes_part.h"
+#include "io_properties.h"
+
+/**
+ * @brief Specifies which b-particle fields to read from a dataset
+ *
+ * @param bparts The b-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void black_holes_read_particles(struct bpart *bparts,
+ struct io_props *list,
+ int *num_fields) {
+
+ /* Say how much we want to read */
+ *num_fields = 6;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, bparts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, bparts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ bparts, mass);
+ list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, bparts, id);
+ list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_LENGTH, bparts, h);
+ list[5] = io_make_input_field("EnergyReservoir", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_ENERGY, bparts, energy_reservoir);
+}
+
+/**
+ * @brief Specifies which b-particle fields to write to a dataset
+ *
+ * @param bparts The b-particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void black_holes_write_particles(const struct bpart *bparts,
+ struct io_props *list,
+ int *num_fields) {
+
+ /* Say how much we want to write */
+ *num_fields = 12;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
+ bparts, x);
+ list[1] =
+ io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, bparts, v);
+ list[2] =
+ io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, bparts, mass);
+ list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
+ bparts, id);
+ list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ bparts, h);
+ list[5] = io_make_output_field("SubgridMasses", FLOAT, 1, UNIT_CONV_MASS,
+ bparts, subgrid_mass);
+ list[6] = io_make_output_field("FormationTime", FLOAT, 1, UNIT_CONV_TIME,
+ bparts, formation_time);
+ list[7] = io_make_output_field("GasDensity", FLOAT, 1, UNIT_CONV_DENSITY,
+ bparts, rho_gas);
+ list[8] = io_make_output_field("GasSoundSpeed", FLOAT, 1, UNIT_CONV_SPEED,
+ bparts, sound_speed_gas);
+ list[9] = io_make_output_field("EnergyReservoir", FLOAT, 1, UNIT_CONV_ENERGY,
+ bparts, energy_reservoir);
+ list[10] = io_make_output_field("AccretionRate", FLOAT, 1,
+ UNIT_CONV_MASS_PER_UNIT_TIME, bparts,
+ accretion_rate);
+ list[11] = io_make_output_field("TotalAccretedMass", FLOAT, 1,
+ UNIT_CONV_MASS_PER_UNIT_TIME, bparts,
+ total_accreted_mass);
+
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+
+ list += *num_fields;
+ *num_fields += 4;
+
+ list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
+ bparts, num_ngb_density);
+ list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
+ bparts, num_ngb_force);
+ list[2] = io_make_output_field("Ids_ngb_density", LONGLONG,
+ MAX_NUM_OF_NEIGHBOURS_BLACK_HOLES,
+ UNIT_CONV_NO_UNITS, bparts, ids_ngbs_density);
+ list[3] = io_make_output_field("Ids_ngb_force", LONGLONG,
+ MAX_NUM_OF_NEIGHBOURS_BLACK_HOLES,
+ UNIT_CONV_NO_UNITS, bparts, ids_ngbs_force);
+#endif
+}
+
+#endif /* SWIFT_EAGLE_BLACK_HOLES_IO_H */
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..b6d15dfa06be1c4572845d1b9656b053439f70d1
--- /dev/null
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -0,0 +1,137 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 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_EAGLE_BLACK_HOLE_PART_H
+#define SWIFT_EAGLE_BLACK_HOLE_PART_H
+
+/**
+ * @brief Particle fields for the black hole particles.
+ *
+ * All quantities related to gravity are stored in the associate #gpart.
+ */
+struct bpart {
+
+ /*! Particle ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+ /*! Particle position. */
+ double x[3];
+
+ /* Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Particle velocity. */
+ float v[3];
+
+ /*! Black hole mass */
+ float mass;
+
+ /* Particle cutoff radius. */
+ float h;
+
+ /*! Particle time bin */
+ timebin_t time_bin;
+
+ struct {
+
+ /* Number of neighbours. */
+ float wcount;
+
+ /* Number of neighbours spatial derivative. */
+ float wcount_dh;
+
+ } density;
+
+ /*! Union for the formation time and formation scale factor */
+ union {
+
+ /*! Formation time */
+ float formation_time;
+
+ /*! Formation scale factor */
+ float formation_scale_factor;
+ };
+
+ /*! Subgrid mass of the black hole */
+ float subgrid_mass;
+
+ /*! Total accreted mass of the black hole */
+ float total_accreted_mass;
+
+ /*! Energy reservoir for feedback */
+ float energy_reservoir;
+
+ /*! Instantaneous accretion rate */
+ float accretion_rate;
+
+ /*! Density of the gas surrounding the black hole. */
+ float rho_gas;
+
+ /*! Smoothed sound speed of the gas surrounding the black hole. */
+ float sound_speed_gas;
+
+ /*! Smoothed velocity (peculiar) of the gas surrounding the black hole */
+ float velocity_gas[3];
+
+ /*! Total mass of the gas neighbours. */
+ float ngb_mass;
+
+ /*! Integer number of neighbours */
+ int num_ngbs;
+
+ /*! Properties used in the feedback loop to distribute to gas neighbours. */
+ struct {
+
+ /*! Probability of heating neighbouring gas particles for AGN feedback */
+ float AGN_heating_probability;
+
+ /*! Change in energy from SNII feedback energy injection */
+ float AGN_delta_u;
+
+ } to_distribute;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+ /*! Number of interactions in the density SELF and PAIR */
+ int num_ngb_density;
+
+ /*! List of interacting particles in the density SELF and PAIR */
+ long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_BLACK_HOLES];
+
+ /*! Number of interactions in the force SELF and PAIR */
+ int num_ngb_force;
+
+ /*! List of interacting particles in the force SELF and PAIR */
+ long long ids_ngbs_force[MAX_NUM_OF_NEIGHBOURS_BLACK_HOLES];
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_EAGLE_BLACK_HOLE_PART_H */
diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h
new file mode 100644
index 0000000000000000000000000000000000000000..01e420509d6eaee7dd07e5af23fc845acbcd330d
--- /dev/null
+++ b/src/black_holes/EAGLE/black_holes_properties.h
@@ -0,0 +1,181 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 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_EAGLE_BLACK_HOLES_PROPERTIES_H
+#define SWIFT_EAGLE_BLACK_HOLES_PROPERTIES_H
+
+#include "chemistry.h"
+#include "hydro_properties.h"
+
+/**
+ * @brief Properties of black holes and AGN feedback in the EAGEL model.
+ */
+struct black_holes_props {
+
+ /* ----- Basic neighbour search properties ------ */
+
+ /*! Resolution parameter */
+ float eta_neighbours;
+
+ /*! Target weightd number of neighbours (for info only)*/
+ float target_neighbours;
+
+ /*! Smoothing length tolerance */
+ float h_tolerance;
+
+ /*! Tolerance on neighbour number (for info only)*/
+ float delta_neighbours;
+
+ /*! Maximal number of iterations to converge h */
+ int max_smoothing_iterations;
+
+ /*! Maximal change of h over one time-step */
+ float log_max_h_change;
+
+ /* ----- Initialisation properties ------ */
+
+ /* ----- Properties of the accretion model ------ */
+
+ /*! Maximal fraction of the Eddington rate allowed. */
+ float f_Edd;
+
+ /*! Radiative efficiency of the black holes. */
+ float epsilon_r;
+
+ /*! Feedback coupling efficiency of the black holes. */
+ float epsilon_f;
+
+ /* ---- Properties of the feedback model ------- */
+
+ /*! Temperature increase induced by AGN feedback (Kelvin) */
+ float AGN_delta_T_desired;
+
+ /*! Number of gas neighbours to heat in a feedback event */
+ float num_ngbs_to_heat;
+
+ /* ---- Common conversion factors --------------- */
+
+ /*! Conversion factor from temperature to internal energy */
+ float temp_to_u_factor;
+};
+
+/**
+ * @brief Initialise the black hole properties from the parameter file.
+ *
+ * For the basic black holes neighbour finding properties we use the
+ * defaults from the hydro scheme if the users did not provide specific
+ * values.
+ *
+ * @param bp The #black_holes_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ * @param hydro_props The already read-in properties of the hydro scheme.
+ * @param cosmo The cosmological model.
+ */
+INLINE static void black_holes_props_init(struct black_holes_props *bp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *hydro_props,
+ const struct cosmology *cosmo) {
+
+ /* Read in the basic neighbour search properties or default to the hydro
+ ones if the user did not provide any different values */
+
+ /* Kernel properties */
+ bp->eta_neighbours = parser_get_opt_param_float(
+ params, "BlackHoles:resolution_eta", hydro_props->eta_neighbours);
+
+ /* Tolerance for the smoothing length Newton-Raphson scheme */
+ bp->h_tolerance = parser_get_opt_param_float(params, "BlackHoles:h_tolerance",
+ hydro_props->h_tolerance);
+
+ /* Get derived properties */
+ bp->target_neighbours = pow_dimension(bp->eta_neighbours) * kernel_norm;
+ const float delta_eta = bp->eta_neighbours * (1.f + bp->h_tolerance);
+ bp->delta_neighbours =
+ (pow_dimension(delta_eta) - pow_dimension(bp->eta_neighbours)) *
+ kernel_norm;
+
+ /* Number of iterations to converge h */
+ bp->max_smoothing_iterations =
+ parser_get_opt_param_int(params, "BlackHoles:max_ghost_iterations",
+ hydro_props->max_smoothing_iterations);
+
+ /* Time integration properties */
+ const float max_volume_change =
+ parser_get_opt_param_float(params, "BlackHoles:max_volume_change", -1);
+ if (max_volume_change == -1)
+ bp->log_max_h_change = hydro_props->log_max_h_change;
+ else
+ bp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+
+ /* Accretion parameters ---------------------------------- */
+
+ bp->f_Edd = parser_get_param_float(params, "EAGLEAGN:max_eddington_fraction");
+ bp->epsilon_r =
+ parser_get_param_float(params, "EAGLEAGN:radiative_efficiency");
+ bp->epsilon_f =
+ parser_get_param_float(params, "EAGLEAGN:coupling_efficiency");
+
+ /* Feedback parameters ---------------------------------- */
+
+ bp->AGN_delta_T_desired =
+ parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_K");
+
+ bp->num_ngbs_to_heat =
+ parser_get_param_float(params, "EAGLEAGN:AGN_num_ngb_to_heat");
+
+ /* Common conversion factors ----------------------------- */
+
+ /* Calculate temperature to internal energy conversion factor (all internal
+ * units) */
+ const double k_B = phys_const->const_boltzmann_k;
+ const double m_p = phys_const->const_proton_mass;
+ const double mu = hydro_props->mu_ionised;
+ bp->temp_to_u_factor = k_B / (mu * hydro_gamma_minus_one * m_p);
+}
+
+/**
+ * @brief Write a black_holes_props struct to the given FILE as a stream of
+ * bytes.
+ *
+ * @param props the black hole properties struct
+ * @param stream the file stream
+ */
+INLINE static void black_holes_struct_dump(
+ const struct black_holes_props *props, FILE *stream) {
+ restart_write_blocks((void *)props, sizeof(struct black_holes_props), 1,
+ stream, "black_holes props", "black holes props");
+}
+
+/**
+ * @brief Restore a black_holes_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param props the black hole properties struct
+ * @param stream the file stream
+ */
+INLINE static void black_holes_struct_restore(
+ const struct black_holes_props *props, FILE *stream) {
+ restart_read_blocks((void *)props, sizeof(struct black_holes_props), 1,
+ stream, NULL, "black holes props");
+}
+
+#endif /* SWIFT_EAGLE_BLACK_HOLES_PROPERTIES_H */
diff --git a/src/black_holes_io.h b/src/black_holes_io.h
index 05e66645b182bbd26907b37d5c5560aa5e8b27a7..bc3b1be3be44588794ababfc9500be7812c14f5e 100644
--- a/src/black_holes_io.h
+++ b/src/black_holes_io.h
@@ -24,6 +24,8 @@
/* Load the correct star type */
#if defined(BLACK_HOLES_NONE)
#include "./black_holes/Default/black_holes_io.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes_io.h"
#else
#error "Invalid choice of star model"
#endif
diff --git a/src/black_holes_properties.h b/src/black_holes_properties.h
new file mode 100644
index 0000000000000000000000000000000000000000..3226b5391a979d27ab125c88d50acb3d59e1edf0
--- /dev/null
+++ b/src/black_holes_properties.h
@@ -0,0 +1,34 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 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_BLACK_HOLES_PROPERTIES_H
+#define SWIFT_BLACK_HOLES_PROPERTIES_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct black_holes model */
+#if defined(BLACK_HOLES_NONE)
+#include "./black_holes/Default/black_holes_properties.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes_properties.h"
+#else
+#error "Invalid choice of black hole model"
+#endif
+
+#endif /* SWIFT_BLACK_HOLES_PROPERTIES_H */
diff --git a/src/cell.c b/src/cell.c
index 86b54823a1860416968ab236e36d4ee43566d6af..4c6bc5b6fdb168c315f777600ea47ae19d021bb1 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -497,12 +497,15 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
/* Start by packing the data of the current cell. */
pc->hydro.h_max = c->hydro.h_max;
pc->stars.h_max = c->stars.h_max;
+ pc->black_holes.h_max = c->black_holes.h_max;
pc->hydro.ti_end_min = c->hydro.ti_end_min;
pc->hydro.ti_end_max = c->hydro.ti_end_max;
pc->grav.ti_end_min = c->grav.ti_end_min;
pc->grav.ti_end_max = c->grav.ti_end_max;
pc->stars.ti_end_min = c->stars.ti_end_min;
pc->stars.ti_end_max = c->stars.ti_end_max;
+ pc->black_holes.ti_end_min = c->black_holes.ti_end_min;
+ pc->black_holes.ti_end_max = c->black_holes.ti_end_max;
pc->hydro.ti_old_part = c->hydro.ti_old_part;
pc->grav.ti_old_part = c->grav.ti_old_part;
pc->grav.ti_old_multipole = c->grav.ti_old_multipole;
@@ -510,6 +513,7 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
pc->hydro.count = c->hydro.count;
pc->grav.count = c->grav.count;
pc->stars.count = c->stars.count;
+ pc->black_holes.count = c->black_holes.count;
pc->maxdepth = c->maxdepth;
/* Copy the Multipole related information */
@@ -602,19 +606,24 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
/* Unpack the current pcell. */
c->hydro.h_max = pc->hydro.h_max;
c->stars.h_max = pc->stars.h_max;
+ c->black_holes.h_max = pc->black_holes.h_max;
c->hydro.ti_end_min = pc->hydro.ti_end_min;
c->hydro.ti_end_max = pc->hydro.ti_end_max;
c->grav.ti_end_min = pc->grav.ti_end_min;
c->grav.ti_end_max = pc->grav.ti_end_max;
c->stars.ti_end_min = pc->stars.ti_end_min;
c->stars.ti_end_max = pc->stars.ti_end_max;
+ c->black_holes.ti_end_min = pc->black_holes.ti_end_min;
+ c->black_holes.ti_end_max = pc->black_holes.ti_end_max;
c->hydro.ti_old_part = pc->hydro.ti_old_part;
c->grav.ti_old_part = pc->grav.ti_old_part;
c->grav.ti_old_multipole = pc->grav.ti_old_multipole;
c->stars.ti_old_part = pc->stars.ti_old_part;
+ c->black_holes.ti_old_part = pc->black_holes.ti_old_part;
c->hydro.count = pc->hydro.count;
c->grav.count = pc->grav.count;
c->stars.count = pc->stars.count;
+ c->black_holes.count = pc->black_holes.count;
c->maxdepth = pc->maxdepth;
#ifdef SWIFT_DEBUG_CHECKS
@@ -664,6 +673,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
temp->hydro.dx_max_sort = 0.f;
temp->stars.dx_max_part = 0.f;
temp->stars.dx_max_sort = 0.f;
+ temp->black_holes.dx_max_part = 0.f;
temp->nodeID = c->nodeID;
temp->parent = c;
c->progeny[k] = temp;
@@ -1786,6 +1796,8 @@ void cell_clean_links(struct cell *c, void *data) {
c->grav.mm = NULL;
c->stars.density = NULL;
c->stars.feedback = NULL;
+ c->black_holes.density = NULL;
+ c->black_holes.feedback = NULL;
}
/**
@@ -2311,6 +2323,44 @@ void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
}
}
+/**
+ * @brief Activate the #bpart drifts on the given cell.
+ */
+void cell_activate_drift_bpart(struct cell *c, struct scheduler *s) {
+
+ /* If this cell is already marked for drift, quit early. */
+ if (cell_get_flag(c, cell_flag_do_bh_drift)) return;
+
+ /* Mark this cell for drifting. */
+ cell_set_flag(c, cell_flag_do_bh_drift);
+
+ /* Set the do_black_holes_sub_drifts all the way up and activate the super
+ drift if this has not yet been done. */
+ if (c == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->black_holes.drift == NULL)
+ error("Trying to activate un-existing c->black_holes.drift");
+#endif
+ scheduler_activate(s, c->black_holes.drift);
+ } else {
+ for (struct cell *parent = c->parent;
+ parent != NULL && !cell_get_flag(parent, cell_flag_do_bh_sub_drift);
+ parent = parent->parent) {
+ /* Mark this cell for drifting */
+ cell_set_flag(parent, cell_flag_do_bh_sub_drift);
+
+ if (parent == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (parent->black_holes.drift == NULL)
+ error("Trying to activate un-existing parent->black_holes.drift");
+#endif
+ scheduler_activate(s, parent->black_holes.drift);
+ break;
+ }
+ }
+ }
+}
+
/**
* @brief Activate the drifts on the given cell.
*/
@@ -2412,7 +2462,6 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
#endif
scheduler_activate(s, c->stars.sorts);
if (c->nodeID == engine_rank) {
- // MATTHIEU: to do: do we actually need both drifts here?
cell_activate_drift_spart(c, s);
}
} else {
@@ -2662,6 +2711,100 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
} /* Otherwise, pair interation */
}
+/**
+ * @brief Traverse a sub-cell task and activate the black_holes drift tasks that
+ * are required by a black_holes task
+ *
+ * @param ci The first #cell we recurse in.
+ * @param cj The second #cell we recurse in.
+ * @param s The task #scheduler.
+ */
+void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
+ struct scheduler *s) {
+ const struct engine *e = s->space->e;
+
+ /* Store the current dx_max and h_max values. */
+ ci->black_holes.dx_max_part_old = ci->black_holes.dx_max_part;
+ ci->black_holes.h_max_old = ci->black_holes.h_max;
+ ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
+ ci->hydro.h_max_old = ci->hydro.h_max;
+
+ if (cj != NULL) {
+ cj->black_holes.dx_max_part_old = cj->black_holes.dx_max_part;
+ cj->black_holes.h_max_old = cj->black_holes.h_max;
+ cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
+ cj->hydro.h_max_old = cj->hydro.h_max;
+ }
+
+ /* Self interaction? */
+ if (cj == NULL) {
+ /* Do anything? */
+ if (!cell_is_active_black_holes(ci, e) || ci->hydro.count == 0 ||
+ ci->black_holes.count == 0)
+ return;
+
+ /* Recurse? */
+ if (cell_can_recurse_in_self_black_holes_task(ci)) {
+ /* Loop over all progenies and pairs of progenies */
+ for (int j = 0; j < 8; j++) {
+ if (ci->progeny[j] != NULL) {
+ cell_activate_subcell_black_holes_tasks(ci->progeny[j], NULL, s);
+ for (int k = j + 1; k < 8; k++)
+ if (ci->progeny[k] != NULL)
+ cell_activate_subcell_black_holes_tasks(ci->progeny[j],
+ ci->progeny[k], s);
+ }
+ }
+ } else {
+ /* We have reached the bottom of the tree: activate drift */
+ cell_activate_drift_bpart(ci, s);
+ cell_activate_drift_part(ci, s);
+ }
+ }
+
+ /* Otherwise, pair interation */
+ else {
+ /* Should we even bother? */
+ if (!cell_is_active_black_holes(ci, e) &&
+ !cell_is_active_black_holes(cj, e))
+ return;
+
+ /* Get the orientation of the pair. */
+ double shift[3];
+ int sid = space_getsid(s->space, &ci, &cj, shift);
+
+ /* recurse? */
+ if (cell_can_recurse_in_pair_black_holes_task(ci, cj) &&
+ cell_can_recurse_in_pair_black_holes_task(cj, ci)) {
+ struct cell_split_pair *csp = &cell_split_pairs[sid];
+ for (int k = 0; k < csp->count; k++) {
+ const int pid = csp->pairs[k].pid;
+ const int pjd = csp->pairs[k].pjd;
+ if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL)
+ cell_activate_subcell_black_holes_tasks(ci->progeny[pid],
+ cj->progeny[pjd], s);
+ }
+ }
+
+ /* Otherwise, activate the drifts. */
+ else {
+ if (cell_is_active_black_holes(ci, e)) {
+
+ /* Activate the drifts if the cells are local. */
+ if (ci->nodeID == engine_rank) cell_activate_drift_bpart(ci, s);
+ if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
+ }
+
+ if (cell_is_active_black_holes(cj, e)) {
+
+ /* Activate the drifts if the cells are local. */
+ if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
+ if (cj->nodeID == engine_rank) cell_activate_drift_bpart(cj, s);
+ }
+ }
+ } /* Otherwise, pair interation */
+}
+
/**
* @brief Traverse a sub-cell task and activate the gravity drift tasks that
* are required by a self gravity task.
@@ -3394,6 +3537,205 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
return rebuild;
}
+/**
+ * @brief Un-skips all the black hole tasks associated with a given cell and
+ * checks if the space needs to be rebuilt.
+ *
+ * @param c the #cell.
+ * @param s the #scheduler.
+ *
+ * @return 1 If the space needs rebuilding. 0 otherwise.
+ */
+int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
+ struct engine *e = s->space->e;
+ const int nodeID = e->nodeID;
+ int rebuild = 0;
+
+ if (c->black_holes.drift != NULL && cell_is_active_black_holes(c, e)) {
+ cell_activate_drift_bpart(c, s);
+ }
+
+ /* Un-skip the density tasks involved with this cell. */
+ for (struct link *l = c->black_holes.density; l != NULL; l = l->next) {
+ struct task *t = l->t;
+ struct cell *ci = t->ci;
+ struct cell *cj = t->cj;
+ const int ci_active = cell_is_active_black_holes(ci, e);
+ const int cj_active = (cj != NULL) ? cell_is_active_black_holes(cj, e) : 0;
+#ifdef WITH_MPI
+ const int ci_nodeID = ci->nodeID;
+ const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
+#else
+ const int ci_nodeID = nodeID;
+ const int cj_nodeID = nodeID;
+#endif
+
+ /* Activate the drifts */
+ if (t->type == task_type_self && ci_active) {
+ cell_activate_drift_part(ci, s);
+ cell_activate_drift_bpart(ci, s);
+ }
+
+ /* Only activate tasks that involve a local active cell. */
+ if ((ci_active || cj_active) &&
+ (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
+ scheduler_activate(s, t);
+
+ if (t->type == task_type_pair) {
+ /* Do ci */
+ if (ci_active) {
+
+ /* Activate the drift tasks. */
+ if (ci_nodeID == nodeID) cell_activate_drift_bpart(ci, s);
+ if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
+ }
+
+ /* Do cj */
+ if (cj_active) {
+
+ /* Activate the drift tasks. */
+ if (cj_nodeID == nodeID) cell_activate_drift_bpart(cj, s);
+ if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
+ }
+ }
+
+ else if (t->type == task_type_sub_pair || t->type == task_type_sub_self) {
+ cell_activate_subcell_black_holes_tasks(ci, cj, s);
+ }
+ }
+
+ /* Only interested in pair interactions as of here. */
+ if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+ /* Check whether there was too much particle motion, i.e. the
+ cell neighbour conditions were violated. */
+ if (cell_need_rebuild_for_black_holes_pair(ci, cj)) rebuild = 1;
+ if (cell_need_rebuild_for_black_holes_pair(cj, ci)) rebuild = 1;
+
+#ifdef WITH_MPI
+ /* Activate the send/recv tasks. */
+ if (ci_nodeID != nodeID) {
+ if (cj_active) {
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
+
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart,
+ ci_nodeID);
+ cell_activate_drift_bpart(cj, s);
+
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_tend_bpart,
+ ci_nodeID);
+ }
+
+ if (ci_active) {
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart);
+
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_tend_bpart);
+
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_xv, ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_rho, ci_nodeID);
+
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(cj, s);
+ }
+
+ } else if (cj_nodeID != nodeID) {
+ /* If the local cell is active, receive data from the foreign cell. */
+ if (ci_active) {
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
+
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart,
+ cj_nodeID);
+ cell_activate_drift_bpart(ci, s);
+
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_bpart,
+ cj_nodeID);
+ }
+
+ if (cj_active) {
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart);
+
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_tend_bpart);
+
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_xv, cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_rho, cj_nodeID);
+
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(ci, s);
+ }
+ }
+#endif
+ }
+ }
+
+ /* Un-skip the feedback tasks involved with this cell. */
+ for (struct link *l = c->black_holes.feedback; l != NULL; l = l->next) {
+ struct task *t = l->t;
+ struct cell *ci = t->ci;
+ struct cell *cj = t->cj;
+ const int ci_active = cell_is_active_black_holes(ci, e);
+ const int cj_active = (cj != NULL) ? cell_is_active_black_holes(cj, e) : 0;
+#ifdef WITH_MPI
+ const int ci_nodeID = ci->nodeID;
+ const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
+#else
+ const int ci_nodeID = nodeID;
+ const int cj_nodeID = nodeID;
+#endif
+
+ if (t->type == task_type_self && ci_active) {
+ scheduler_activate(s, t);
+ }
+
+ else if (t->type == task_type_sub_self && ci_active) {
+ scheduler_activate(s, t);
+ }
+
+ else if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+ /* We only want to activate the task if the cell is active and is
+ going to update some gas on the *local* node */
+ if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
+ (ci_active || cj_active)) {
+ scheduler_activate(s, t);
+
+ } else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) && (cj_active)) {
+ scheduler_activate(s, t);
+
+ } else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) && (ci_active)) {
+ scheduler_activate(s, t);
+ }
+ }
+
+ /* Nothing more to do here, all drifts activated above */
+ }
+
+ /* Unskip all the other task types. */
+ if (c->nodeID == nodeID && cell_is_active_black_holes(c, e)) {
+ if (c->black_holes.ghost != NULL)
+ scheduler_activate(s, c->black_holes.ghost);
+ if (c->black_holes.black_holes_in != NULL)
+ scheduler_activate(s, c->black_holes.black_holes_in);
+ if (c->black_holes.black_holes_out != NULL)
+ scheduler_activate(s, c->black_holes.black_holes_out);
+ if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
+ if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
+ if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+ if (c->logger != NULL) scheduler_activate(s, c->logger);
+ }
+
+ return rebuild;
+}
+
/**
* @brief Set the super-cell pointers for all cells in a hierarchy.
*
@@ -4236,7 +4578,8 @@ void cell_check_timesteps(struct cell *c) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->hydro.ti_end_min == 0 && c->grav.ti_end_min == 0 &&
- c->stars.ti_end_min == 0 && c->nr_tasks > 0)
+ c->stars.ti_end_min == 0 && c->black_holes.ti_end_min == 0 &&
+ c->nr_tasks > 0)
error("Cell without assigned time-step");
if (c->split) {
diff --git a/src/cell.h b/src/cell.h
index 0c2a1bb5be414530e4ac42af76229c49843ee544..123b7fc5fda2f4b9d6b94c50eeaee991c849a23d 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -172,6 +172,26 @@ struct pcell {
} stars;
+ /*! Black hole variables */
+ struct {
+
+ /*! Number of #spart in this cell. */
+ int count;
+
+ /*! Maximal smoothing length. */
+ double h_max;
+
+ /*! Minimal integer end-of-timestep in this cell for black hole tasks */
+ integertime_t ti_end_min;
+
+ /*! Maximal integer end-of-timestep in this cell for black hole tasks */
+ integertime_t ti_end_max;
+
+ /*! Integer time of the last drift of the #spart in this cell */
+ integertime_t ti_old_part;
+
+ } black_holes;
+
/*! Maximal depth in that part of the tree */
int maxdepth;
@@ -606,6 +626,22 @@ struct cell {
/*! The drift task for bparts */
struct task *drift;
+ /*! Implicit tasks marking the entry of the BH physics block of tasks
+ */
+ struct task *black_holes_in;
+
+ /*! Implicit tasks marking the exit of the BH physics block of tasks */
+ struct task *black_holes_out;
+
+ /*! The star ghost task itself */
+ struct task *ghost;
+
+ /*! Linked list of the tasks computing this cell's star density. */
+ struct link *density;
+
+ /*! Linked list of the tasks computing this cell's star feedback. */
+ struct link *feedback;
+
/*! Max smoothing length in this cell. */
double h_max;
@@ -627,6 +663,9 @@ struct cell {
/*! Maximum part movement in this cell since last construction. */
float dx_max_part;
+ /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+ float dx_max_part_old;
+
/*! Maximum end of (integer) time step in this cell for black tasks. */
integertime_t ti_end_min;
@@ -784,6 +823,7 @@ void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s);
int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s);
+int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s);
int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s);
void cell_drift_part(struct cell *c, const struct engine *e, int force);
void cell_drift_gpart(struct cell *c, const struct engine *e, int force);
@@ -799,12 +839,15 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s);
void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s);
+void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
+ struct scheduler *s);
void cell_activate_subcell_external_grav_tasks(struct cell *ci,
struct scheduler *s);
void cell_activate_super_spart_drifts(struct cell *c, struct scheduler *s);
void cell_activate_drift_part(struct cell *c, struct scheduler *s);
void cell_activate_drift_gpart(struct cell *c, struct scheduler *s);
void cell_activate_drift_spart(struct cell *c, struct scheduler *s);
+void cell_activate_drift_bpart(struct cell *c, struct scheduler *s);
void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s);
void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s);
void cell_activate_limiter(struct cell *c, struct scheduler *s);
@@ -969,6 +1012,43 @@ cell_can_recurse_in_self_stars_task(const struct cell *c) {
(kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
}
+/**
+ * @brief Can a sub-pair black hole task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param ci The #cell with black holes.
+ * @param cj The #cell with hydro parts.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_pair_black_holes_task(const struct cell *ci,
+ const struct cell *cj) {
+
+ /* Is the cell split ? */
+ /* If so, is the cut-off radius plus the max distance the parts have moved */
+ /* smaller than the sub-cell sizes ? */
+ /* Note: We use the _old values as these might have been updated by a drift */
+ return ci->split && cj->split &&
+ ((kernel_gamma * ci->black_holes.h_max_old +
+ ci->black_holes.dx_max_part_old) < 0.5f * ci->dmin) &&
+ ((kernel_gamma * cj->hydro.h_max_old + cj->hydro.dx_max_part_old) <
+ 0.5f * cj->dmin);
+}
+
+/**
+ * @brief Can a sub-self black hole task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param c The #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_self_black_holes_task(const struct cell *c) {
+
+ /* Is the cell split and not smaller than the smoothing length? */
+ return c->split &&
+ (kernel_gamma * c->black_holes.h_max_old < 0.5f * c->dmin) &&
+ (kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
+}
+
/**
* @brief Can a pair hydro task associated with a cell be split into smaller
* sub-tasks.
@@ -985,7 +1065,8 @@ __attribute__((always_inline)) INLINE static int cell_can_split_pair_hydro_task(
/* into account any part motion (i.e. dx_max == 0 here) */
return c->split &&
(space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin) &&
- (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
+ (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin) &&
+ (space_stretch * kernel_gamma * c->black_holes.h_max < 0.5f * c->dmin);
}
/**
@@ -1004,7 +1085,8 @@ __attribute__((always_inline)) INLINE static int cell_can_split_self_hydro_task(
/* tasks will be created. So no need to check for h_max */
return c->split &&
(space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin) &&
- (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
+ (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin) &&
+ (space_stretch * kernel_gamma * c->black_holes.h_max < 0.5f * c->dmin);
}
/**
@@ -1054,6 +1136,7 @@ cell_need_rebuild_for_hydro_pair(const struct cell *ci, const struct cell *cj) {
}
return 0;
}
+
/**
* @brief Have star particles in a pair of cells moved too much and require a
* rebuild?
@@ -1075,6 +1158,28 @@ cell_need_rebuild_for_stars_pair(const struct cell *ci, const struct cell *cj) {
return 0;
}
+/**
+ * @brief Have star particles in a pair of cells moved too much and require a
+ * rebuild?
+ *
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_need_rebuild_for_black_holes_pair(const struct cell *ci,
+ const struct cell *cj) {
+
+ /* Is the cut-off radius plus the max distance the parts in both cells have */
+ /* moved larger than the cell size ? */
+ /* Note ci->dmin == cj->dmin */
+ if (kernel_gamma * max(ci->black_holes.h_max, cj->hydro.h_max) +
+ ci->black_holes.dx_max_part + cj->hydro.dx_max_part >
+ cj->dmin) {
+ return 1;
+ }
+ return 0;
+}
+
/**
* @brief Add a unique tag to a cell, mostly for MPI communications.
*
diff --git a/src/engine.c b/src/engine.c
index afa859e9aeff53e7d03eae5eeaa24a9d017766d1..a82e8db6d6b73ef95c5705b8b6a4adb3b74cec05 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2873,13 +2873,10 @@ void engine_collect_end_of_step_recurse_stars(struct cell *c,
void engine_collect_end_of_step_recurse_black_holes(struct cell *c,
const struct engine *e) {
-/* Skip super-cells (Their values are already set) */
-#ifdef WITH_MPI
- // MATTHIEU
- if (c->timestep != NULL)
- return; // || c->mpi.black_holes.recv_ti != NULL) return;
-#else
+ /* Skip super-cells (Their values are already set) */
if (c->timestep != NULL) return;
+#ifdef WITH_MPI
+ if (cell_get_recv(c, task_subtype_tend_bpart) != NULL) return;
#endif /* WITH_MPI */
#ifdef SWIFT_DEBUG_CHECKS
@@ -3305,8 +3302,9 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Skip everything that updates the particles */
if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
- t->type == task_type_drift_spart || t->type == task_type_kick1 ||
- t->type == task_type_kick2 || t->type == task_type_timestep ||
+ t->type == task_type_drift_spart || t->type == task_type_drift_bpart ||
+ t->type == task_type_kick1 || t->type == task_type_kick2 ||
+ t->type == task_type_timestep ||
t->type == task_type_timestep_limiter ||
t->subtype == task_subtype_force ||
t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
@@ -3320,9 +3318,11 @@ void engine_skip_force_and_kick(struct engine *e) {
t->type == task_type_extra_ghost ||
t->subtype == task_subtype_gradient ||
t->subtype == task_subtype_stars_feedback ||
+ t->subtype == task_subtype_bh_feedback ||
t->subtype == task_subtype_tend_part ||
t->subtype == task_subtype_tend_gpart ||
t->subtype == task_subtype_tend_spart ||
+ t->subtype == task_subtype_tend_bpart ||
t->subtype == task_subtype_rho || t->subtype == task_subtype_gpart)
t->skip = 1;
}
@@ -4115,6 +4115,7 @@ void engine_unskip(struct engine *e) {
const int with_ext_grav = e->policy & engine_policy_external_gravity;
const int with_stars = e->policy & engine_policy_stars;
const int with_feedback = e->policy & engine_policy_feedback;
+ const int with_black_holes = e->policy & engine_policy_black_holes;
#ifdef WITH_PROFILER
static int count = 0;
@@ -4134,7 +4135,8 @@ void engine_unskip(struct engine *e) {
(with_ext_grav && c->nodeID == nodeID &&
cell_is_active_gravity(c, e)) ||
(with_feedback && cell_is_active_stars(c, e)) ||
- (with_stars && c->nodeID == nodeID && cell_is_active_stars(c, e))) {
+ (with_stars && c->nodeID == nodeID && cell_is_active_stars(c, e)) ||
+ (with_black_holes && cell_is_active_black_holes(c, e))) {
if (num_active_cells != k)
memswap(&local_cells[k], &local_cells[num_active_cells], sizeof(int));
@@ -4857,6 +4859,7 @@ void engine_unpin(void) {
* @param entropy_floor The #entropy_floor_properties for this run.
* @param gravity The #gravity_props used for this run.
* @param stars The #stars_props used for this run.
+ * @param black_holes The #black_holes_props used for this run.
* @param feedback The #feedback_props used for this run.
* @param mesh The #pm_mesh used for the long-range periodic forces.
* @param potential The properties of the external potential.
@@ -4872,6 +4875,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cosmology *cosmo, struct hydro_props *hydro,
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
+ const struct black_holes_props *black_holes,
const struct feedback_props *feedback, struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
@@ -4939,6 +4943,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->entropy_floor = entropy_floor;
e->gravity_properties = gravity;
e->stars_properties = stars;
+ e->black_holes_properties = black_holes;
e->mesh = mesh;
e->external_potential = potential;
e->cooling_func = cooling_func;
@@ -6083,6 +6088,7 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
cooling_struct_dump(e->cooling_func, stream);
starformation_struct_dump(e->star_formation, stream);
feedback_struct_dump(e->feedback_props, stream);
+ black_holes_struct_dump(e->black_holes_properties, stream);
chemistry_struct_dump(e->chemistry, stream);
parser_struct_dump(e->parameter_file, stream);
if (e->output_list_snapshots)
@@ -6190,6 +6196,11 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
feedback_struct_restore(feedback_properties, stream);
e->feedback_props = feedback_properties;
+ struct black_holes_props *black_holes_properties =
+ (struct black_holes_props *)malloc(sizeof(struct black_holes_props));
+ black_holes_struct_restore(black_holes_properties, stream);
+ e->black_holes_properties = black_holes_properties;
+
struct chemistry_global_data *chemistry =
(struct chemistry_global_data *)malloc(
sizeof(struct chemistry_global_data));
diff --git a/src/engine.h b/src/engine.h
index 3afa8c221d2d1690dec45d2017a18fc0d3aee8ca..1a653fdfd40c88775b7bbb67265922c9db68a5fa 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -34,6 +34,7 @@
/* Includes. */
#include "barrier.h"
+#include "black_holes_properties.h"
#include "chemistry_struct.h"
#include "clocks.h"
#include "collectgroup.h"
@@ -399,6 +400,9 @@ struct engine {
/* Properties of the star model */
const struct stars_props *stars_properties;
+ /* Properties of the black hole model */
+ const struct black_holes_props *black_holes_properties;
+
/* Properties of the self-gravity scheme */
struct gravity_props *gravity_properties;
@@ -473,6 +477,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cosmology *cosmo, struct hydro_props *hydro,
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
+ const struct black_holes_props *black_holes,
const struct feedback_props *feedback, struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index ca69599b8e64460484e174aca9fc786dd45d8349..201bc357d7e85c89712847804e0ba3d111c8661f 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -292,6 +292,75 @@ void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
#endif
}
+/**
+ * @brief Add send tasks for the black holes pairs to a hierarchy of cells.
+ *
+ * @param e The #engine.
+ * @param ci The sending #cell.
+ * @param cj Dummy cell containing the nodeID of the receiving node.
+ * @param t_feedback The send_feed #task, if it has already been created.
+ * @param t_ti The recv_ti_end #task, if it has already been created.
+ */
+void engine_addtasks_send_black_holes(struct engine *e, struct cell *ci,
+ struct cell *cj, struct task *t_feedback,
+ struct task *t_ti) {
+
+#ifdef WITH_MPI
+
+ struct link *l = NULL;
+ struct scheduler *s = &e->sched;
+ const int nodeID = cj->nodeID;
+
+ /* Check if any of the density tasks are for the target node. */
+ for (l = ci->black_holes.density; l != NULL; l = l->next)
+ if (l->t->ci->nodeID == nodeID ||
+ (l->t->cj != NULL && l->t->cj->nodeID == nodeID))
+ break;
+
+ /* If so, attach send tasks. */
+ if (l != NULL) {
+
+ if (t_feedback == NULL) {
+
+ /* Make sure this cell is tagged. */
+ cell_ensure_tagged(ci);
+
+ /* Create the tasks and their dependencies? */
+ t_feedback = scheduler_addtask(s, task_type_send, task_subtype_bpart,
+ ci->mpi.tag, 0, ci, cj);
+
+ t_ti = scheduler_addtask(s, task_type_send, task_subtype_tend_bpart,
+ ci->mpi.tag, 0, ci, cj);
+
+ /* The send_black_holes task should unlock the super_cell's kick task. */
+ scheduler_addunlock(s, t_feedback,
+ ci->hydro.super->black_holes.black_holes_out);
+
+ /* Ghost before you send */
+ scheduler_addunlock(s, ci->hydro.super->black_holes.ghost, t_feedback);
+
+ /* Drift before you send */
+ scheduler_addunlock(s, ci->hydro.super->black_holes.drift, t_feedback);
+
+ scheduler_addunlock(s, ci->super->timestep, t_ti);
+ }
+
+ engine_addlink(e, &ci->mpi.send, t_feedback);
+ engine_addlink(e, &ci->mpi.send, t_ti);
+ }
+
+ /* Recurse? */
+ if (ci->split)
+ for (int k = 0; k < 8; k++)
+ if (ci->progeny[k] != NULL)
+ engine_addtasks_send_black_holes(e, ci->progeny[k], cj, t_feedback,
+ t_ti);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
/**
* @brief Add recv tasks for hydro pairs to a hierarchy of cells.
*
@@ -445,6 +514,66 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
#endif
}
+/**
+ * @brief Add recv tasks for black_holes pairs to a hierarchy of cells.
+ *
+ * @param e The #engine.
+ * @param c The foreign #cell.
+ * @param t_feedback The recv_feed #task, if it has already been created.
+ * @param t_ti The recv_ti_end #task, if it has already been created.
+ */
+void engine_addtasks_recv_black_holes(struct engine *e, struct cell *c,
+ struct task *t_feedback,
+ struct task *t_ti) {
+
+#ifdef WITH_MPI
+ struct scheduler *s = &e->sched;
+
+ /* Have we reached a level where there are any black_holes tasks ? */
+ if (t_feedback == NULL && c->black_holes.density != NULL) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure this cell has a valid tag. */
+ if (c->mpi.tag < 0) error("Trying to receive from untagged cell.");
+#endif // SWIFT_DEBUG_CHECKS
+
+ /* Create the tasks. */
+ t_feedback = scheduler_addtask(s, task_type_recv, task_subtype_bpart,
+ c->mpi.tag, 0, c, NULL);
+
+ t_ti = scheduler_addtask(s, task_type_recv, task_subtype_tend_bpart,
+ c->mpi.tag, 0, c, NULL);
+ }
+
+ if (t_feedback != NULL) {
+ engine_addlink(e, &c->mpi.recv, t_feedback);
+ engine_addlink(e, &c->mpi.recv, t_ti);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID == e->nodeID) error("Local cell!");
+#endif
+
+ for (struct link *l = c->black_holes.density; l != NULL; l = l->next) {
+ scheduler_addunlock(s, l->t, t_feedback);
+ }
+
+ for (struct link *l = c->black_holes.feedback; l != NULL; l = l->next) {
+ scheduler_addunlock(s, t_feedback, l->t);
+ scheduler_addunlock(s, l->t, t_ti);
+ }
+ }
+
+ /* Recurse? */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_addtasks_recv_black_holes(e, c->progeny[k], t_feedback, t_ti);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
/**
* @brief Add recv tasks for gravity pairs to a hierarchy of cells.
*
@@ -729,6 +858,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
const int with_feedback = (e->policy & engine_policy_feedback);
const int with_cooling = (e->policy & engine_policy_cooling);
const int with_star_formation = (e->policy & engine_policy_star_formation);
+ const int with_black_holes = (e->policy & engine_policy_black_holes);
/* Are we in a super-cell ? */
if (c->hydro.super == c) {
@@ -775,6 +905,13 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
scheduler_addunlock(s, c->stars.drift, c->super->kick2);
}
+ /* Black holes */
+ if (with_black_holes) {
+ c->black_holes.drift = scheduler_addtask(
+ s, task_type_drift_bpart, task_subtype_none, 0, 0, c, NULL);
+ scheduler_addunlock(s, c->black_holes.drift, c->super->kick2);
+ }
+
/* Subgrid tasks: cooling */
if (with_cooling) {
@@ -810,6 +947,25 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
c->stars.stars_in);
}
}
+
+ /* Subgrid tasks: black hole feedback */
+ if (with_black_holes) {
+
+ c->black_holes.black_holes_in =
+ scheduler_addtask(s, task_type_bh_in, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+
+ c->black_holes.black_holes_out =
+ scheduler_addtask(s, task_type_bh_out, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+
+ c->black_holes.ghost = scheduler_addtask(
+ s, task_type_bh_ghost, task_subtype_none, 0, 0, c, NULL);
+
+ scheduler_addunlock(s, c->super->kick2, c->black_holes.black_holes_in);
+ scheduler_addunlock(s, c->black_holes.black_holes_out,
+ c->super->timestep);
+ }
}
} else { /* We are above the super-cell so need to go deeper */
@@ -1085,6 +1241,10 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->stars.density, t);
} else if (t->subtype == task_subtype_stars_feedback) {
engine_addlink(e, &ci->stars.feedback, t);
+ } else if (t->subtype == task_subtype_bh_density) {
+ engine_addlink(e, &ci->black_holes.density, t);
+ } else if (t->subtype == task_subtype_bh_feedback) {
+ engine_addlink(e, &ci->black_holes.feedback, t);
}
/* Link pair tasks to cells. */
@@ -1104,6 +1264,12 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
} else if (t->subtype == task_subtype_stars_feedback) {
engine_addlink(e, &ci->stars.feedback, t);
engine_addlink(e, &cj->stars.feedback, t);
+ } else if (t->subtype == task_subtype_bh_density) {
+ engine_addlink(e, &ci->black_holes.density, t);
+ engine_addlink(e, &cj->black_holes.density, t);
+ } else if (t->subtype == task_subtype_bh_feedback) {
+ engine_addlink(e, &ci->black_holes.feedback, t);
+ engine_addlink(e, &cj->black_holes.feedback, t);
}
#ifdef SWIFT_DEBUG_CHECKS
else if (t_subtype == task_subtype_external_grav) {
@@ -1125,6 +1291,10 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->stars.density, t);
} else if (t->subtype == task_subtype_stars_feedback) {
engine_addlink(e, &ci->stars.feedback, t);
+ } else if (t->subtype == task_subtype_bh_density) {
+ engine_addlink(e, &ci->black_holes.density, t);
+ } else if (t->subtype == task_subtype_bh_feedback) {
+ engine_addlink(e, &ci->black_holes.feedback, t);
}
/* Link sub-pair tasks to cells. */
@@ -1144,6 +1314,12 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
} else if (t->subtype == task_subtype_stars_feedback) {
engine_addlink(e, &ci->stars.feedback, t);
engine_addlink(e, &cj->stars.feedback, t);
+ } else if (t->subtype == task_subtype_bh_density) {
+ engine_addlink(e, &ci->black_holes.density, t);
+ engine_addlink(e, &cj->black_holes.density, t);
+ } else if (t->subtype == task_subtype_bh_feedback) {
+ engine_addlink(e, &ci->black_holes.feedback, t);
+ engine_addlink(e, &cj->black_holes.feedback, t);
}
#ifdef SWIFT_DEBUG_CHECKS
else if (t_subtype == task_subtype_external_grav) {
@@ -1398,6 +1574,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
const int with_cooling = (e->policy & engine_policy_cooling);
const int with_limiter = (e->policy & engine_policy_limiter);
const int with_feedback = (e->policy & engine_policy_feedback);
+ const int with_black_holes = (e->policy & engine_policy_black_holes);
#ifdef EXTRA_HYDRO_LOOP
struct task *t_gradient = NULL;
#endif
@@ -1405,6 +1582,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
struct task *t_limiter = NULL;
struct task *t_star_density = NULL;
struct task *t_star_feedback = NULL;
+ struct task *t_bh_density = NULL;
+ struct task *t_bh_feedback = NULL;
for (int ind = 0; ind < num_elements; ind++) {
@@ -1451,6 +1630,15 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
task_subtype_stars_feedback, flags, 0, ci, NULL);
}
+ /* The black hole feedback tasks */
+ if (with_black_holes) {
+ t_bh_density = scheduler_addtask(
+ sched, task_type_self, task_subtype_bh_density, flags, 0, ci, NULL);
+ t_bh_feedback =
+ scheduler_addtask(sched, task_type_self, task_subtype_bh_feedback,
+ flags, 0, ci, NULL);
+ }
+
/* Link the tasks to the cells */
engine_addlink(e, &ci->hydro.force, t_force);
if (with_limiter) {
@@ -1460,6 +1648,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->stars.density, t_star_density);
engine_addlink(e, &ci->stars.feedback, t_star_feedback);
}
+ if (with_black_holes) {
+ engine_addlink(e, &ci->black_holes.density, t_bh_density);
+ engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
+ }
#ifdef EXTRA_HYDRO_LOOP
@@ -1500,6 +1692,21 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
ci->hydro.super->stars.stars_out);
}
+ if (with_black_holes) {
+
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
+ t_bh_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_bh_density);
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.black_holes_in,
+ t_bh_density);
+ scheduler_addunlock(sched, t_bh_density,
+ ci->hydro.super->black_holes.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.ghost,
+ t_bh_feedback);
+ scheduler_addunlock(sched, t_bh_feedback,
+ ci->hydro.super->black_holes.black_holes_out);
+ }
+
if (with_limiter) {
scheduler_addunlock(sched, ci->super->kick2, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
@@ -1544,6 +1751,14 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
task_subtype_stars_feedback, flags, 0, ci, cj);
}
+ /* The black hole feedback tasks */
+ if (with_black_holes) {
+ t_bh_density = scheduler_addtask(
+ sched, task_type_pair, task_subtype_bh_density, flags, 0, ci, cj);
+ t_bh_feedback = scheduler_addtask(
+ sched, task_type_pair, task_subtype_bh_feedback, flags, 0, ci, cj);
+ }
+
engine_addlink(e, &ci->hydro.force, t_force);
engine_addlink(e, &cj->hydro.force, t_force);
if (with_limiter) {
@@ -1556,6 +1771,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->stars.feedback, t_star_feedback);
engine_addlink(e, &cj->stars.feedback, t_star_feedback);
}
+ if (with_black_holes) {
+ engine_addlink(e, &ci->black_holes.density, t_bh_density);
+ engine_addlink(e, &cj->black_holes.density, t_bh_density);
+ engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
+ engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
+ }
#ifdef EXTRA_HYDRO_LOOP
@@ -1624,6 +1845,22 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
ci->hydro.super->stars.stars_out);
}
+ if (with_black_holes) {
+
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
+ t_bh_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_bh_density);
+ scheduler_addunlock(
+ sched, ci->hydro.super->black_holes.black_holes_in, t_bh_density);
+ scheduler_addunlock(sched, t_bh_density,
+ ci->hydro.super->black_holes.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.ghost,
+ t_bh_feedback);
+ scheduler_addunlock(sched, t_bh_feedback,
+ ci->hydro.super->black_holes.black_holes_out);
+ }
+
if (with_limiter) {
scheduler_addunlock(sched, ci->super->kick2, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
@@ -1660,6 +1897,23 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
cj->hydro.super->stars.stars_out);
}
+ if (with_black_holes) {
+
+ scheduler_addunlock(sched, cj->hydro.super->black_holes.drift,
+ t_bh_density);
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+ t_bh_density);
+ scheduler_addunlock(sched,
+ cj->hydro.super->black_holes.black_holes_in,
+ t_bh_density);
+ scheduler_addunlock(sched, t_bh_density,
+ cj->hydro.super->black_holes.ghost);
+ scheduler_addunlock(sched, cj->hydro.super->black_holes.ghost,
+ t_bh_feedback);
+ scheduler_addunlock(sched, t_bh_feedback,
+ cj->hydro.super->black_holes.black_holes_out);
+ }
+
if (with_limiter) {
scheduler_addunlock(sched, cj->super->kick2, t_limiter);
scheduler_addunlock(sched, t_limiter, cj->super->timestep);
@@ -1702,6 +1956,16 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
task_subtype_stars_feedback, flags, 0, ci, NULL);
}
+ /* The black hole feedback tasks */
+ if (with_black_holes) {
+ t_bh_density =
+ scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_bh_density, flags, 0, ci, NULL);
+ t_bh_feedback =
+ scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_bh_feedback, flags, 0, ci, NULL);
+ }
+
/* Add the link between the new loop and the cell */
engine_addlink(e, &ci->hydro.force, t_force);
if (with_limiter) {
@@ -1711,6 +1975,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->stars.density, t_star_density);
engine_addlink(e, &ci->stars.feedback, t_star_feedback);
}
+ if (with_black_holes) {
+ engine_addlink(e, &ci->black_holes.density, t_bh_density);
+ engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
+ }
#ifdef EXTRA_HYDRO_LOOP
@@ -1757,6 +2025,21 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
ci->hydro.super->stars.stars_out);
}
+ if (with_black_holes) {
+
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
+ t_bh_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_bh_density);
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.black_holes_in,
+ t_bh_density);
+ scheduler_addunlock(sched, t_bh_density,
+ ci->hydro.super->black_holes.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.ghost,
+ t_bh_feedback);
+ scheduler_addunlock(sched, t_bh_feedback,
+ ci->hydro.super->black_holes.black_holes_out);
+ }
+
if (with_limiter) {
scheduler_addunlock(sched, ci->super->kick2, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
@@ -1803,6 +2086,16 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
task_subtype_stars_feedback, flags, 0, ci, cj);
}
+ /* The black hole feedback tasks */
+ if (with_black_holes) {
+ t_bh_density =
+ scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_bh_density, flags, 0, ci, cj);
+ t_bh_feedback =
+ scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_bh_feedback, flags, 0, ci, cj);
+ }
+
engine_addlink(e, &ci->hydro.force, t_force);
engine_addlink(e, &cj->hydro.force, t_force);
if (with_limiter) {
@@ -1815,6 +2108,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->stars.feedback, t_star_feedback);
engine_addlink(e, &cj->stars.feedback, t_star_feedback);
}
+ if (with_black_holes) {
+ engine_addlink(e, &ci->black_holes.density, t_bh_density);
+ engine_addlink(e, &cj->black_holes.density, t_bh_density);
+ engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
+ engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
+ }
#ifdef EXTRA_HYDRO_LOOP
@@ -1882,6 +2181,22 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
ci->hydro.super->stars.stars_out);
}
+ if (with_black_holes) {
+
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
+ t_bh_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_bh_density);
+ scheduler_addunlock(
+ sched, ci->hydro.super->black_holes.black_holes_in, t_bh_density);
+ scheduler_addunlock(sched, t_bh_density,
+ ci->hydro.super->black_holes.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->black_holes.ghost,
+ t_bh_feedback);
+ scheduler_addunlock(sched, t_bh_feedback,
+ ci->hydro.super->black_holes.black_holes_out);
+ }
+
if (with_limiter) {
scheduler_addunlock(sched, ci->super->kick2, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
@@ -1920,6 +2235,23 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
cj->hydro.super->stars.stars_out);
}
+ if (with_black_holes) {
+
+ scheduler_addunlock(sched, cj->hydro.super->black_holes.drift,
+ t_bh_density);
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+ t_bh_density);
+ scheduler_addunlock(sched,
+ cj->hydro.super->black_holes.black_holes_in,
+ t_bh_density);
+ scheduler_addunlock(sched, t_bh_density,
+ cj->hydro.super->black_holes.ghost);
+ scheduler_addunlock(sched, cj->hydro.super->black_holes.ghost,
+ t_bh_feedback);
+ scheduler_addunlock(sched, t_bh_feedback,
+ cj->hydro.super->black_holes.black_holes_out);
+ }
+
if (with_limiter) {
scheduler_addunlock(sched, cj->super->kick2, t_limiter);
scheduler_addunlock(sched, t_limiter, cj->super->timestep);
@@ -1956,6 +2288,7 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
const int periodic = e->s->periodic;
const int with_feedback = (e->policy & engine_policy_feedback);
const int with_stars = (e->policy & engine_policy_stars);
+ const int with_black_holes = (e->policy & engine_policy_black_holes);
struct space *s = e->s;
struct scheduler *sched = &e->sched;
@@ -1978,7 +2311,8 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
struct cell *ci = &cells[cid];
/* Skip cells without hydro or star particles */
- if ((ci->hydro.count == 0) && (!with_stars || ci->stars.count == 0))
+ if ((ci->hydro.count == 0) && (!with_stars || ci->stars.count == 0) &&
+ (!with_black_holes || ci->black_holes.count == 0))
continue;
/* If the cell is local build a self-interaction */
@@ -2008,7 +2342,8 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* Is that neighbour local and does it have gas or star particles ? */
if ((cid >= cjd) ||
((cj->hydro.count == 0) &&
- (!with_feedback || cj->stars.count == 0)) ||
+ (!with_feedback || cj->stars.count == 0) &&
+ (!with_black_holes || cj->black_holes.count == 0)) ||
(ci->nodeID != nodeID && cj->nodeID != nodeID))
continue;
@@ -2083,9 +2418,6 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
struct cell *cj = cell_type_pairs[k].cj;
const int type = cell_type_pairs[k].type;
- /* Add the send task for the particle timesteps. */
- // engine_addtasks_send_timestep(e, ci, cj, NULL, NULL, with_limiter);
-
/* Add the send tasks for the cells in the proxy that have a hydro
* connection. */
if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
@@ -2098,6 +2430,13 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
engine_addtasks_send_stars(e, ci, cj, /*t_feedback=*/NULL, /*t_ti=*/NULL);
+ /* Add the send tasks for the cells in the proxy that have a black holes
+ * connection. */
+ if ((e->policy & engine_policy_black_holes) &&
+ (type & proxy_cell_type_hydro))
+ engine_addtasks_send_black_holes(e, ci, cj, /*t_feedback=*/NULL,
+ /*t_ti=*/NULL);
+
/* Add the send tasks for the cells in the proxy that have a gravity
* connection. */
if ((e->policy & engine_policy_self_gravity) &&
@@ -2117,9 +2456,6 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
struct cell *ci = cell_type_pairs[k].ci;
const int type = cell_type_pairs[k].type;
- /* Add the recv task for the particle timesteps. */
- // engine_addtasks_recv_timestep(e, ci, NULL, NULL, with_limiter);
-
/* Add the recv tasks for the cells in the proxy that have a hydro
* connection. */
if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
@@ -2130,6 +2466,12 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
engine_addtasks_recv_stars(e, ci, NULL, NULL);
+ /* Add the recv tasks for the cells in the proxy that have a black holes
+ * connection. */
+ if ((e->policy & engine_policy_black_holes) &&
+ (type & proxy_cell_type_hydro))
+ engine_addtasks_recv_black_holes(e, ci, NULL, NULL);
+
/* Add the recv tasks for the cells in the proxy that have a gravity
* connection. */
if ((e->policy & engine_policy_self_gravity) &&
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index fd0d99432ef184643245696ed8770db215560386..ba52dd48de43b1388735a440f2a6f842dfbc22c6 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -164,6 +164,37 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_stars(ci, e)) scheduler_activate(s, t);
}
+ /* Activate the black hole density */
+ else if (t_type == task_type_self &&
+ t_subtype == task_subtype_bh_density) {
+ if (cell_is_active_black_holes(ci, e)) {
+ scheduler_activate(s, t);
+ cell_activate_drift_part(ci, s);
+ cell_activate_drift_bpart(ci, s);
+ }
+ }
+
+ /* Store current values of dx_max and h_max. */
+ else if (t_type == task_type_sub_self &&
+ t_subtype == task_subtype_bh_density) {
+ if (cell_is_active_black_holes(ci, e)) {
+ scheduler_activate(s, t);
+ cell_activate_subcell_black_holes_tasks(ci, NULL, s);
+ }
+ }
+
+ else if (t_type == task_type_self &&
+ t_subtype == task_subtype_bh_feedback) {
+ if (cell_is_active_black_holes(ci, e)) {
+ scheduler_activate(s, t);
+ }
+ }
+
+ else if (t_type == task_type_sub_self &&
+ t_subtype == task_subtype_bh_feedback) {
+ if (cell_is_active_black_holes(ci, e)) scheduler_activate(s, t);
+ }
+
/* Activate the gravity drift */
else if (t_type == task_type_self && t_subtype == task_subtype_grav) {
if (cell_is_active_gravity(ci, e)) {
@@ -210,6 +241,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
const int ci_active_stars = cell_is_active_stars(ci, e);
const int cj_active_stars = cell_is_active_stars(cj, e);
+ const int ci_active_black_holes = cell_is_active_black_holes(ci, e);
+ const int cj_active_black_holes = cell_is_active_black_holes(cj, e);
+
/* Only activate tasks that involve a local active cell. */
if ((t_subtype == task_subtype_density ||
t_subtype == task_subtype_gradient ||
@@ -330,6 +364,62 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
}
+ /* Black_Holes density */
+ else if ((t_subtype == task_subtype_bh_density) &&
+ (ci_active_black_holes || cj_active_black_holes) &&
+ (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
+
+ scheduler_activate(s, t);
+
+ /* Set the correct sorting flags */
+ if (t_type == task_type_pair) {
+
+ /* Do ci */
+ if (ci_active_black_holes) {
+
+ /* Activate the drift tasks. */
+ if (ci_nodeID == nodeID) cell_activate_drift_bpart(ci, s);
+ if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
+ }
+
+ /* Do cj */
+ if (cj_active_black_holes) {
+
+ /* Activate the drift tasks. */
+ if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
+ if (cj_nodeID == nodeID) cell_activate_drift_bpart(cj, s);
+ }
+ }
+
+ /* Store current values of dx_max and h_max. */
+ else if (t_type == task_type_sub_pair &&
+ t_subtype == task_subtype_bh_density) {
+ cell_activate_subcell_black_holes_tasks(ci, cj, s);
+ }
+ }
+
+ /* Black_Holes feedback */
+ else if (t_subtype == task_subtype_bh_feedback) {
+
+ /* We only want to activate the task if the cell is active and is
+ going to update some gas on the *local* node */
+ if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
+ (ci_active_black_holes || cj_active_black_holes)) {
+
+ scheduler_activate(s, t);
+
+ } else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
+ (cj_active_black_holes)) {
+
+ scheduler_activate(s, t);
+
+ } else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
+ (ci_active_black_holes)) {
+
+ scheduler_activate(s, t);
+ }
+ }
+
/* Gravity */
else if ((t_subtype == task_subtype_grav) &&
((ci_active_gravity && ci_nodeID == nodeID) ||
@@ -531,6 +621,85 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
#endif
}
+ /* Only interested in black hole density tasks as of here. */
+ else if (t->subtype == task_subtype_bh_density) {
+
+ /* Too much particle movement? */
+ if (cell_need_rebuild_for_black_holes_pair(ci, cj)) *rebuild_space = 1;
+ if (cell_need_rebuild_for_black_holes_pair(cj, ci)) *rebuild_space = 1;
+
+#ifdef WITH_MPI
+ /* Activate the send/recv tasks. */
+ if (ci_nodeID != nodeID) {
+
+ if (cj_active_black_holes) {
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
+
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart,
+ ci_nodeID);
+ cell_activate_drift_bpart(cj, s);
+
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_tend_bpart,
+ ci_nodeID);
+ }
+
+ if (ci_active_black_holes) {
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart);
+
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_tend_bpart);
+
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_xv,
+ ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_rho,
+ ci_nodeID);
+
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(cj, s);
+ }
+
+ } else if (cj_nodeID != nodeID) {
+
+ /* If the local cell is active, receive data from the foreign cell. */
+ if (ci_active_black_holes) {
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
+
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart,
+ cj_nodeID);
+ cell_activate_drift_bpart(ci, s);
+
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_bpart,
+ cj_nodeID);
+ }
+
+ if (cj_active_black_holes) {
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart);
+
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_tend_bpart);
+
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_xv,
+ cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_rho,
+ cj_nodeID);
+
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(ci, s);
+ }
+ }
+#endif
+ }
+
/* Only interested in gravity tasks as of here. */
else if (t_subtype == task_subtype_grav) {
@@ -609,7 +778,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Kick ? */
else if (t_type == task_type_kick1 || t_type == task_type_kick2) {
- if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e))
+ if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
+ cell_is_active_stars(t->ci, e) ||
+ cell_is_active_black_holes(t->ci, e))
scheduler_activate(s, t);
}
@@ -667,12 +838,25 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_stars(t->ci, e)) scheduler_activate(s, t);
}
+ /* Black hole ghost tasks ? */
+ else if (t_type == task_type_bh_ghost) {
+ if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
+ }
+
+ /* Black holes implicit tasks? */
+ else if (t_type == task_type_bh_in || t_type == task_type_bh_out) {
+ if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
+ }
+
/* Time-step? */
else if (t_type == task_type_timestep) {
t->ci->hydro.updated = 0;
t->ci->grav.updated = 0;
t->ci->stars.updated = 0;
- if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e))
+ t->ci->black_holes.updated = 0;
+ if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
+ cell_is_active_stars(t->ci, e) ||
+ cell_is_active_black_holes(t->ci, e))
scheduler_activate(s, t);
}
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index 37d028b64f944cd59cd931cb013b04bcde81c88d..c6bffce2cae2ca6b536cf7046076d38ef3ab544a 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -794,8 +794,6 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
/**
* @brief Initialize the global properties of the feedback scheme.
*
- * By default, takes the values provided by the hydro.
- *
* @param fp The #feedback_props.
* @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system.
diff --git a/src/parser.h b/src/parser.h
index 9bbc1abcebd052a1160d4fdb7bbef7b11fd2bba5..87e618503dca567f7cf381475fcc91f1575b3913 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 256
+#define PARSER_MAX_NO_OF_PARAMS 512
#define PARSER_MAX_NO_OF_SECTIONS 64
/* A parameter in the input file */
diff --git a/src/part.h b/src/part.h
index 629f6dbfa58692fc067546b2eb32c89519dc93d7..a9fdb7889368ca10c2e6f157da1564b679370264 100644
--- a/src/part.h
+++ b/src/part.h
@@ -109,6 +109,8 @@
/* Import the right black hole particle definition */
#if defined(BLACK_HOLES_NONE)
#include "./black_holes/Default/black_holes_part.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes_part.h"
#else
#error "Invalid choice of black hole particle"
#endif
diff --git a/src/runner.c b/src/runner.c
index 3ff1c1274c09936c7b77e939a5c0b6b88750858d..6c782536dbe4202e4919e46aabeac81524d2deb3 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -41,6 +41,8 @@
#include "active.h"
#include "approx_math.h"
#include "atomic.h"
+#include "black_holes.h"
+#include "black_holes_properties.h"
#include "cell.h"
#include "chemistry.h"
#include "const.h"
@@ -126,6 +128,20 @@
#undef FUNCTION_TASK_LOOP
#undef FUNCTION
+/* Import the black hole density loop functions. */
+#define FUNCTION density
+#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
+#include "runner_doiact_black_holes.h"
+#undef FUNCTION_TASK_LOOP
+#undef FUNCTION
+
+/* Import the black hole feedback loop functions. */
+#define FUNCTION feedback
+#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
+#include "runner_doiact_black_holes.h"
+#undef FUNCTION_TASK_LOOP
+#undef FUNCTION
+
/**
* @brief Intermediate task after the density to check that the smoothing
* lengths are correct.
@@ -485,6 +501,352 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_do_stars_ghost);
}
+/**
+ * @brief Intermediate task after the density to check that the smoothing
+ * lengths are correct.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
+
+ struct bpart *restrict bparts = c->black_holes.parts;
+ const struct engine *e = r->e;
+ const int with_cosmology = e->policy & engine_policy_cosmology;
+ const struct cosmology *cosmo = e->cosmology;
+ const float black_holes_h_max = e->hydro_properties->h_max;
+ const float black_holes_h_min = e->hydro_properties->h_min;
+ const float eps = e->black_holes_properties->h_tolerance;
+ const float black_holes_eta_dim =
+ pow_dimension(e->black_holes_properties->eta_neighbours);
+ const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
+ int redo = 0, bcount = 0;
+
+ /* Running value of the maximal smoothing length */
+ double h_max = c->black_holes.h_max;
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if (c->black_holes.count == 0) return;
+ if (!cell_is_active_black_holes(c, e)) return;
+
+ /* Recurse? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ runner_do_black_holes_ghost(r, c->progeny[k], 0);
+
+ /* Update h_max */
+ h_max = max(h_max, c->progeny[k]->black_holes.h_max);
+ }
+ }
+ } else {
+
+ /* Init the list of active particles that have to be updated. */
+ int *sid = NULL;
+ float *h_0 = NULL;
+ float *left = NULL;
+ float *right = NULL;
+ if ((sid = (int *)malloc(sizeof(int) * c->black_holes.count)) == NULL)
+ error("Can't allocate memory for sid.");
+ if ((h_0 = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
+ error("Can't allocate memory for h_0.");
+ if ((left = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
+ error("Can't allocate memory for left.");
+ if ((right = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
+ error("Can't allocate memory for right.");
+ for (int k = 0; k < c->black_holes.count; k++)
+ if (bpart_is_active(&bparts[k], e)) {
+ sid[bcount] = k;
+ h_0[bcount] = bparts[k].h;
+ left[bcount] = 0.f;
+ right[bcount] = black_holes_h_max;
+ ++bcount;
+ }
+
+ /* While there are particles that need to be updated... */
+ for (int num_reruns = 0; bcount > 0 && num_reruns < max_smoothing_iter;
+ num_reruns++) {
+
+ /* Reset the redo-count. */
+ redo = 0;
+
+ /* Loop over the remaining active parts in this cell. */
+ for (int i = 0; i < bcount; i++) {
+
+ /* Get a direct pointer on the part. */
+ struct bpart *bp = &bparts[sid[i]];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Is this part within the timestep? */
+ if (!bpart_is_active(bp, e))
+ error("Ghost applied to inactive particle");
+#endif
+
+ /* Get some useful values */
+ const float h_init = h_0[i];
+ const float h_old = bp->h;
+ const float h_old_dim = pow_dimension(h_old);
+ const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
+
+ float h_new;
+ int has_no_neighbours = 0;
+
+ if (bp->density.wcount == 0.f) { /* No neighbours case */
+
+ /* Flag that there were no neighbours */
+ has_no_neighbours = 1;
+
+ /* Double h and try again */
+ h_new = 2.f * h_old;
+
+ } else {
+
+ /* Finish the density calculation */
+ black_holes_end_density(bp, cosmo);
+
+ /* Compute one step of the Newton-Raphson scheme */
+ const float n_sum = bp->density.wcount * h_old_dim;
+ const float n_target = black_holes_eta_dim;
+ const float f = n_sum - n_target;
+ const float f_prime =
+ bp->density.wcount_dh * h_old_dim +
+ hydro_dimension * bp->density.wcount * h_old_dim_minus_one;
+
+ /* Improve the bisection bounds */
+ if (n_sum < n_target)
+ left[i] = max(left[i], h_old);
+ else if (n_sum > n_target)
+ right[i] = min(right[i], h_old);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check the validity of the left and right bounds */
+ if (left[i] > right[i])
+ error("Invalid left (%e) and right (%e)", left[i], right[i]);
+#endif
+
+ /* Skip if h is already h_max and we don't have enough neighbours */
+ /* Same if we are below h_min */
+ if (((bp->h >= black_holes_h_max) && (f < 0.f)) ||
+ ((bp->h <= black_holes_h_min) && (f > 0.f))) {
+
+ /* Get particle time-step */
+ double dt;
+ if (with_cosmology) {
+ const integertime_t ti_step = get_integer_timestep(bp->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(e->ti_current - 1, bp->time_bin);
+
+ dt = cosmology_get_delta_time(e->cosmology, ti_begin,
+ ti_begin + ti_step);
+ } else {
+ dt = get_timestep(bp->time_bin, e->time_base);
+ }
+
+ /* Compute variables required for the feedback loop */
+ black_holes_prepare_feedback(bp, e->black_holes_properties,
+ e->physical_constants, e->cosmology,
+ dt);
+
+ /* Reset quantities computed by the feedback loop */
+ black_holes_reset_feedback(bp);
+
+ /* Ok, we are done with this particle */
+ continue;
+ }
+
+ /* Normal case: Use Newton-Raphson to get a better value of h */
+
+ /* Avoid floating point exception from f_prime = 0 */
+ h_new = h_old - f / (f_prime + FLT_MIN);
+
+ /* Be verbose about the particles that struggle to converge */
+ if (num_reruns > max_smoothing_iter - 10) {
+
+ message(
+ "Smoothing length convergence problem: iter=%d p->id=%lld "
+ "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
+ "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
+ num_reruns, bp->id, h_init, h_old, h_new, f, f_prime, n_sum,
+ n_target, left[i], right[i]);
+ }
+
+ /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
+ h_new = min(h_new, 2.f * h_old);
+ h_new = max(h_new, 0.5f * h_old);
+
+ /* Verify that we are actually progrssing towards the answer */
+ h_new = max(h_new, left[i]);
+ h_new = min(h_new, right[i]);
+ }
+
+ /* Check whether the particle has an inappropriate smoothing length */
+ if (fabsf(h_new - h_old) > eps * h_old) {
+
+ /* Ok, correct then */
+
+ /* Case where we have been oscillating around the solution */
+ if ((h_new == left[i] && h_old == right[i]) ||
+ (h_old == left[i] && h_new == right[i])) {
+
+ /* Bissect the remaining interval */
+ bp->h = pow_inv_dimension(
+ 0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));
+
+ } else {
+
+ /* Normal case */
+ bp->h = h_new;
+ }
+
+ /* If below the absolute maximum, try again */
+ if (bp->h < black_holes_h_max && bp->h > black_holes_h_min) {
+
+ /* Flag for another round of fun */
+ sid[redo] = sid[i];
+ h_0[redo] = h_0[i];
+ left[redo] = left[i];
+ right[redo] = right[i];
+ redo += 1;
+
+ /* Re-initialise everything */
+ black_holes_init_bpart(bp);
+
+ /* Off we go ! */
+ continue;
+
+ } else if (bp->h <= black_holes_h_min) {
+
+ /* Ok, this particle is a lost cause... */
+ bp->h = black_holes_h_min;
+
+ } else if (bp->h >= black_holes_h_max) {
+
+ /* Ok, this particle is a lost cause... */
+ bp->h = black_holes_h_max;
+
+ /* Do some damage control if no neighbours at all were found */
+ if (has_no_neighbours) {
+ black_holes_bpart_has_no_neighbours(bp, cosmo);
+ }
+
+ } else {
+ error(
+ "Fundamental problem with the smoothing length iteration "
+ "logic.");
+ }
+ }
+
+ /* We now have a particle whose smoothing length has converged */
+
+ /* Check if h_max has increased */
+ h_max = max(h_max, bp->h);
+
+ /* Get particle time-step */
+ double dt;
+ if (with_cosmology) {
+ const integertime_t ti_step = get_integer_timestep(bp->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(e->ti_current - 1, bp->time_bin);
+
+ dt = cosmology_get_delta_time(e->cosmology, ti_begin,
+ ti_begin + ti_step);
+ } else {
+ dt = get_timestep(bp->time_bin, e->time_base);
+ }
+
+ /* Compute variables required for the feedback loop */
+ black_holes_prepare_feedback(bp, e->black_holes_properties,
+ e->physical_constants, e->cosmology, dt);
+
+ /* Reset quantities computed by the feedback loop */
+ black_holes_reset_feedback(bp);
+ }
+
+ /* We now need to treat the particles whose smoothing length had not
+ * converged again */
+
+ /* Re-set the counter for the next loop (potentially). */
+ bcount = redo;
+ if (bcount > 0) {
+
+ /* Climb up the cell hierarchy. */
+ for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
+
+ /* Run through this cell's density interactions. */
+ for (struct link *l = finger->black_holes.density; l != NULL;
+ l = l->next) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (l->t->ti_run < r->e->ti_current)
+ error("Density task should have been run.");
+#endif
+
+ /* Self-interaction? */
+ if (l->t->type == task_type_self)
+ runner_doself_subset_branch_bh_density(r, finger, bparts, sid,
+ bcount);
+
+ /* Otherwise, pair interaction? */
+ else if (l->t->type == task_type_pair) {
+
+ /* Left or right? */
+ if (l->t->ci == finger)
+ runner_dopair_subset_branch_bh_density(r, finger, bparts, sid,
+ bcount, l->t->cj);
+ else
+ runner_dopair_subset_branch_bh_density(r, finger, bparts, sid,
+ bcount, l->t->ci);
+ }
+
+ /* Otherwise, sub-self interaction? */
+ else if (l->t->type == task_type_sub_self)
+ runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
+ NULL, 1);
+
+ /* Otherwise, sub-pair interaction? */
+ else if (l->t->type == task_type_sub_pair) {
+
+ /* Left or right? */
+ if (l->t->ci == finger)
+ runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
+ l->t->cj, 1);
+ else
+ runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
+ l->t->ci, 1);
+ }
+ }
+ }
+ }
+ }
+
+ if (bcount) {
+ error("Smoothing length failed to converge on %i particles.", bcount);
+ }
+
+ /* Be clean */
+ free(left);
+ free(right);
+ free(sid);
+ free(h_0);
+ }
+
+ /* Update h_max */
+ c->black_holes.h_max = h_max;
+
+ /* The ghost may not always be at the top level.
+ * Therefore we need to update h_max between the super- and top-levels */
+ if (c->black_holes.ghost) {
+ for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
+ atomic_max_d(&tmp->black_holes.h_max, h_max);
+ }
+ }
+
+ if (timer) TIMER_TOC(timer_do_black_holes_ghost);
+}
+
/**
* @brief Calculate gravity acceleration from external potential
*
@@ -1999,6 +2361,35 @@ static void runner_do_unskip_stars(struct cell *c, struct engine *e) {
if (forcerebuild) atomic_inc(&e->forcerebuild);
}
+/**
+ * @brief Unskip any black hole tasks associated with active cells.
+ *
+ * @param c The cell.
+ * @param e The engine.
+ */
+static void runner_do_unskip_black_holes(struct cell *c, struct engine *e) {
+
+ /* Ignore empty cells. */
+ if (c->black_holes.count == 0) return;
+
+ /* Skip inactive cells. */
+ if (!cell_is_active_black_holes(c, e)) return;
+
+ /* Recurse */
+ if (c->split) {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+ runner_do_unskip_black_holes(cp, e);
+ }
+ }
+ }
+
+ /* Unskip any active tasks. */
+ const int forcerebuild = cell_unskip_black_holes_tasks(c, &e->sched);
+ if (forcerebuild) atomic_inc(&e->forcerebuild);
+}
+
/**
* @brief Unskip any gravity tasks associated with active cells.
*
@@ -2056,6 +2447,10 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
/* Stars tasks */
if (e->policy & engine_policy_stars) runner_do_unskip_stars(c, e);
+
+ /* Black hole tasks */
+ if (e->policy & engine_policy_black_holes)
+ runner_do_unskip_black_holes(c, e);
}
}
}
@@ -3450,6 +3845,92 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
#endif
}
+/**
+ * @brief Construct the cell properties from the received #bpart.
+ *
+ * Note that we do not need to clear the sorts since we do not sort
+ * the black holes.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param clear_sorts Should we clear the sort flag and hence trigger a sort ?
+ * @param timer Are we timing this ?
+ */
+void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
+ int timer) {
+
+#ifdef WITH_MPI
+
+ struct bpart *restrict bparts = c->black_holes.parts;
+ const size_t nr_bparts = c->black_holes.count;
+ const integertime_t ti_current = r->e->ti_current;
+
+ TIMER_TIC;
+
+ integertime_t ti_black_holes_end_min = max_nr_timesteps;
+ integertime_t ti_black_holes_end_max = 0;
+ timebin_t time_bin_min = num_time_bins;
+ timebin_t time_bin_max = 0;
+ float h_max = 0.f;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID == engine_rank) error("Updating a local cell!");
+#endif
+
+ /* If this cell is a leaf, collect the particle data. */
+ if (!c->split) {
+
+ /* Collect everything... */
+ for (size_t k = 0; k < nr_bparts; k++) {
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+ bparts[k].num_ngb_force = 0;
+#endif
+ if (bparts[k].time_bin == time_bin_inhibited) continue;
+ time_bin_min = min(time_bin_min, bparts[k].time_bin);
+ time_bin_max = max(time_bin_max, bparts[k].time_bin);
+ h_max = max(h_max, bparts[k].h);
+ }
+
+ /* Convert into a time */
+ ti_black_holes_end_min = get_integer_time_end(ti_current, time_bin_min);
+ ti_black_holes_end_max = get_integer_time_end(ti_current, time_bin_max);
+ }
+
+ /* Otherwise, recurse and collect. */
+ else {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL && c->progeny[k]->black_holes.count > 0) {
+ runner_do_recv_bpart(r, c->progeny[k], clear_sorts, 0);
+ ti_black_holes_end_min =
+ min(ti_black_holes_end_min, c->progeny[k]->black_holes.ti_end_min);
+ ti_black_holes_end_max =
+ max(ti_black_holes_end_max, c->progeny[k]->grav.ti_end_max);
+ h_max = max(h_max, c->progeny[k]->black_holes.h_max);
+ }
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ti_black_holes_end_min < ti_current)
+ error(
+ "Received a cell at an incorrect time c->ti_end_min=%lld, "
+ "e->ti_current=%lld.",
+ ti_black_holes_end_min, ti_current);
+#endif
+
+ /* ... and store. */
+ // c->grav.ti_end_min = ti_gravity_end_min;
+ // c->grav.ti_end_max = ti_gravity_end_max;
+ c->black_holes.ti_old_part = ti_current;
+ c->black_holes.h_max = h_max;
+
+ if (timer) TIMER_TOC(timer_dorecv_bpart);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
/**
* @brief The #runner main thread routine.
*
@@ -3532,6 +4013,10 @@ void *runner_main(void *data) {
runner_doself_branch_stars_density(r, ci);
else if (t->subtype == task_subtype_stars_feedback)
runner_doself_branch_stars_feedback(r, ci);
+ else if (t->subtype == task_subtype_bh_density)
+ runner_doself_branch_bh_density(r, ci);
+ else if (t->subtype == task_subtype_bh_feedback)
+ runner_doself_branch_bh_feedback(r, ci);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -3553,6 +4038,10 @@ void *runner_main(void *data) {
runner_dopair_branch_stars_density(r, ci, cj);
else if (t->subtype == task_subtype_stars_feedback)
runner_dopair_branch_stars_feedback(r, ci, cj);
+ else if (t->subtype == task_subtype_bh_density)
+ runner_dopair_branch_bh_density(r, ci, cj);
+ else if (t->subtype == task_subtype_bh_feedback)
+ runner_dopair_branch_bh_feedback(r, ci, cj);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -3572,6 +4061,10 @@ void *runner_main(void *data) {
runner_dosub_self_stars_density(r, ci, 1);
else if (t->subtype == task_subtype_stars_feedback)
runner_dosub_self_stars_feedback(r, ci, 1);
+ else if (t->subtype == task_subtype_bh_density)
+ runner_dosub_self_bh_density(r, ci, 1);
+ else if (t->subtype == task_subtype_bh_feedback)
+ runner_dosub_self_bh_feedback(r, ci, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -3591,6 +4084,10 @@ void *runner_main(void *data) {
runner_dosub_pair_stars_density(r, ci, cj, 1);
else if (t->subtype == task_subtype_stars_feedback)
runner_dosub_pair_stars_feedback(r, ci, cj, 1);
+ else if (t->subtype == task_subtype_bh_density)
+ runner_dosub_pair_bh_density(r, ci, cj, 1);
+ else if (t->subtype == task_subtype_bh_feedback)
+ runner_dosub_pair_bh_feedback(r, ci, cj, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -3625,6 +4122,9 @@ void *runner_main(void *data) {
case task_type_stars_ghost:
runner_do_stars_ghost(r, ci, 1);
break;
+ case task_type_bh_ghost:
+ runner_do_black_holes_ghost(r, ci, 1);
+ break;
case task_type_drift_part:
runner_do_drift_part(r, ci, 1);
break;
@@ -3699,6 +4199,8 @@ void *runner_main(void *data) {
runner_do_recv_gpart(r, ci, 1);
} else if (t->subtype == task_subtype_spart) {
runner_do_recv_spart(r, ci, 1, 1);
+ } else if (t->subtype == task_subtype_bpart) {
+ runner_do_recv_bpart(r, ci, 1, 1);
} else if (t->subtype == task_subtype_multipole) {
cell_unpack_multipoles(ci, (struct gravity_tensors *)t->buff);
free(t->buff);
diff --git a/src/runner_doiact_black_holes.h b/src/runner_doiact_black_holes.h
new file mode 100644
index 0000000000000000000000000000000000000000..efd2d02b3aca3c530b4e1af13e38328368c99f7c
--- /dev/null
+++ b/src/runner_doiact_black_holes.h
@@ -0,0 +1,765 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Before including this file, define FUNCTION, which is the
+ name of the interaction function. This creates the interaction functions
+ runner_dopair_FUNCTION, runner_doself_FUNCTION and runner_dosub_FUNCTION
+ calling the pairwise interaction function runner_iact_FUNCTION. */
+
+#define PASTE(x, y) x##_##y
+
+#define _DOSELF1_BH(f) PASTE(runner_doself_bh, f)
+#define DOSELF1_BH _DOSELF1_BH(FUNCTION)
+
+#define _DO_SYM_PAIR1_BH(f) PASTE(runner_do_sym_pair_bh, f)
+#define DO_SYM_PAIR1_BH _DO_SYM_PAIR1_BH(FUNCTION)
+
+#define _DO_NONSYM_PAIR1_BH_NAIVE(f) PASTE(runner_do_nonsym_pair_bh_naive, f)
+#define DO_NONSYM_PAIR1_BH_NAIVE _DO_NONSYM_PAIR1_BH_NAIVE(FUNCTION)
+
+#define _DOPAIR1_BH_NAIVE(f) PASTE(runner_dopair_bh_naive, f)
+#define DOPAIR1_BH_NAIVE _DOPAIR1_BH_NAIVE(FUNCTION)
+
+#define _DOPAIR1_SUBSET_BH(f) PASTE(runner_dopair_subset_bh, f)
+#define DOPAIR1_SUBSET_BH _DOPAIR1_SUBSET_BH(FUNCTION)
+
+#define _DOPAIR1_SUBSET_BH_NAIVE(f) PASTE(runner_dopair_subset_bh_naive, f)
+#define DOPAIR1_SUBSET_BH_NAIVE _DOPAIR1_SUBSET_BH_NAIVE(FUNCTION)
+
+#define _DOSELF1_SUBSET_BH(f) PASTE(runner_doself_subset_bh, f)
+#define DOSELF1_SUBSET_BH _DOSELF1_SUBSET_BH(FUNCTION)
+
+#define _DOSELF1_SUBSET_BRANCH_BH(f) PASTE(runner_doself_subset_branch_bh, f)
+#define DOSELF1_SUBSET_BRANCH_BH _DOSELF1_SUBSET_BRANCH_BH(FUNCTION)
+
+#define _DOPAIR1_SUBSET_BRANCH_BH(f) PASTE(runner_dopair_subset_branch_bh, f)
+#define DOPAIR1_SUBSET_BRANCH_BH _DOPAIR1_SUBSET_BRANCH_BH(FUNCTION)
+
+#define _DOSUB_SUBSET_BH(f) PASTE(runner_dosub_subset_bh, f)
+#define DOSUB_SUBSET_BH _DOSUB_SUBSET_BH(FUNCTION)
+
+#define _DOSELF1_BRANCH_BH(f) PASTE(runner_doself_branch_bh, f)
+#define DOSELF1_BRANCH_BH _DOSELF1_BRANCH_BH(FUNCTION)
+
+#define _DOPAIR1_BRANCH_BH(f) PASTE(runner_dopair_branch_bh, f)
+#define DOPAIR1_BRANCH_BH _DOPAIR1_BRANCH_BH(FUNCTION)
+
+#define _DOSUB_PAIR1_BH(f) PASTE(runner_dosub_pair_bh, f)
+#define DOSUB_PAIR1_BH _DOSUB_PAIR1_BH(FUNCTION)
+
+#define _DOSUB_SELF1_BH(f) PASTE(runner_dosub_self_bh, f)
+#define DOSUB_SELF1_BH _DOSUB_SELF1_BH(FUNCTION)
+
+#define _TIMER_DOSELF_BH(f) PASTE(timer_doself_bh, f)
+#define TIMER_DOSELF_BH _TIMER_DOSELF_BH(FUNCTION)
+
+#define _TIMER_DOPAIR_BH(f) PASTE(timer_dopair_bh, f)
+#define TIMER_DOPAIR_BH _TIMER_DOPAIR_BH(FUNCTION)
+
+#define _TIMER_DOSUB_SELF_BH(f) PASTE(timer_dosub_self_bh, f)
+#define TIMER_DOSUB_SELF_BH _TIMER_DOSUB_SELF_BH(FUNCTION)
+
+#define _TIMER_DOSUB_PAIR_BH(f) PASTE(timer_dosub_pair_bh, f)
+#define TIMER_DOSUB_PAIR_BH _TIMER_DOSUB_PAIR_BH(FUNCTION)
+
+#define _IACT_BH(f) PASTE(runner_iact_nonsym_bh, f)
+#define IACT_BH _IACT_BH(FUNCTION)
+
+/**
+ * @brief Calculate the number density of #part around the #bpart
+ *
+ * @param r runner task
+ * @param c cell
+ * @param timer 1 if the time is to be recorded.
+ */
+void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID != engine_rank) error("Should be run on a different node");
+#endif
+
+ TIMER_TIC;
+
+ const struct engine *e = r->e;
+ const integertime_t ti_current = e->ti_current;
+ const struct cosmology *cosmo = e->cosmology;
+
+ /* Anything to do here? */
+ if (c->hydro.count == 0 || c->black_holes.count == 0) return;
+ if (!cell_is_active_black_holes(c, e)) return;
+
+ const int bcount = c->black_holes.count;
+ const int count = c->hydro.count;
+ struct bpart *restrict bparts = c->black_holes.parts;
+ struct part *restrict parts = c->hydro.parts;
+ struct xpart *restrict xparts = c->hydro.xparts;
+
+ /* Loop over the bparts in ci. */
+ for (int sid = 0; sid < bcount; sid++) {
+
+ /* Get a hold of the ith bpart in ci. */
+ struct bpart *restrict si = &bparts[sid];
+
+ /* Skip inactive particles */
+ if (!bpart_is_active(si, e)) continue;
+
+ const float hi = si->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+ const float six[3] = {(float)(si->x[0] - c->loc[0]),
+ (float)(si->x[1] - c->loc[1]),
+ (float)(si->x[2] - c->loc[2])};
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts[pjd];
+ struct xpart *restrict xpj = &xparts[pjd];
+ const float hj = pj->h;
+
+ /* Early abort? */
+ if (part_is_inhibited(pj, e)) continue;
+
+ /* Compute the pairwise distance. */
+ const float pjx[3] = {(float)(pj->x[0] - c->loc[0]),
+ (float)(pj->x[1] - c->loc[1]),
+ (float)(pj->x[2] - c->loc[2])};
+ float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
+ if (r2 < hig2) {
+ IACT_BH(r2, dx, hi, hj, si, pj, xpj, cosmo, ti_current);
+ }
+ } /* loop over the parts in ci. */
+ } /* loop over the bparts in ci. */
+
+ TIMER_TOC(TIMER_DOSELF_BH);
+}
+
+/**
+ * @brief Calculate the number density of cj #part around the ci #bpart
+ *
+ * @param r runner task
+ * @param ci The first #cell
+ * @param cj The second #cell
+ */
+void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
+ struct cell *restrict cj) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ if (ci->nodeID != engine_rank) error("Should be run on a different node");
+#else
+ if (cj->nodeID != engine_rank) error("Should be run on a different node");
+#endif
+#endif
+
+ const struct engine *e = r->e;
+ const integertime_t ti_current = e->ti_current;
+ const struct cosmology *cosmo = e->cosmology;
+
+ /* Anything to do here? */
+ if (cj->hydro.count == 0 || ci->black_holes.count == 0) return;
+ if (!cell_is_active_black_holes(ci, e)) return;
+
+ const int bcount_i = ci->black_holes.count;
+ const int count_j = cj->hydro.count;
+ struct bpart *restrict bparts_i = ci->black_holes.parts;
+ struct part *restrict parts_j = cj->hydro.parts;
+ struct xpart *restrict xparts_j = cj->hydro.xparts;
+
+ /* Get the relative distance between the pairs, wrapping. */
+ double shift[3] = {0.0, 0.0, 0.0};
+ for (int k = 0; k < 3; k++) {
+ if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2)
+ shift[k] = e->s->dim[k];
+ else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2)
+ shift[k] = -e->s->dim[k];
+ }
+
+ /* Loop over the bparts in ci. */
+ for (int sid = 0; sid < bcount_i; sid++) {
+
+ /* Get a hold of the ith bpart in ci. */
+ struct bpart *restrict si = &bparts_i[sid];
+
+ /* Skip inactive particles */
+ if (!bpart_is_active(si, e)) continue;
+
+ const float hi = si->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+ const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])),
+ (float)(si->x[1] - (cj->loc[1] + shift[1])),
+ (float)(si->x[2] - (cj->loc[2] + shift[2]))};
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[pjd];
+ struct xpart *restrict xpj = &xparts_j[pjd];
+ const float hj = pj->h;
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
+
+ /* Compute the pairwise distance. */
+ const float pjx[3] = {(float)(pj->x[0] - cj->loc[0]),
+ (float)(pj->x[1] - cj->loc[1]),
+ (float)(pj->x[2] - cj->loc[2])};
+ float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
+ if (r2 < hig2) {
+ IACT_BH(r2, dx, hi, hj, si, pj, xpj, cosmo, ti_current);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+}
+
+void DOPAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
+ struct cell *restrict cj, int timer) {
+
+ TIMER_TIC;
+
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ const int do_ci_bh = ci->nodeID == r->e->nodeID;
+ const int do_cj_bh = cj->nodeID == r->e->nodeID;
+#else
+ /* here we are updating the hydro -> switch ci, cj */
+ const int do_ci_bh = cj->nodeID == r->e->nodeID;
+ const int do_cj_bh = ci->nodeID == r->e->nodeID;
+#endif
+ if (do_ci_bh && ci->black_holes.count != 0 && cj->hydro.count != 0)
+ DO_NONSYM_PAIR1_BH_NAIVE(r, ci, cj);
+ if (do_cj_bh && cj->black_holes.count != 0 && ci->hydro.count != 0)
+ DO_NONSYM_PAIR1_BH_NAIVE(r, cj, ci);
+
+ TIMER_TOC(TIMER_DOPAIR_BH);
+}
+
+/**
+ * @brief Compute the interactions between a cell pair, but only for the
+ * given indices in ci.
+ *
+ * Version using a brute-force algorithm.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param bparts_i The #bpart to interact with @c cj.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param bcount The number of particles in @c ind.
+ * @param cj The second #cell.
+ * @param shift The shift vector to apply to the particles in ci.
+ */
+void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
+ struct bpart *restrict bparts_i, int *restrict ind,
+ const int bcount, struct cell *restrict cj,
+ const double *shift) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->nodeID != engine_rank) error("Should be run on a different node");
+#endif
+
+ const struct engine *e = r->e;
+ const integertime_t ti_current = e->ti_current;
+ const struct cosmology *cosmo = e->cosmology;
+
+ const int count_j = cj->hydro.count;
+ struct part *restrict parts_j = cj->hydro.parts;
+ struct xpart *restrict xparts_j = cj->hydro.xparts;
+
+ /* Early abort? */
+ if (count_j == 0) return;
+
+ /* Loop over the parts_i. */
+ for (int pid = 0; pid < bcount; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct bpart *restrict bpi = &bparts_i[ind[pid]];
+
+ const double pix = bpi->x[0] - (shift[0]);
+ const double piy = bpi->x[1] - (shift[1]);
+ const double piz = bpi->x[2] - (shift[2]);
+ const float hi = bpi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (!bpart_is_active(bpi, e))
+ error("Trying to correct smoothing length of inactive particle !");
+#endif
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[pjd];
+ struct xpart *restrict xpj = &xparts_j[pjd];
+
+ /* Skip inhibited particles */
+ if (part_is_inhibited(pj, e)) continue;
+
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+ const float hj = pj->h;
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy),
+ (float)(piz - pjz)};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+ /* Hit or miss? */
+ if (r2 < hig2) {
+ IACT_BH(r2, dx, hi, hj, bpi, pj, xpj, cosmo, ti_current);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+}
+
+/**
+ * @brief Compute the interactions between a cell pair, but only for the
+ * given indices in ci.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param bparts The #bpart to interact.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param bcount The number of particles in @c ind.
+ */
+void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
+ struct bpart *restrict bparts, int *restrict ind,
+ const int bcount) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->nodeID != engine_rank) error("Should be run on a different node");
+#endif
+
+ const struct engine *e = r->e;
+ const integertime_t ti_current = e->ti_current;
+ const struct cosmology *cosmo = e->cosmology;
+
+ const int count_i = ci->hydro.count;
+ struct part *restrict parts_j = ci->hydro.parts;
+ struct xpart *restrict xparts_j = ci->hydro.xparts;
+
+ /* Early abort? */
+ if (count_i == 0) return;
+
+ /* Loop over the parts in ci. */
+ for (int bpid = 0; bpid < bcount; bpid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct bpart *bpi = &bparts[ind[bpid]];
+ const float bpix[3] = {(float)(bpi->x[0] - ci->loc[0]),
+ (float)(bpi->x[1] - ci->loc[1]),
+ (float)(bpi->x[2] - ci->loc[2])};
+ const float hi = bpi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (!bpart_is_active(bpi, e))
+ error("Inactive particle in subset function!");
+#endif
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_i; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[pjd];
+ struct xpart *restrict xpj = &xparts_j[pjd];
+
+ /* Early abort? */
+ if (part_is_inhibited(pj, e)) continue;
+
+ /* Compute the pairwise distance. */
+ const float pjx[3] = {(float)(pj->x[0] - ci->loc[0]),
+ (float)(pj->x[1] - ci->loc[1]),
+ (float)(pj->x[2] - ci->loc[2])};
+ float dx[3] = {bpix[0] - pjx[0], bpix[1] - pjx[1], bpix[2] - pjx[2]};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
+ /* Hit or miss? */
+ if (r2 < hig2) {
+ IACT_BH(r2, dx, hi, pj->h, bpi, pj, xpj, cosmo, ti_current);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+}
+
+/**
+ * @brief Determine which version of DOSELF1_SUBSET_BH needs to be called
+ * depending on the optimisation level.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param bparts The #bpart to interact.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param bcount The number of particles in @c ind.
+ */
+void DOSELF1_SUBSET_BRANCH_BH(struct runner *r, struct cell *restrict ci,
+ struct bpart *restrict bparts, int *restrict ind,
+ const int bcount) {
+
+ DOSELF1_SUBSET_BH(r, ci, bparts, ind, bcount);
+}
+
+/**
+ * @brief Determine which version of DOPAIR1_SUBSET_BH needs to be called
+ * depending on the orientation of the cells or whether DOPAIR1_SUBSET_BH
+ * needs to be called at all.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param bparts_i The #bpart to interact with @c cj.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param bcount The number of particles in @c ind.
+ * @param cj The second #cell.
+ */
+void DOPAIR1_SUBSET_BRANCH_BH(struct runner *r, struct cell *restrict ci,
+ struct bpart *restrict bparts_i,
+ int *restrict ind, int const bcount,
+ struct cell *restrict cj) {
+
+ const struct engine *e = r->e;
+
+ /* Anything to do here? */
+ if (cj->hydro.count == 0) return;
+
+ /* Get the relative distance between the pairs, wrapping. */
+ double shift[3] = {0.0, 0.0, 0.0};
+ for (int k = 0; k < 3; k++) {
+ if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2)
+ shift[k] = e->s->dim[k];
+ else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2)
+ shift[k] = -e->s->dim[k];
+ }
+
+ DOPAIR1_SUBSET_BH_NAIVE(r, ci, bparts_i, ind, bcount, cj, shift);
+}
+
+void DOSUB_SUBSET_BH(struct runner *r, struct cell *ci, struct bpart *bparts,
+ int *ind, const int bcount, struct cell *cj,
+ int gettimer) {
+
+ const struct engine *e = r->e;
+ struct space *s = e->s;
+
+ /* Should we even bother? */
+ if (!cell_is_active_black_holes(ci, e) &&
+ (cj == NULL || !cell_is_active_black_holes(cj, e)))
+ return;
+
+ /* Find out in which sub-cell of ci the parts are. */
+ struct cell *sub = NULL;
+ if (ci->split) {
+ for (int k = 0; k < 8; k++) {
+ if (ci->progeny[k] != NULL) {
+ if (&bparts[ind[0]] >= &ci->progeny[k]->black_holes.parts[0] &&
+ &bparts[ind[0]] <
+ &ci->progeny[k]
+ ->black_holes.parts[ci->progeny[k]->black_holes.count]) {
+ sub = ci->progeny[k];
+ break;
+ }
+ }
+ }
+ }
+
+ /* Is this a single cell? */
+ if (cj == NULL) {
+
+ /* Recurse? */
+ if (cell_can_recurse_in_self_black_holes_task(ci)) {
+
+ /* Loop over all progeny. */
+ DOSUB_SUBSET_BH(r, sub, bparts, ind, bcount, NULL, 0);
+ for (int j = 0; j < 8; j++)
+ if (ci->progeny[j] != sub && ci->progeny[j] != NULL)
+ DOSUB_SUBSET_BH(r, sub, bparts, ind, bcount, ci->progeny[j], 0);
+
+ }
+
+ /* Otherwise, compute self-interaction. */
+ else
+ DOSELF1_SUBSET_BRANCH_BH(r, ci, bparts, ind, bcount);
+ } /* self-interaction. */
+
+ /* Otherwise, it's a pair interaction. */
+ else {
+
+ /* Recurse? */
+ if (cell_can_recurse_in_pair_black_holes_task(ci, cj) &&
+ cell_can_recurse_in_pair_black_holes_task(cj, ci)) {
+
+ /* Get the type of pair and flip ci/cj if needed. */
+ double shift[3] = {0.0, 0.0, 0.0};
+ const int sid = space_getsid(s, &ci, &cj, shift);
+
+ struct cell_split_pair *csp = &cell_split_pairs[sid];
+ for (int k = 0; k < csp->count; k++) {
+ const int pid = csp->pairs[k].pid;
+ const int pjd = csp->pairs[k].pjd;
+ if (ci->progeny[pid] == sub && cj->progeny[pjd] != NULL)
+ DOSUB_SUBSET_BH(r, ci->progeny[pid], bparts, ind, bcount,
+ cj->progeny[pjd], 0);
+ if (ci->progeny[pid] != NULL && cj->progeny[pjd] == sub)
+ DOSUB_SUBSET_BH(r, cj->progeny[pjd], bparts, ind, bcount,
+ ci->progeny[pid], 0);
+ }
+ }
+
+ /* Otherwise, compute the pair directly. */
+ else if (cell_is_active_black_holes(ci, e) && cj->hydro.count > 0) {
+
+ /* Do any of the cells need to be drifted first? */
+ if (cell_is_active_black_holes(ci, e)) {
+ if (!cell_are_bpart_drifted(ci, e)) error("Cell should be drifted!");
+ if (!cell_are_part_drifted(cj, e)) error("Cell should be drifted!");
+ }
+
+ DOPAIR1_SUBSET_BRANCH_BH(r, ci, bparts, ind, bcount, cj);
+ }
+
+ } /* otherwise, pair interaction. */
+}
+
+/**
+ * @brief Determine which version of DOSELF1_BH needs to be called depending
+ * on the optimisation level.
+ *
+ * @param r #runner
+ * @param c #cell c
+ *
+ */
+void DOSELF1_BRANCH_BH(struct runner *r, struct cell *c) {
+
+ const struct engine *restrict e = r->e;
+
+ /* Anything to do here? */
+ if (c->black_holes.count == 0) return;
+
+ /* Anything to do here? */
+ if (!cell_is_active_black_holes(c, e)) return;
+
+ /* Did we mess up the recursion? */
+ if (c->black_holes.h_max_old * kernel_gamma > c->dmin)
+ error("Cell smaller than smoothing length");
+
+ DOSELF1_BH(r, c, 1);
+}
+
+/**
+ * @brief Determine which version of DOPAIR1_BH needs to be called depending
+ * on the orientation of the cells or whether DOPAIR1_BH needs to be called
+ * at all.
+ *
+ * @param r #runner
+ * @param ci #cell ci
+ * @param cj #cell cj
+ *
+ */
+void DOPAIR1_BRANCH_BH(struct runner *r, struct cell *ci, struct cell *cj) {
+
+ const struct engine *restrict e = r->e;
+
+ const int ci_active = cell_is_active_black_holes(ci, e);
+ const int cj_active = cell_is_active_black_holes(cj, e);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ const int do_ci_bh = ci->nodeID == e->nodeID;
+ const int do_cj_bh = cj->nodeID == e->nodeID;
+#else
+ /* here we are updating the hydro -> switch ci, cj */
+ const int do_ci_bh = cj->nodeID == e->nodeID;
+ const int do_cj_bh = ci->nodeID == e->nodeID;
+#endif
+ const int do_ci = (ci->black_holes.count != 0 && cj->hydro.count != 0 &&
+ ci_active && do_ci_bh);
+ const int do_cj = (cj->black_holes.count != 0 && ci->hydro.count != 0 &&
+ cj_active && do_cj_bh);
+
+ /* Anything to do here? */
+ if (!do_ci && !do_cj) return;
+
+ /* Check that cells are drifted. */
+ if (do_ci &&
+ (!cell_are_bpart_drifted(ci, e) || !cell_are_part_drifted(cj, e)))
+ error("Interacting undrifted cells.");
+
+ if (do_cj &&
+ (!cell_are_part_drifted(ci, e) || !cell_are_bpart_drifted(cj, e)))
+ error("Interacting undrifted cells.");
+
+ /* No sorted intreactions here -> use the naive ones */
+ DOPAIR1_BH_NAIVE(r, ci, cj, 1);
+}
+
+/**
+ * @brief Compute grouped sub-cell interactions for pairs
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ * @param gettimer Do we have a timer ?
+ *
+ * @todo Hard-code the sid on the recursive calls to avoid the
+ * redundant computations to find the sid on-the-fly.
+ */
+void DOSUB_PAIR1_BH(struct runner *r, struct cell *ci, struct cell *cj,
+ int gettimer) {
+
+ TIMER_TIC;
+
+ struct space *s = r->e->s;
+ const struct engine *e = r->e;
+
+ /* Should we even bother? */
+ const int should_do_ci = ci->black_holes.count != 0 && cj->hydro.count != 0 &&
+ cell_is_active_black_holes(ci, e);
+ const int should_do_cj = cj->black_holes.count != 0 && ci->hydro.count != 0 &&
+ cell_is_active_black_holes(cj, e);
+ if (!should_do_ci && !should_do_cj) return;
+
+ /* Get the type of pair and flip ci/cj if needed. */
+ double shift[3];
+ const int sid = space_getsid(s, &ci, &cj, shift);
+
+ /* Recurse? */
+ if (cell_can_recurse_in_pair_black_holes_task(ci, cj) &&
+ cell_can_recurse_in_pair_black_holes_task(cj, ci)) {
+ struct cell_split_pair *csp = &cell_split_pairs[sid];
+ for (int k = 0; k < csp->count; k++) {
+ const int pid = csp->pairs[k].pid;
+ const int pjd = csp->pairs[k].pjd;
+ if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL)
+ DOSUB_PAIR1_BH(r, ci->progeny[pid], cj->progeny[pjd], 0);
+ }
+ }
+
+ /* Otherwise, compute the pair directly. */
+ else {
+
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ const int do_ci_bh = ci->nodeID == e->nodeID;
+ const int do_cj_bh = cj->nodeID == e->nodeID;
+#else
+ /* here we are updating the hydro -> switch ci, cj */
+ const int do_ci_bh = cj->nodeID == e->nodeID;
+ const int do_cj_bh = ci->nodeID == e->nodeID;
+#endif
+ const int do_ci = ci->black_holes.count != 0 && cj->hydro.count != 0 &&
+ cell_is_active_black_holes(ci, e) && do_ci_bh;
+ const int do_cj = cj->black_holes.count != 0 && ci->hydro.count != 0 &&
+ cell_is_active_black_holes(cj, e) && do_cj_bh;
+
+ if (do_ci) {
+
+ /* Make sure both cells are drifted to the current timestep. */
+ if (!cell_are_bpart_drifted(ci, e))
+ error("Interacting undrifted cells (bparts).");
+
+ if (!cell_are_part_drifted(cj, e))
+ error("Interacting undrifted cells (parts).");
+ }
+
+ if (do_cj) {
+
+ /* Make sure both cells are drifted to the current timestep. */
+ if (!cell_are_part_drifted(ci, e))
+ error("Interacting undrifted cells (parts).");
+
+ if (!cell_are_bpart_drifted(cj, e))
+ error("Interacting undrifted cells (bparts).");
+ }
+
+ if (do_ci || do_cj) DOPAIR1_BRANCH_BH(r, ci, cj);
+ }
+
+ TIMER_TOC(TIMER_DOSUB_PAIR_BH);
+}
+
+/**
+ * @brief Compute grouped sub-cell interactions for self tasks
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param gettimer Do we have a timer ?
+ */
+void DOSUB_SELF1_BH(struct runner *r, struct cell *ci, int gettimer) {
+
+ TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->nodeID != engine_rank)
+ error("This function should not be called on foreign cells");
+#endif
+
+ /* Should we even bother? */
+ if (ci->hydro.count == 0 || ci->black_holes.count == 0 ||
+ !cell_is_active_black_holes(ci, r->e))
+ return;
+
+ /* Recurse? */
+ if (cell_can_recurse_in_self_black_holes_task(ci)) {
+
+ /* Loop over all progeny. */
+ for (int k = 0; k < 8; k++)
+ if (ci->progeny[k] != NULL) {
+ DOSUB_SELF1_BH(r, ci->progeny[k], 0);
+ for (int j = k + 1; j < 8; j++)
+ if (ci->progeny[j] != NULL)
+ DOSUB_PAIR1_BH(r, ci->progeny[k], ci->progeny[j], 0);
+ }
+ }
+
+ /* Otherwise, compute self-interaction. */
+ else {
+
+ /* Drift the cell to the current timestep if needed. */
+ if (!cell_are_bpart_drifted(ci, r->e)) error("Interacting undrifted cell.");
+
+ DOSELF1_BRANCH_BH(r, ci);
+ }
+
+ TIMER_TOC(TIMER_DOSUB_SELF_BH);
+}
diff --git a/src/scheduler.c b/src/scheduler.c
index 6dba5f1593c8bf050a260ded5061ebc54fc63310..af01e12f66540a039d5c485c266ac8f92c5d5a8d 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -539,6 +539,8 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
access the engine... */
const int with_feedback = (s->space->e->policy & engine_policy_feedback);
const int with_stars = (s->space->e->policy & engine_policy_stars);
+ const int with_black_holes =
+ (s->space->e->policy & engine_policy_black_holes);
/* Iterate on this task until we're done with it. */
int redo = 1;
@@ -549,14 +551,18 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
/* Is this a non-empty self-task? */
const int is_self =
(t->type == task_type_self) && (t->ci != NULL) &&
- ((t->ci->hydro.count > 0) || (with_stars && t->ci->stars.count > 0));
+ ((t->ci->hydro.count > 0) || (with_stars && t->ci->stars.count > 0) ||
+ (with_black_holes && t->ci->black_holes.count > 0));
/* Is this a non-empty pair-task? */
- const int is_pair =
- (t->type == task_type_pair) && (t->ci != NULL) && (t->cj != NULL) &&
- ((t->ci->hydro.count > 0) ||
- (with_feedback && t->ci->stars.count > 0)) &&
- ((t->cj->hydro.count > 0) || (with_feedback && t->cj->stars.count > 0));
+ const int is_pair = (t->type == task_type_pair) && (t->ci != NULL) &&
+ (t->cj != NULL) &&
+ ((t->ci->hydro.count > 0) ||
+ (with_feedback && t->ci->stars.count > 0) ||
+ (with_black_holes && t->ci->black_holes.count > 0)) &&
+ ((t->cj->hydro.count > 0) ||
+ (with_feedback && t->cj->stars.count > 0) ||
+ (with_black_holes && t->cj->black_holes.count > 0));
/* Empty task? */
if (!is_self && !is_pair) {
@@ -1175,7 +1181,7 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
const float scount_i = (t->ci != NULL) ? t->ci->stars.count : 0.f;
const float scount_j = (t->cj != NULL) ? t->cj->stars.count : 0.f;
const float bcount_i = (t->ci != NULL) ? t->ci->black_holes.count : 0.f;
- // const float bcount_j = (t->cj != NULL) ? t->cj->black_holes.count : 0.f;
+ const float bcount_j = (t->cj != NULL) ? t->cj->black_holes.count : 0.f;
switch (t->type) {
case task_type_sort:
@@ -1195,7 +1201,13 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
cost = 1.f * wscale * gcount_i;
else if (t->subtype == task_subtype_stars_density)
cost = 1.f * wscale * scount_i * count_i;
- else
+ else if (t->subtype == task_subtype_stars_feedback)
+ cost = 1.f * wscale * scount_i * count_i;
+ else if (t->subtype == task_subtype_bh_density)
+ cost = 1.f * wscale * bcount_i * count_i;
+ else if (t->subtype == task_subtype_bh_feedback)
+ cost = 1.f * wscale * bcount_i * count_i;
+ else // hydro loops
cost = 1.f * (wscale * count_i) * count_i;
break;
@@ -1205,7 +1217,9 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
cost = 3.f * (wscale * gcount_i) * gcount_j;
else
cost = 2.f * (wscale * gcount_i) * gcount_j;
- } else if (t->subtype == task_subtype_stars_density) {
+
+ } else if (t->subtype == task_subtype_stars_density ||
+ t->subtype == task_subtype_stars_feedback) {
if (t->ci->nodeID != nodeID)
cost = 3.f * wscale * count_i * scount_j * sid_scale[t->flags];
else if (t->cj->nodeID != nodeID)
@@ -1213,7 +1227,18 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
else
cost = 2.f * wscale * (scount_i * count_j + scount_j * count_i) *
sid_scale[t->flags];
- } else {
+
+ } else if (t->subtype == task_subtype_bh_density ||
+ t->subtype == task_subtype_bh_feedback) {
+ if (t->ci->nodeID != nodeID)
+ cost = 3.f * wscale * count_i * bcount_j * sid_scale[t->flags];
+ else if (t->cj->nodeID != nodeID)
+ cost = 3.f * wscale * bcount_i * count_j * sid_scale[t->flags];
+ else
+ cost = 2.f * wscale * (bcount_i * count_j + bcount_j * count_i) *
+ sid_scale[t->flags];
+
+ } else { // hydro loops
if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID)
cost = 3.f * (wscale * count_i) * count_j * sid_scale[t->flags];
else
@@ -1225,7 +1250,8 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
#ifdef SWIFT_DEBUG_CHECKS
if (t->flags < 0) error("Negative flag value!");
#endif
- if (t->subtype == task_subtype_stars_density) {
+ if (t->subtype == task_subtype_stars_density ||
+ t->subtype == task_subtype_stars_feedback) {
if (t->ci->nodeID != nodeID) {
cost = 3.f * (wscale * count_i) * scount_j * sid_scale[t->flags];
} else if (t->cj->nodeID != nodeID) {
@@ -1235,7 +1261,18 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
sid_scale[t->flags];
}
- } else {
+ } else if (t->subtype == task_subtype_bh_density ||
+ t->subtype == task_subtype_bh_feedback) {
+ if (t->ci->nodeID != nodeID) {
+ cost = 3.f * (wscale * count_i) * bcount_j * sid_scale[t->flags];
+ } else if (t->cj->nodeID != nodeID) {
+ cost = 3.f * (wscale * bcount_i) * count_j * sid_scale[t->flags];
+ } else {
+ cost = 2.f * wscale * (bcount_i * count_j + bcount_j * count_i) *
+ sid_scale[t->flags];
+ }
+
+ } else { // hydro loops
if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID) {
cost = 3.f * (wscale * count_i) * count_j * sid_scale[t->flags];
} else {
@@ -1247,6 +1284,12 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_sub_self:
if (t->subtype == task_subtype_stars_density) {
cost = 1.f * (wscale * scount_i) * count_i;
+ } else if (t->subtype == task_subtype_stars_feedback) {
+ cost = 1.f * (wscale * scount_i) * count_i;
+ } else if (t->subtype == task_subtype_bh_density) {
+ cost = 1.f * (wscale * bcount_i) * count_i;
+ } else if (t->subtype == task_subtype_bh_feedback) {
+ cost = 1.f * (wscale * bcount_i) * count_i;
} else {
cost = 1.f * (wscale * count_i) * count_i;
}
@@ -1260,6 +1303,9 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_stars_ghost:
if (t->ci == t->ci->hydro.super) cost = wscale * scount_i;
break;
+ case task_type_bh_ghost:
+ if (t->ci == t->ci->hydro.super) cost = wscale * bcount_i;
+ break;
case task_type_drift_part:
cost = wscale * count_i;
break;
@@ -1535,6 +1581,10 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
err = MPI_Irecv(t->ci->stars.parts, t->ci->stars.count,
spart_mpi_type, t->ci->nodeID, t->flags,
subtaskMPI_comms[t->subtype], &t->req);
+ } else if (t->subtype == task_subtype_bpart) {
+ err = MPI_Irecv(t->ci->black_holes.parts, t->ci->black_holes.count,
+ bpart_mpi_type, t->ci->nodeID, t->flags,
+ subtaskMPI_comms[t->subtype], &t->req);
} else if (t->subtype == task_subtype_multipole) {
t->buff = (struct gravity_tensors *)malloc(
sizeof(struct gravity_tensors) * t->ci->mpi.pcell_size);
@@ -1659,6 +1709,16 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
err = MPI_Issend(t->ci->stars.parts, t->ci->stars.count,
spart_mpi_type, t->cj->nodeID, t->flags,
subtaskMPI_comms[t->subtype], &t->req);
+ } else if (t->subtype == task_subtype_bpart) {
+ if ((t->ci->black_holes.count * sizeof(struct bpart)) >
+ s->mpi_message_limit)
+ err = MPI_Isend(t->ci->black_holes.parts, t->ci->black_holes.count,
+ bpart_mpi_type, t->cj->nodeID, t->flags,
+ subtaskMPI_comms[t->subtype], &t->req);
+ else
+ err = MPI_Issend(t->ci->black_holes.parts, t->ci->black_holes.count,
+ bpart_mpi_type, t->cj->nodeID, t->flags,
+ subtaskMPI_comms[t->subtype], &t->req);
} else if (t->subtype == task_subtype_multipole) {
t->buff = (struct gravity_tensors *)malloc(
sizeof(struct gravity_tensors) * t->ci->mpi.pcell_size);
diff --git a/src/space.c b/src/space.c
index 562766d9793c33743167f8c84490c9edbea1e661..1712e512e5b65bfbe71231fc2cbaba20447a5a8e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -235,6 +235,9 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->stars.ghost = NULL;
c->stars.density = NULL;
c->stars.feedback = NULL;
+ c->black_holes.ghost = NULL;
+ c->black_holes.density = NULL;
+ c->black_holes.feedback = NULL;
c->kick1 = NULL;
c->kick2 = NULL;
c->timestep = NULL;
@@ -245,6 +248,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->stars.stars_in = NULL;
c->stars.stars_out = NULL;
c->black_holes.drift = NULL;
+ c->black_holes.black_holes_in = NULL;
+ c->black_holes.black_holes_out = NULL;
c->grav.drift = NULL;
c->grav.drift_out = NULL;
c->hydro.cooling = NULL;
@@ -347,6 +352,7 @@ void space_regrid(struct space *s, int verbose) {
const size_t nr_parts = s->nr_parts;
const size_t nr_sparts = s->nr_sparts;
+ const size_t nr_bparts = s->nr_bparts;
const ticks tic = getticks();
const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
@@ -366,6 +372,9 @@ void space_regrid(struct space *s, int verbose) {
if (c->stars.h_max > h_max) {
h_max = c->stars.h_max;
}
+ if (c->black_holes.h_max > h_max) {
+ h_max = c->black_holes.h_max;
+ }
}
/* Can we instead use all the top-level cells? */
@@ -378,6 +387,9 @@ void space_regrid(struct space *s, int verbose) {
if (c->nodeID == engine_rank && c->stars.h_max > h_max) {
h_max = c->stars.h_max;
}
+ if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) {
+ h_max = c->black_holes.h_max;
+ }
}
/* Last option: run through the particles */
@@ -388,6 +400,9 @@ void space_regrid(struct space *s, int verbose) {
for (size_t k = 0; k < nr_sparts; k++) {
if (s->sparts[k].h > h_max) h_max = s->sparts[k].h;
}
+ for (size_t k = 0; k < nr_bparts; k++) {
+ if (s->bparts[k].h > h_max) h_max = s->bparts[k].h;
+ }
}
}
@@ -555,6 +570,8 @@ void space_regrid(struct space *s, int verbose) {
error("Failed to init spinlock for multipoles.");
if (lock_init(&s->cells_top[k].stars.lock) != 0)
error("Failed to init spinlock for stars.");
+ if (lock_init(&s->cells_top[k].black_holes.lock) != 0)
+ error("Failed to init spinlock for black holes.");
if (lock_init(&s->cells_top[k].stars.star_formation_lock) != 0)
error("Failed to init spinlock for star formation.");
}
@@ -584,6 +601,7 @@ void space_regrid(struct space *s, int verbose) {
c->hydro.ti_old_part = ti_current;
c->grav.ti_old_part = ti_current;
c->stars.ti_old_part = ti_current;
+ c->black_holes.ti_old_part = ti_current;
c->grav.ti_old_multipole = ti_current;
#ifdef WITH_MPI
c->mpi.tag = -1;
@@ -3331,6 +3349,7 @@ void space_split_recursive(struct space *s, struct cell *c,
/* Update the cell-wide properties */
h_max = max(h_max, cp->hydro.h_max);
stars_h_max = max(stars_h_max, cp->stars.h_max);
+ black_holes_h_max = max(black_holes_h_max, cp->black_holes.h_max);
ti_hydro_end_min = min(ti_hydro_end_min, cp->hydro.ti_end_min);
ti_hydro_end_max = max(ti_hydro_end_max, cp->hydro.ti_end_max);
ti_hydro_beg_max = max(ti_hydro_beg_max, cp->hydro.ti_beg_max);
@@ -3683,6 +3702,7 @@ void space_recycle(struct space *s, struct cell *c) {
/* Clear the cell. */
if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->grav.plock) != 0 ||
lock_destroy(&c->mlock) != 0 || lock_destroy(&c->stars.lock) != 0 ||
+ lock_destroy(&c->black_holes.lock) != 0 ||
lock_destroy(&c->stars.star_formation_lock))
error("Failed to destroy spinlocks.");
@@ -3733,6 +3753,7 @@ void space_recycle_list(struct space *s, struct cell *cell_list_begin,
/* Clear the cell. */
if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->grav.plock) != 0 ||
lock_destroy(&c->mlock) != 0 || lock_destroy(&c->stars.lock) != 0 ||
+ lock_destroy(&c->black_holes.lock) != 0 ||
lock_destroy(&c->stars.star_formation_lock))
error("Failed to destroy spinlocks.");
@@ -3833,6 +3854,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
lock_init(&cells[j]->grav.plock) != 0 ||
lock_init(&cells[j]->grav.mlock) != 0 ||
lock_init(&cells[j]->stars.lock) != 0 ||
+ lock_init(&cells[j]->black_holes.lock) != 0 ||
lock_init(&cells[j]->stars.star_formation_lock) != 0)
error("Failed to initialize cell spinlocks.");
}
@@ -4227,6 +4249,7 @@ void space_first_init_bparts_mapper(void *restrict map_data, int count,
struct bpart *restrict bp = (struct bpart *)map_data;
const struct space *restrict s = (struct space *)extra_data;
const struct engine *e = s->e;
+ const struct black_holes_props *props = e->black_holes_properties;
#ifdef SWIFT_DEBUG_CHECKS
const ptrdiff_t delta = bp - s->bparts;
@@ -4260,10 +4283,17 @@ void space_first_init_bparts_mapper(void *restrict map_data, int count,
#endif
}
+ /* Check that the smoothing lengths are non-zero */
+ for (int k = 0; k < count; k++) {
+ if (bp[k].h <= 0.)
+ error("Invalid value of smoothing length for bpart %lld h=%e", bp[k].id,
+ bp[k].h);
+ }
+
/* Initialise the rest */
for (int k = 0; k < count; k++) {
- black_holes_first_init_bpart(&bp[k]);
+ black_holes_first_init_bpart(&bp[k], props);
#ifdef SWIFT_DEBUG_CHECKS
if (bp[k].gpart && bp[k].gpart->id_or_neg_offset != -(k + delta))
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 4d1013c0f6dd13552cb533aa942f47dd588a469c..aaec0197100ab84003180f77ceb2570b5d342cbe 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -201,7 +201,7 @@ INLINE static double EOS_pressure(const double n_H,
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cooling The cooling data struct.
- * @param entropy_floor The entropy floor assumed in this run.
+ * @param entropy_floor_props The entropy floor assumed in this run.
*/
INLINE static int star_formation_is_star_forming(
const struct part* restrict p, const struct xpart* restrict xp,
@@ -376,6 +376,10 @@ INLINE static void star_formation_update_part_not_SFR(
* @param starform the star formation law properties to use.
* @param cosmo the cosmological parameters and properties.
* @param with_cosmology if we run with cosmology.
+ * @param phys_const the physical constants in internal units.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cooling The cooling data struct.
*/
INLINE static void star_formation_copy_properties(
const struct part* p, const struct xpart* xp, struct spart* sp,
diff --git a/src/swift.h b/src/swift.h
index cd7159e63d129e97f8aac5c8d426fd041bd298b0..562a1fd5f49e3599693940aafd4cb305743eed7e 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -25,6 +25,7 @@
/* Local headers. */
#include "active.h"
#include "atomic.h"
+#include "black_holes_properties.h"
#include "cache.h"
#include "cell.h"
#include "chemistry.h"
diff --git a/src/task.c b/src/task.c
index 005bd229ab8c88cfff1d800025091db6991a91cb..ad287157b892391bace086ead1ada1981d027ce9 100644
--- a/src/task.c
+++ b/src/task.c
@@ -88,15 +88,19 @@ const char *taskID_names[task_type_count] = {"none",
"stars_ghost_in",
"stars_ghost",
"stars_ghost_out",
- "stars_sort"};
+ "stars_sort",
+ "bh_in",
+ "bh_out",
+ "bh_ghost"};
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
- "none", "density", "gradient", "force",
- "limiter", "grav", "external_grav", "tend_part",
- "tend_gpart", "tend_spart", "tend_bpart", "xv",
- "rho", "gpart", "multipole", "spart",
- "stars_density", "stars_feedback"};
+ "none", "density", "gradient", "force",
+ "limiter", "grav", "external_grav", "tend_part",
+ "tend_gpart", "tend_spart", "tend_bpart", "xv",
+ "rho", "gpart", "multipole", "spart",
+ "stars_density", "stars_feedback", "bpart", "bh_density",
+ "bh_feedback"};
#ifdef WITH_MPI
/* MPI communicators for the subtypes. */
@@ -165,6 +169,11 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_spart;
break;
+ case task_type_drift_bpart:
+ case task_type_bh_ghost:
+ return task_action_bpart;
+ break;
+
case task_type_self:
case task_type_pair:
case task_type_sub_self:
@@ -183,6 +192,11 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_all;
break;
+ case task_subtype_bh_density:
+ case task_subtype_bh_feedback:
+ return task_action_all;
+ break;
+
case task_subtype_grav:
case task_subtype_external_grav:
return task_action_gpart;
@@ -706,6 +720,12 @@ void task_get_group_name(int type, int subtype, char *cluster) {
case task_subtype_stars_feedback:
strcpy(cluster, "StarsFeedback");
break;
+ case task_subtype_bh_density:
+ strcpy(cluster, "BHDensity");
+ break;
+ case task_subtype_bh_feedback:
+ strcpy(cluster, "BHFeedback");
+ break;
default:
strcpy(cluster, "None");
break;
diff --git a/src/task.h b/src/task.h
index b62b2d531b5cf4ae98331281249199388ee0070c..33c15372bb3af30c74e36fb572a2f12b8d29e6d4 100644
--- a/src/task.h
+++ b/src/task.h
@@ -84,6 +84,9 @@ enum task_types {
task_type_stars_ghost,
task_type_stars_ghost_out, /* Implicit */
task_type_stars_sort,
+ task_type_bh_in, /* Implicit */
+ task_type_bh_out, /* Implicit */
+ task_type_bh_ghost,
task_type_count
} __attribute__((packed));
@@ -109,6 +112,9 @@ enum task_subtypes {
task_subtype_spart,
task_subtype_stars_density,
task_subtype_stars_feedback,
+ task_subtype_bpart,
+ task_subtype_bh_density,
+ task_subtype_bh_feedback,
task_subtype_count
} __attribute__((packed));
@@ -120,6 +126,7 @@ enum task_actions {
task_action_part,
task_action_gpart,
task_action_spart,
+ task_action_bpart,
task_action_all,
task_action_multipole,
task_action_count
diff --git a/src/units.c b/src/units.c
index 1194735aafd80d51897e4c4cb3e4b0976478145a..55a1c9e1bdf62fce97b22189b4741be0fea2abb4 100644
--- a/src/units.c
+++ b/src/units.c
@@ -377,6 +377,7 @@ void units_get_base_unit_exponants_array(float baseUnitsExp[5],
break;
case UNIT_CONV_SFR:
+ case UNIT_CONV_MASS_PER_UNIT_TIME:
baseUnitsExp[UNIT_MASS] = 1.f;
baseUnitsExp[UNIT_TIME] = -1.f;
break;
diff --git a/src/units.h b/src/units.h
index 62669425e52c4e39800330a4150259856d8fc0bb..cc87696ccfdeb40b7cf4f23ccfcd853f8c8c1e56 100644
--- a/src/units.h
+++ b/src/units.h
@@ -96,7 +96,8 @@ enum unit_conversion_factor {
UNIT_CONV_VOLUME,
UNIT_CONV_INV_VOLUME,
UNIT_CONV_SFR,
- UNIT_CONV_SSFR
+ UNIT_CONV_SSFR,
+ UNIT_CONV_MASS_PER_UNIT_TIME
};
void units_init_cgs(struct unit_system*);
diff --git a/tools/plot_task_dependencies.py b/tools/plot_task_dependencies.py
index 14ba7c99f621c0c62c3c9cb8e8d3b0c78e491401..9a8b63ab315215c2973187a026fc1a5e182cd47b 100644
--- a/tools/plot_task_dependencies.py
+++ b/tools/plot_task_dependencies.py
@@ -136,6 +136,20 @@ def appendData(data):
return data[0]
+def taskIsBlackHoles(name):
+ """
+ Does the task concern black holes?
+
+ Parameters
+ ----------
+
+ name: str
+ Task name
+ """
+ if "bh" in name or "bpart" in name:
+ return True
+ return False
+
def taskIsStars(name):
"""
Does the task concern stars?
@@ -150,7 +164,6 @@ def taskIsStars(name):
return True
return False
-
def taskIsHydro(name):
"""
Does the task concern the hydro?
@@ -163,7 +176,7 @@ def taskIsHydro(name):
"""
if "_part" in name:
return True
- if "density" in name and "stars" not in name:
+ if "density" in name and "stars" not in name and "bh" not in name:
return True
if "rho" in name:
return True
@@ -291,6 +304,9 @@ def writeTask(f, name, implicit, mpi, with_calls):
if mpi:
txt += "shape=diamond,"
+ if taskIsBlackHoles(name):
+ txt += "color=forestgreen,"
+
if taskIsStars(name):
txt += "color=darkorange1,"
diff --git a/tools/task_plots/analyse_tasks.py b/tools/task_plots/analyse_tasks.py
index e938d0610ab878efbd6463909bbf75fe2ec60bc2..798efdfdaa7fb982de34917488ab372b2f0d350a 100755
--- a/tools/task_plots/analyse_tasks.py
+++ b/tools/task_plots/analyse_tasks.py
@@ -77,6 +77,7 @@ TASKTYPES = [
"extra_ghost",
"drift_part",
"drift_spart",
+ "drift_bpart",
"drift_gpart",
"drift_gpart_out",
"hydro_end_force",
@@ -103,6 +104,9 @@ TASKTYPES = [
"stars_ghost",
"stars_ghost_out",
"stars_sort",
+ "bh_in",
+ "bh_out",
+ "bh_ghost",
"count",
]
@@ -117,6 +121,7 @@ SUBTYPES = [
"tend_part",
"tend_gpart",
"tend_spart",
+ "tend_bpart",
"xv",
"rho",
"gpart",
@@ -124,6 +129,9 @@ SUBTYPES = [
"spart",
"stars_density",
"stars_feedback",
+ "bpart",
+ "bh_density",
+ "bh_feedback",
"count",
]
diff --git a/tools/task_plots/plot_tasks.py b/tools/task_plots/plot_tasks.py
index b692ee963784b9487f5e19ca7274c0ecaa9b7a89..a76fcbb3ac4092daafa2c2334a810dbada49c2a4 100755
--- a/tools/task_plots/plot_tasks.py
+++ b/tools/task_plots/plot_tasks.py
@@ -162,6 +162,7 @@ TASKTYPES = [
"extra_ghost",
"drift_part",
"drift_spart",
+ "drift_bpart",
"drift_gpart",
"drift_gpart_out",
"hydro_end_force",
@@ -188,6 +189,9 @@ TASKTYPES = [
"stars_ghost",
"stars_ghost_out",
"stars_sort",
+ "bh_in",
+ "bh_out",
+ "bh_ghost",
"count",
]
@@ -202,6 +206,7 @@ SUBTYPES = [
"tend_part",
"tend_gpart",
"tend_spart",
+ "tend_bpart",
"xv",
"rho",
"gpart",
@@ -209,6 +214,9 @@ SUBTYPES = [
"spart",
"stars_density",
"stars_feedback",
+ "bpart",
+ "bh_density",
+ "bh_feedback",
"count",
]
@@ -242,10 +250,14 @@ FULLTYPES = [
"send/tend_gpart",
"recv/tend_spart",
"send/tend_spart",
+ "recv/tend_bpart",
+ "send/tend_bpart",
"recv/gpart",
"send/gpart",
"recv/spart",
"send/spart",
+ "recv/bpart",
+ "send/bpart",
"self/stars_density",
"pair/stars_density",
"sub_self/stars_density",
@@ -254,6 +266,14 @@ FULLTYPES = [
"pair/stars_feedback",
"sub_self/stars_feedback",
"sub_pair/stars_feedback",
+ "self/bh_density",
+ "pair/bh_density",
+ "sub_self/bh_density",
+ "sub_pair/bh_density",
+ "self/bh_feedback",
+ "pair/bh_feedback",
+ "sub_self/bh_feedback",
+ "sub_pair/bh_feedback",
]
# A number of colours for the various types. Recycled when there are