diff --git a/README b/README
index 7060589401d80e205733fb5770f258708263d966..2e26f9d997f7396109638a5f44369922ebc4b275 100644
--- a/README
+++ b/README
@@ -35,6 +35,7 @@ Parameters:
-M, --multipole-reconstruction Reconstruct the multipoles every time-step.
-s, --hydro Run with hydrodynamics.
-S, --stars Run with stars.
+ -B, --black-holes Run with black holes.
-x, --velociraptor Run with structure finding.
--limiter Run with time-step limiter.
diff --git a/configure.ac b/configure.ac
index 50619f589b9377001538f3c62870c72249b3f102..59a9826767fd0ad82ceda23657b143c377e14472 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1330,6 +1330,7 @@ case "$with_subgrid" in
with_subgrid_stars=GEAR
with_subgrid_star_formation=GEAR
with_subgrid_feedback=thermal
+ with_subgrid_black_holes=none
;;
EAGLE)
with_subgrid_cooling=EAGLE
@@ -1339,6 +1340,7 @@ case "$with_subgrid" in
with_subgrid_stars=EAGLE
with_subgrid_star_formation=EAGLE
with_subgrid_feedback=none
+ with_subgrid_black_holes=none
;;
*)
AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
@@ -1765,6 +1767,33 @@ case "$with_feedback" in
;;
esac
+# Stellar model.
+AC_ARG_WITH([black-holes],
+ [AS_HELP_STRING([--with-black-holes=<model>],
+ [Black holes model to use @<:@none, default: none@:>@]
+ )],
+ [with_black_holes="$withval"],
+ [with_black_holes="none"]
+)
+
+if test "$with_subgrid" != "none"; then
+ if test "$with_black_holes" != "none"; then
+ AC_MSG_ERROR([Cannot provide with-subgrid and with-black-holes together])
+ else
+ with_black_holes="$with_subgrid_black_holes"
+ fi
+fi
+
+case "$with_black_holes" in
+ none)
+ AC_DEFINE([BLACK_HOLES_NONE], [1], [No black hole model])
+ ;;
+ *)
+ AC_MSG_ERROR([Unknown stellar model: $with_black_holes])
+ ;;
+esac
+
+
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
@@ -1962,7 +1991,8 @@ AC_MSG_RESULT([
Tracers : $with_tracers
Stellar model : $with_stars
Star formation model : $with_star_formation
- Feedback model : $with_feedback
+ Star feedback model : $with_feedback
+ Black holes model : $with_black_holes
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
diff --git a/examples/main.c b/examples/main.c
index 52eed321a5cd80f581eeea11f315bc766457b4ab..ba0b11eae89431ac663a5f1a6cb6baa45306bfaa 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -101,6 +101,7 @@ int main(int argc, char *argv[]) {
struct phys_const prog_const;
struct space s;
struct spart *sparts = NULL;
+ struct bpart *bparts = NULL;
struct unit_system us;
int nr_nodes = 1, myrank = 0;
@@ -156,6 +157,7 @@ int main(int argc, char *argv[]) {
int with_stars = 0;
int with_star_formation = 0;
int with_feedback = 0;
+ int with_black_holes = 0;
int with_limiter = 0;
int with_fp_exceptions = 0;
int with_drift_all = 0;
@@ -204,6 +206,8 @@ int main(int argc, char *argv[]) {
OPT_BOOLEAN('s', "hydro", &with_hydro, "Run with hydrodynamics.", NULL, 0,
0),
OPT_BOOLEAN('S', "stars", &with_stars, "Run with stars.", NULL, 0, 0),
+ OPT_BOOLEAN('B', "black-holes", &with_black_holes,
+ "Run with black holes.", NULL, 0, 0),
OPT_BOOLEAN('x', "velociraptor", &with_structure_finding,
"Run with structure finding.", NULL, 0, 0),
OPT_BOOLEAN(0, "limiter", &with_limiter, "Run with time-step limiter.",
@@ -344,6 +348,16 @@ int main(int argc, char *argv[]) {
return 1;
}
+ if (with_black_holes && !with_external_gravity && !with_self_gravity) {
+ if (myrank == 0) {
+ argparse_usage(&argparse);
+ printf(
+ "\nError: Cannot process black holes without gravity, -g or -G "
+ "must be chosen.\n");
+ }
+ return 1;
+ }
+
if (!with_stars && with_star_formation) {
if (myrank == 0) {
argparse_usage(&argparse);
@@ -448,6 +462,7 @@ int main(int argc, char *argv[]) {
message("sizeof(part) is %4zi bytes.", sizeof(struct part));
message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
+ message("sizeof(bpart) is %4zi bytes.", sizeof(struct bpart));
message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
message("sizeof(grav_tensor) is %4zi bytes.", sizeof(struct grav_tensor));
@@ -484,6 +499,7 @@ int main(int argc, char *argv[]) {
if (with_star_formation && with_feedback)
error("Can't run with star formation and feedback over MPI (yet)");
if (with_limiter) error("Can't run with time-step limiter over MPI (yet)");
+ if (with_black_holes) error("Can't run with black holes over MPI (yet)");
#endif
/* Temporary early aborts for modes not supported with hand-vec. */
@@ -777,7 +793,7 @@ int main(int argc, char *argv[]) {
fflush(stdout);
/* Get ready to read particles of all kinds */
- size_t Ngas = 0, Ngpart = 0, Nspart = 0;
+ size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
double dim[3] = {0., 0., 0.};
if (myrank == 0) clocks_gettime(&tic);
#if defined(HAVE_HDF5)
@@ -798,11 +814,11 @@ int main(int argc, char *argv[]) {
dry_run);
#endif
#else
- read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
- &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
- (with_external_gravity || with_self_gravity), with_stars,
- cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
- dry_run);
+ read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
+ &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
+ with_hydro, (with_external_gravity || with_self_gravity),
+ with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
+ cosmo.h, cosmo.a, nr_threads, dry_run);
#endif
#endif
if (myrank == 0) {
@@ -825,35 +841,39 @@ int main(int argc, char *argv[]) {
/* Check that the other links are correctly set */
if (!dry_run)
- part_verify_links(parts, gparts, sparts, Ngas, Ngpart, Nspart, 1);
+ part_verify_links(parts, gparts, sparts, bparts, Ngas, Ngpart, Nspart,
+ Nbpart, 1);
#endif
/* Get the total number of particles across all nodes. */
- long long N_total[3] = {0, 0, 0};
+ long long N_total[4] = {0, 0, 0, 0};
#if defined(WITH_MPI)
- long long N_long[3] = {Ngas, Ngpart, Nspart};
+ long long N_long[4] = {Ngas, Ngpart, Nspart, Nbpart};
MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
MPI_COMM_WORLD);
#else
N_total[0] = Ngas;
N_total[1] = Ngpart;
N_total[2] = Nspart;
+ N_total[3] = Nbpart;
#endif
if (myrank == 0)
message(
- "Read %lld gas particles, %lld stars particles and %lld gparts from "
- "the ICs.",
- N_total[0], N_total[2], N_total[1]);
+ "Read %lld gas particles, %lld stars particles, %lld black hole "
+ "particles"
+ " and %lld gparts from the ICs.",
+ N_total[0], N_total[2], N_total[3], N_total[1]);
/* Verify that the fields to dump actually exist */
if (myrank == 0) io_check_output_fields(params, N_total);
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
- space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
- Nspart, periodic, replicate, generate_gas_in_ics, with_hydro,
- with_self_gravity, with_star_formation, talking, dry_run);
+ space_init(&s, params, &cosmo, dim, parts, gparts, sparts, bparts, Ngas,
+ Ngpart, Nspart, Nbpart, periodic, replicate, generate_gas_in_ics,
+ with_hydro, with_self_gravity, with_star_formation, talking,
+ dry_run);
if (myrank == 0) {
clocks_gettime(&toc);
@@ -885,12 +905,14 @@ int main(int argc, char *argv[]) {
N_long[0] = s.nr_parts;
N_long[1] = s.nr_gparts;
N_long[2] = s.nr_sparts;
- MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
+ N_long[3] = s.nr_bparts;
+ MPI_Allreduce(&N_long, &N_total, 4, MPI_LONG_LONG_INT, MPI_SUM,
MPI_COMM_WORLD);
#else
N_total[0] = s.nr_parts;
N_total[1] = s.nr_gparts;
N_total[2] = s.nr_sparts;
+ N_total[3] = s.nr_bparts;
#endif
/* Say a few nice things about the space we just created. */
@@ -903,6 +925,7 @@ int main(int argc, char *argv[]) {
message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
+ message("%zi bparts in %i cells.", s.nr_bparts, s.tot_cells);
message("maximum depth is %d.", s.maxdepth);
fflush(stdout);
}
@@ -949,6 +972,7 @@ int main(int argc, char *argv[]) {
if (with_stars) engine_policies |= engine_policy_stars;
if (with_star_formation) engine_policies |= engine_policy_star_formation;
if (with_feedback) engine_policies |= engine_policy_feedback;
+ if (with_black_holes) engine_policies |= engine_policy_black_holes;
if (with_structure_finding)
engine_policies |= engine_policy_structure_finding;
@@ -971,11 +995,12 @@ int main(int argc, char *argv[]) {
/* Get some info to the user. */
if (myrank == 0) {
- long long N_DM = N_total[1] - N_total[2] - N_total[0];
+ long long N_DM = N_total[1] - N_total[2] - N_total[3] - N_total[0];
message(
- "Running on %lld gas particles, %lld stars particles and %lld DM "
- "particles (%lld gravity particles)",
- N_total[0], N_total[2], N_total[1] > 0 ? N_DM : 0, N_total[1]);
+ "Running on %lld gas particles, %lld stars particles %lld black "
+ "hole particles and %lld DM particles (%lld gravity particles)",
+ N_total[0], N_total[2], N_total[3], N_total[1] > 0 ? N_DM : 0,
+ N_total[1]);
message(
"from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
"task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
diff --git a/src/active.h b/src/active.h
index 6466cd314fdc18ad324bf01a1ff4e73e214e35d5..cef22d3bf3bf94866c8c852500a19bd237abd81e 100644
--- a/src/active.h
+++ b/src/active.h
@@ -97,6 +97,29 @@ __attribute__((always_inline)) INLINE static int cell_are_spart_drifted(
return (c->stars.ti_old_part == e->ti_current);
}
+/**
+ * @brief Check that the #bpart in a #cell have been drifted to the current
+ * time.
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell has been drifted to the current time, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_are_bpart_drifted(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->black_holes.ti_old_part > e->ti_current)
+ error(
+ "Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
+ "and e->ti_current=%lld (t=%e)",
+ c->black_holes.ti_old_part, c->black_holes.ti_old_part * e->time_base,
+ e->ti_current, e->ti_current * e->time_base);
+#endif
+
+ return (c->black_holes.ti_old_part == e->ti_current);
+}
+
/* Are cells / particles active for regular tasks ? */
/**
@@ -220,6 +243,28 @@ __attribute__((always_inline)) INLINE static int cell_is_active_stars(
return (c->stars.ti_end_min == e->ti_current);
}
+/**
+ * @brief Does a cell contain any b-particle finishing their time-step now ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell contains at least an active particle, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_active_black_holes(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->black_holes.ti_end_min < e->ti_current)
+ error(
+ "cell in an impossible time-zone! c->ti_end_min=%lld (t=%e) and "
+ "e->ti_current=%lld (t=%e, a=%e)",
+ c->black_holes.ti_end_min, c->black_holes.ti_end_min * e->time_base, e->ti_current,
+ e->ti_current * e->time_base, e->cosmology->a);
+#endif
+
+ return (c->black_holes.ti_end_min == e->ti_current);
+}
+
/**
* @brief Is this particle finishing its time-step now ?
*
@@ -308,6 +353,33 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
return (spart_bin <= max_active_bin);
}
+/**
+ * @brief Is this b-particle finishing its time-step now ?
+ *
+ * @param bp The #bpart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #bpart is active, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int bpart_is_active(
+ const struct bpart *bp, const struct engine *e) {
+
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t bpart_bin = bp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_end = get_integer_time_end(ti_current, bp->time_bin);
+
+ if (ti_end < ti_current)
+ error(
+ "s-particle in an impossible time-zone! bp->ti_end=%lld "
+ "e->ti_current=%lld",
+ ti_end, ti_current);
+#endif
+
+ return (bpart_bin <= max_active_bin);
+}
+
/**
* @brief Has this particle been inhibited?
*
@@ -344,6 +416,18 @@ __attribute__((always_inline)) INLINE static int spart_is_inhibited(
return sp->time_bin == time_bin_inhibited;
}
+/**
+ * @brief Has this black hole particle been inhibited?
+ *
+ * @param sp The #bpart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #part is inhibited, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int bpart_is_inhibited(
+ const struct bpart *bp, const struct engine *e) {
+ return bp->time_bin == time_bin_inhibited;
+}
+
/* Are cells / particles active for kick1 tasks ? */
/**
@@ -496,4 +580,32 @@ __attribute__((always_inline)) INLINE static int spart_is_starting(
return (spart_bin <= max_active_bin);
}
+/**
+ * @brief Is this b-particle starting its time-step now ?
+ *
+ * @param sp The #bpart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #bpart is active, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int bpart_is_starting(
+ const struct bpart *bp, const struct engine *e) {
+
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t bpart_bin = bp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, bp->time_bin);
+
+ if (ti_beg > ti_current)
+ error(
+ "s-particle in an impossible time-zone! bp->ti_beg=%lld "
+ "e->ti_current=%lld",
+ ti_beg, ti_current);
+#endif
+
+ return (bpart_bin <= max_active_bin);
+}
+
#endif /* SWIFT_ACTIVE_H */
diff --git a/src/black_holes.h b/src/black_holes.h
new file mode 100644
index 0000000000000000000000000000000000000000..0decc3fbf5eddd30cc5b5a15b73d0a4c2447e5ee
--- /dev/null
+++ b/src/black_holes.h
@@ -0,0 +1,32 @@
+/*******************************************************************************
+ * 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_BLACK_HOLES_H
+#define SWIFT_BLACK_HOLES_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct star model */
+#if defined(BLACK_HOLES_NONE)
+#include "./black_holes/Default/black_holes.h"
+#else
+#error "Invalid choice of star model"
+#endif
+
+#endif
diff --git a/src/black_holes/Default/black_holes.h b/src/black_holes/Default/black_holes.h
new file mode 100644
index 0000000000000000000000000000000000000000..2617cb1e9f14209c0ec4f568cbe91e12065522ab
--- /dev/null
+++ b/src/black_holes/Default/black_holes.h
@@ -0,0 +1,163 @@
+/*******************************************************************************
+ * 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_DEFAULT_BLACK_HOLES_H
+#define SWIFT_DEFAULT_BLACK_HOLES_H
+
+#include <float.h>
+#include "dimension.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+
+/**
+ * @brief Computes the gravity time-step of a given star particle.
+ *
+ * @param sp 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 s-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
+ struct bpart* bp) {
+
+ bp->time_bin = 0;
+}
+
+/**
+ * @brief Prepares a s-particle for its interactions
+ *
+ * @param sp 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)
+ sp->ids_ngbs_density[i] = -1;
+ sp->num_ngb_density = 0;
+#endif
+
+ bp->density.wcount = 0.f;
+ bp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param sp 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 sp The particle.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_reset_predicted_values(
+ struct bpart* restrict bp) {}
+
+/**
+ * @brief Finishes the calculation of (non-gravity) forces acting on stars
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_end_feedback(
+ struct bpart* bp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param sp 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 stars
+ *
+ * @param sp 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;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #spart has 0
+ * ngbs.
+ *
+ * @param sp 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 Reset acceleration fields of a particle
+ *
+ * This is the equivalent of hydro_reset_acceleration.
+ * We do not compute the acceleration on star, therefore no need to use it.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_reset_feedback(
+ struct bpart* restrict bp) {
+
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ p->ids_ngbs_force[i] = -1;
+ p->num_ngb_force = 0;
+#endif
+}
+
+#endif /* SWIFT_DEFAULT_BLACK_HOLES_H */
diff --git a/src/black_holes/Default/black_holes_io.h b/src/black_holes/Default/black_holes_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..a92086b2c99bed7138f4a298473377bc3af7c026
--- /dev/null
+++ b/src/black_holes/Default/black_holes_io.h
@@ -0,0 +1,96 @@
+/*******************************************************************************
+ * 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_DEFAULT_BLACK_HOLES_IO_H
+#define SWIFT_DEFAULT_BLACK_HOLES_IO_H
+
+#include "io_properties.h"
+#include "black_holes_part.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 = 5;
+
+ /* 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);
+}
+
+/**
+ * @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 = 5;
+
+ /* 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);
+
+#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_DEFAULT_BLACK_HOLES_IO_H */
diff --git a/src/black_holes/Default/black_holes_part.h b/src/black_holes/Default/black_holes_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..360eb45d9e2ccb14e8ff55b6d286e7bb91ef89cc
--- /dev/null
+++ b/src/black_holes/Default/black_holes_part.h
@@ -0,0 +1,92 @@
+/*******************************************************************************
+ * 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_DEFAULT_BLACK_HOLE_PART_H
+#define SWIFT_DEFAULT_BLACK_HOLE_PART_H
+
+/* Some standard headers. */
+#include <stdlib.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;
+
+#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_DEFAULT_BLACK_HOLE_PART_H */
diff --git a/src/black_holes_io.h b/src/black_holes_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..05e66645b182bbd26907b37d5c5560aa5e8b27a7
--- /dev/null
+++ b/src/black_holes_io.h
@@ -0,0 +1,31 @@
+/*******************************************************************************
+ * 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_BLACK_HOLES_IO_H
+#define SWIFT_BLACK_HOLES_IO_H
+
+#include "../config.h"
+
+/* Load the correct star type */
+#if defined(BLACK_HOLES_NONE)
+#include "./black_holes/Default/black_holes_io.h"
+#else
+#error "Invalid choice of star model"
+#endif
+
+#endif /* SWIFT_BLACK_HOLES_IO_H */
diff --git a/src/cell.c b/src/cell.c
index c4b738e112cbe2a11c2de0189aa04152476eacc0..07d896db2103a6ed7549cc7737c5377c9e577680 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -49,6 +49,7 @@
/* Local headers. */
#include "active.h"
#include "atomic.h"
+#include "black_holes.h"
#include "chemistry.h"
#include "drift.h"
#include "engine.h"
@@ -808,6 +809,74 @@ int cell_unpack_end_step_stars(struct cell *restrict c,
#endif
}
+/**
+ * @brief Pack the time information of the given cell and all it's sub-cells.
+ *
+ * @param c The #cell.
+ * @param pcells (output) The end-of-timestep information we pack into
+ *
+ * @return The number of packed cells.
+ */
+int cell_pack_end_step_black_holes(struct cell *restrict c,
+ struct pcell_step_black_holes *restrict pcells) {
+
+#ifdef WITH_MPI
+
+ /* Pack this cell's data. */
+ pcells[0].ti_end_min = c->black_holes.ti_end_min;
+ pcells[0].ti_end_max = c->black_holes.ti_end_max;
+ pcells[0].dx_max_part = c->black_holes.dx_max_part;
+
+ /* Fill in the progeny, depth-first recursion. */
+ int count = 1;
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ count += cell_pack_end_step_black_holes(c->progeny[k], &pcells[count]);
+ }
+
+ /* Return the number of packed values. */
+ return count;
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+ return 0;
+#endif
+}
+
+/**
+ * @brief Unpack the time information of a given cell and its sub-cells.
+ *
+ * @param c The #cell
+ * @param pcells The end-of-timestep information to unpack
+ *
+ * @return The number of cells created.
+ */
+int cell_unpack_end_step_black_holes(struct cell *restrict c,
+ struct pcell_step_black_holes *restrict pcells) {
+
+#ifdef WITH_MPI
+
+ /* Unpack this cell's data. */
+ c->black_holes.ti_end_min = pcells[0].ti_end_min;
+ c->black_holes.ti_end_max = pcells[0].ti_end_max;
+ c->black_holes.dx_max_part = pcells[0].dx_max_part;
+
+ /* Fill in the progeny, depth-first recursion. */
+ int count = 1;
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ count += cell_unpack_end_step_black_holes(c->progeny[k], &pcells[count]);
+ }
+
+ /* Return the number of packed values. */
+ return count;
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+ return 0;
+#endif
+}
+
/**
* @brief Pack the multipole information of the given cell and all it's
* sub-cells.
@@ -1213,23 +1282,30 @@ void cell_sunlocktree(struct cell *c) {
* space's parts array, i.e. c->hydro.parts - s->parts.
* @param sparts_offset Offset of the cell sparts array relative to the
* space's sparts array, i.e. c->stars.parts - s->stars.parts.
+ * @param bparts_offset Offset of the cell bparts array relative to the
+ * space's bparts array, i.e. c->black_holes.parts -
+ * s->black_holes.parts.
* @param buff A buffer with at least max(c->hydro.count, c->grav.count)
* entries, used for sorting indices.
* @param sbuff A buffer with at least max(c->stars.count, c->grav.count)
* entries, used for sorting indices for the sparts.
+ * @param bbuff A buffer with at least max(c->black_holes.count, c->grav.count)
+ * entries, used for sorting indices for the sparts.
* @param gbuff A buffer with at least max(c->hydro.count, c->grav.count)
* entries, used for sorting indices for the gparts.
*/
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
- struct cell_buff *buff, struct cell_buff *sbuff,
+ ptrdiff_t bparts_offset, struct cell_buff *buff,
+ struct cell_buff *sbuff, struct cell_buff *bbuff,
struct cell_buff *gbuff) {
const int count = c->hydro.count, gcount = c->grav.count,
- scount = c->stars.count;
+ scount = c->stars.count, bcount = c->black_holes.count;
struct part *parts = c->hydro.parts;
struct xpart *xparts = c->hydro.xparts;
struct gpart *gparts = c->grav.parts;
struct spart *sparts = c->stars.parts;
+ struct bpart *bparts = c->black_holes.parts;
const double pivot[3] = {c->loc[0] + c->width[0] / 2,
c->loc[1] + c->width[1] / 2,
c->loc[2] + c->width[2] / 2};
@@ -1253,6 +1329,11 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
sbuff[k].x[2] != sparts[k].x[2])
error("Inconsistent sbuff contents.");
}
+ for (int k = 0; k < bcount; k++) {
+ if (bbuff[k].x[0] != bparts[k].x[0] || bbuff[k].x[1] != bparts[k].x[1] ||
+ bbuff[k].x[2] != bparts[k].x[2])
+ error("Inconsistent bbuff contents.");
+ }
#endif /* SWIFT_DEBUG_CHECKS */
/* Fill the buffer with the indices. */
@@ -1427,6 +1508,60 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
c->progeny[k]->stars.parts = &c->stars.parts[bucket_offset[k]];
}
+ /* Now do the same song and dance for the bparts. */
+ for (int k = 0; k < 8; k++) bucket_count[k] = 0;
+
+ /* Fill the buffer with the indices. */
+ for (int k = 0; k < bcount; k++) {
+ const int bid = (bbuff[k].x[0] > pivot[0]) * 4 +
+ (bbuff[k].x[1] > pivot[1]) * 2 + (bbuff[k].x[2] > pivot[2]);
+ bucket_count[bid]++;
+ bbuff[k].ind = bid;
+ }
+
+ /* Set the buffer offsets. */
+ bucket_offset[0] = 0;
+ for (int k = 1; k <= 8; k++) {
+ bucket_offset[k] = bucket_offset[k - 1] + bucket_count[k - 1];
+ bucket_count[k - 1] = 0;
+ }
+
+ /* Run through the buckets, and swap particles to their correct spot. */
+ for (int bucket = 0; bucket < 8; bucket++) {
+ for (int k = bucket_offset[bucket] + bucket_count[bucket];
+ k < bucket_offset[bucket + 1]; k++) {
+ int bid = bbuff[k].ind;
+ if (bid != bucket) {
+ struct bpart bpart = bparts[k];
+ struct cell_buff temp_buff = bbuff[k];
+ while (bid != bucket) {
+ int j = bucket_offset[bid] + bucket_count[bid]++;
+ while (bbuff[j].ind == bid) {
+ j++;
+ bucket_count[bid]++;
+ }
+ memswap(&bparts[j], &bpart, sizeof(struct bpart));
+ memswap(&bbuff[j], &temp_buff, sizeof(struct cell_buff));
+ if (bparts[j].gpart)
+ bparts[j].gpart->id_or_neg_offset = -(j + bparts_offset);
+ bid = temp_buff.ind;
+ }
+ bparts[k] = bpart;
+ bbuff[k] = temp_buff;
+ if (bparts[k].gpart)
+ bparts[k].gpart->id_or_neg_offset = -(k + bparts_offset);
+ }
+ bucket_count[bid]++;
+ }
+ }
+
+ /* Store the counts and offsets. */
+ for (int k = 0; k < 8; k++) {
+ c->progeny[k]->black_holes.count = bucket_count[k];
+ c->progeny[k]->black_holes.count_total = c->progeny[k]->black_holes.count;
+ c->progeny[k]->black_holes.parts = &c->black_holes.parts[bucket_offset[k]];
+ }
+
/* Finally, do the same song and dance for the gparts. */
for (int k = 0; k < 8; k++) bucket_count[k] = 0;
@@ -1467,6 +1602,9 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
} else if (gparts[j].type == swift_type_stars) {
sparts[-gparts[j].id_or_neg_offset - sparts_offset].gpart =
&gparts[j];
+ } else if (gparts[j].type == swift_type_black_hole) {
+ bparts[-gparts[j].id_or_neg_offset - bparts_offset].gpart =
+ &gparts[j];
}
bid = temp_buff.ind;
}
@@ -1477,6 +1615,9 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
} else if (gparts[k].type == swift_type_stars) {
sparts[-gparts[k].id_or_neg_offset - sparts_offset].gpart =
&gparts[k];
+ } else if (gparts[k].type == swift_type_black_hole) {
+ bparts[-gparts[k].id_or_neg_offset - bparts_offset].gpart =
+ &gparts[k];
}
}
bucket_count[bid]++;
@@ -4381,6 +4522,164 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
c->stars.do_sub_drift = 0;
}
+/**
+ * @brief Recursively drifts the #bpart in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ * @param force Drift the particles irrespective of the #cell flags.
+ */
+void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
+
+ const int periodic = e->s->periodic;
+ const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const float black_holes_h_max = e->hydro_properties->h_max;
+ const float black_holes_h_min = e->hydro_properties->h_min;
+ const integertime_t ti_old_bpart = c->black_holes.ti_old_part;
+ const integertime_t ti_current = e->ti_current;
+ struct bpart *const bparts = c->black_holes.parts;
+
+ float dx_max = 0.f, dx2_max = 0.f;
+ float cell_h_max = 0.f;
+
+ /* Drift irrespective of cell flags? */
+ force |= c->black_holes.do_drift;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that we only drift local cells. */
+ if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_bpart) error("Attempt to drift to the past");
+#endif
+
+ /* Early abort? */
+ if (c->black_holes.count == 0) {
+
+ /* Clear the drift flags. */
+ c->black_holes.do_drift = 0;
+ c->black_holes.do_sub_drift = 0;
+
+ /* Update the time of the last drift */
+ c->black_holes.ti_old_part = ti_current;
+
+ return;
+ }
+
+ /* Ok, we have some particles somewhere in the hierarchy to drift */
+
+ /* Are we not in a leaf ? */
+ if (c->split && (force || c->black_holes.do_sub_drift)) {
+
+ /* Loop over the progeny and collect their data. */
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+
+ /* Recurse */
+ cell_drift_bpart(cp, e, force);
+
+ /* Update */
+ dx_max = max(dx_max, cp->black_holes.dx_max_part);
+ cell_h_max = max(cell_h_max, cp->black_holes.h_max);
+ }
+ }
+
+ /* Store the values */
+ c->black_holes.h_max = cell_h_max;
+ c->black_holes.dx_max_part = dx_max;
+
+ /* Update the time of the last drift */
+ c->black_holes.ti_old_part = ti_current;
+
+ } else if (!c->split && force && ti_current > ti_old_bpart) {
+
+ /* Drift from the last time the cell was drifted to the current time */
+ double dt_drift;
+ if (with_cosmology) {
+ dt_drift =
+ cosmology_get_drift_factor(e->cosmology, ti_old_bpart, ti_current);
+ } else {
+ dt_drift = (ti_current - ti_old_bpart) * e->time_base;
+ }
+
+ /* Loop over all the star particles in the cell */
+ const size_t nr_bparts = c->black_holes.count;
+ for (size_t k = 0; k < nr_bparts; k++) {
+
+ /* Get a handle on the bpart. */
+ struct bpart *const bp = &bparts[k];
+
+ /* Ignore inhibited particles */
+ if (bpart_is_inhibited(bp, e)) continue;
+
+ /* Drift... */
+ drift_bpart(bp, dt_drift, ti_old_bpart, ti_current);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure the particle does not drift by more than a box length. */
+ if (fabs(bp->v[0] * dt_drift) > e->s->dim[0] ||
+ fabs(bp->v[1] * dt_drift) > e->s->dim[1] ||
+ fabs(bp->v[2] * dt_drift) > e->s->dim[2]) {
+ error("Particle drifts by more than a box length!");
+ }
+#endif
+
+ /* In non-periodic BC runs, remove particles that crossed the border */
+ if (!periodic) {
+
+ /* Did the particle leave the box? */
+ if ((bp->x[0] > dim[0]) || (bp->x[0] < 0.) || // x
+ (bp->x[1] > dim[1]) || (bp->x[1] < 0.) || // y
+ (bp->x[2] > dim[2]) || (bp->x[2] < 0.)) { // z
+
+ /* Remove the particle entirely */
+ struct gpart *gp = bp->gpart;
+ cell_remove_bpart(e, c, bp);
+
+ /* and it's gravity friend */
+ cell_remove_gpart(e, c, gp);
+
+ continue;
+ }
+ }
+
+ /* Limit h to within the allowed range */
+ bp->h = min(bp->h, black_holes_h_max);
+ bp->h = max(bp->h, black_holes_h_min);
+
+ /* Compute (square of) motion since last cell construction */
+ const float dx2 = bp->x_diff[0] * bp->x_diff[0] +
+ bp->x_diff[1] * bp->x_diff[1] +
+ bp->x_diff[2] * bp->x_diff[2];
+ dx2_max = max(dx2_max, dx2);
+
+ /* Maximal smoothing length */
+ cell_h_max = max(cell_h_max, bp->h);
+
+ /* Get ready for a density calculation */
+ if (bpart_is_active(bp, e)) {
+ black_holes_init_bpart(bp);
+ }
+ }
+
+ /* Now, get the maximal particle motion from its square */
+ dx_max = sqrtf(dx2_max);
+
+ /* Store the values */
+ c->black_holes.h_max = cell_h_max;
+ c->black_holes.dx_max_part = dx_max;
+
+ /* Update the time of the last drift */
+ c->black_holes.ti_old_part = ti_current;
+ }
+
+ /* Clear the drift flags. */
+ c->black_holes.do_drift = 0;
+ c->black_holes.do_sub_drift = 0;
+}
+
/**
* @brief Recursively drifts all multipoles in a cell hierarchy.
*
@@ -4790,6 +5089,34 @@ void cell_remove_spart(const struct engine *e, struct cell *c,
sp->gpart = NULL;
}
+/**
+ * @brief "Remove" a black hole particle from the calculation.
+ *
+ * The particle is inhibited and will officially be removed at the next rebuild.
+ *
+ * @param e The #engine running on this node.
+ * @param c The #cell from which to remove the particle.
+ * @param sp The #bpart to remove.
+ */
+void cell_remove_bpart(const struct engine *e, struct cell *c,
+ struct bpart *bp) {
+
+ /* Quick cross-check */
+ if (c->nodeID != e->nodeID)
+ error("Can't remove a particle in a foreign cell.");
+
+ /* Mark the particle as inhibited and stand-alone */
+ bp->time_bin = time_bin_inhibited;
+ if (bp->gpart) {
+ bp->gpart->time_bin = time_bin_inhibited;
+ bp->gpart->id_or_neg_offset = bp->id;
+ bp->gpart->type = swift_type_dark_matter;
+ }
+
+ /* Un-link the bpart */
+ bp->gpart = NULL;
+}
+
/**
* @brief "Remove" a gas particle from the calculation and convert its gpart
* friend to a dark matter particle.
diff --git a/src/cell.h b/src/cell.h
index b9b77a002060e73c2fd384bf82471a5a73a61874..6795e1ffa19a3f5cc651da5ef375fb7151c7c731 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -208,6 +208,18 @@ struct pcell_step_stars {
float dx_max_part;
};
+struct pcell_step_black_holes {
+
+ /*! Minimal integer end-of-timestep in this cell (black_holes) */
+ integertime_t ti_end_min;
+
+ /*! Maximal integer end-of-timestep in this cell (black_holes) */
+ integertime_t ti_end_max;
+
+ /*! Maximal distance any #part has travelled since last rebuild */
+ float dx_max_part;
+};
+
/**
* @brief Cell within the tree structure.
*
@@ -548,7 +560,7 @@ struct cell {
/*! Do any of this cell's sub-cells need to be sorted? */
char do_sub_sort;
- /*! Maximum end of (integer) time step in this cell for gravity tasks. */
+ /*! Maximum end of (integer) time step in this cell for star tasks. */
integertime_t ti_end_min;
/*! Maximum end of (integer) time step in this cell for star tasks. */
@@ -580,6 +592,64 @@ struct cell {
} stars;
+ /*! Black hole variables */
+ struct {
+
+ /*! Pointer to the #bpart data. */
+ struct bpart *parts;
+
+ /*! The drift task for bparts */
+ struct task *drift;
+
+ /*! Max smoothing length in this cell. */
+ double h_max;
+
+ /*! Last (integer) time the cell's bpart were drifted forward in time. */
+ integertime_t ti_old_part;
+
+ /*! Spin lock for various uses (#bpart case). */
+ swift_lock_type lock;
+
+ /*! Nr of #bpart in this cell. */
+ int count;
+
+ /*! Nr of #bpart this cell can hold after addition of new #bpart. */
+ int count_total;
+
+ /*! Values of h_max before the drifts, used for sub-cell tasks. */
+ float h_max_old;
+
+ /*! Maximum part movement in this cell since last construction. */
+ float dx_max_part;
+
+ /*! Maximum end of (integer) time step in this cell for black tasks. */
+ integertime_t ti_end_min;
+
+ /*! Maximum end of (integer) time step in this cell for black hole tasks. */
+ integertime_t ti_end_max;
+
+ /*! Maximum beginning of (integer) time step in this cell for black hole
+ * tasks.
+ */
+ integertime_t ti_beg_max;
+
+ /*! Number of #bpart updated in this cell. */
+ int updated;
+
+ /*! Number of #bpart inhibited in this cell. */
+ int inhibited;
+
+ /*! Is the #bpart data of this cell being used in a sub-cell? */
+ int hold;
+
+ /*! Does this cell need to be drifted (black holes)? */
+ char do_drift;
+
+ /*! Do any of this cell's sub-cells need to be drifted (black holes)? */
+ char do_sub_drift;
+
+ } black_holes;
+
#ifdef WITH_MPI
/*! MPI variables */
struct {
@@ -723,7 +793,8 @@ struct cell {
/* Function prototypes. */
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
- struct cell_buff *buff, struct cell_buff *sbuff,
+ ptrdiff_t bparts_offset, struct cell_buff *buff,
+ struct cell_buff *sbuff, struct cell_buff *bbuff,
struct cell_buff *gbuff);
void cell_sanitize(struct cell *c, int treated);
int cell_locktree(struct cell *c);
@@ -745,6 +816,8 @@ int cell_pack_end_step_grav(struct cell *c, struct pcell_step_grav *pcell);
int cell_unpack_end_step_grav(struct cell *c, struct pcell_step_grav *pcell);
int cell_pack_end_step_stars(struct cell *c, struct pcell_step_stars *pcell);
int cell_unpack_end_step_stars(struct cell *c, struct pcell_step_stars *pcell);
+int cell_pack_end_step_black_holes(struct cell *c, struct pcell_step_black_holes *pcell);
+int cell_unpack_end_step_black_holes(struct cell *c, struct pcell_step_black_holes *pcell);
int cell_pack_multipoles(struct cell *c, struct gravity_tensors *m);
int cell_unpack_multipoles(struct cell *c, struct gravity_tensors *m);
int cell_getsize(struct cell *c);
@@ -771,6 +844,7 @@ 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);
void cell_drift_spart(struct cell *c, const struct engine *e, int force);
+void cell_drift_bpart(struct cell *c, const struct engine *e, int force);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
@@ -803,6 +877,8 @@ void cell_remove_gpart(const struct engine *e, struct cell *c,
struct gpart *gp);
void cell_remove_spart(const struct engine *e, struct cell *c,
struct spart *sp);
+void cell_remove_bpart(const struct engine *e, struct cell *c,
+ struct bpart *bp);
struct spart *cell_add_spart(struct engine *e, struct cell *c);
struct gpart *cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
struct part *p, struct xpart *xp);
diff --git a/src/common_io.c b/src/common_io.c
index 2b8a55f0eba047a6076a3023569d05ddf49b376b..1b56fb43b4c774f87cf59c296b31e66a1d0806d3 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -950,6 +950,66 @@ void io_convert_spart_l_mapper(void* restrict temp, int N,
props.convert_spart_l(e, sparts + delta + i, &temp_l[i * dim]);
}
+/**
+ * @brief Mapper function to copy #bpart into a buffer of floats using a
+ * conversion function.
+ */
+void io_convert_bpart_f_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct bpart* restrict bparts = props.bparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ float* restrict temp_f = (float*)temp;
+ const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_bpart_f(e, bparts + delta + i, &temp_f[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #bpart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_bpart_d_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct bpart* restrict bparts = props.bparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ double* restrict temp_d = (double*)temp;
+ const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_bpart_d(e, bparts + delta + i, &temp_d[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #bpart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_bpart_l_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct bpart* restrict bparts = props.bparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ long long* restrict temp_l = (long long*)temp;
+ const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_bpart_l(e, bparts + delta + i, &temp_l[i * dim]);
+}
+
/**
* @brief Copy the particle data into a temporary buffer ready for i/o.
*
@@ -1091,6 +1151,42 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
io_convert_spart_l_mapper, temp_l, N, copySize, 0,
(void*)&props);
+ } else if (props.convert_bpart_f != NULL) {
+
+ /* Prepare some parameters */
+ float* temp_f = (float*)temp;
+ props.start_temp_f = (float*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_bpart_f_mapper, temp_f, N, copySize, 0,
+ (void*)&props);
+
+ } else if (props.convert_bpart_d != NULL) {
+
+ /* Prepare some parameters */
+ double* temp_d = (double*)temp;
+ props.start_temp_d = (double*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_bpart_d_mapper, temp_d, N, copySize, 0,
+ (void*)&props);
+
+ } else if (props.convert_bpart_l != NULL) {
+
+ /* Prepare some parameters */
+ long long* temp_l = (long long*)temp;
+ props.start_temp_l = (long long*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_bpart_l_mapper, temp_l, N, copySize, 0,
+ (void*)&props);
+
} else {
error("Missing conversion function");
}
@@ -1155,9 +1251,11 @@ struct duplication_data {
struct part* parts;
struct gpart* gparts;
struct spart* sparts;
+ struct bpart* bparts;
int Ndm;
int Ngas;
int Nstars;
+ int Nblackholes;
};
void io_duplicate_hydro_gparts_mapper(void* restrict data, int Ngas,
@@ -1216,7 +1314,7 @@ void io_duplicate_hydro_gparts(struct threadpool* tp, struct part* const parts,
sizeof(struct part), 0, &data);
}
-void io_duplicate_hydro_sparts_mapper(void* restrict data, int Nstars,
+void io_duplicate_stars_gparts_mapper(void* restrict data, int Nstars,
void* restrict extra_data) {
struct duplication_data* temp = (struct duplication_data*)extra_data;
@@ -1270,10 +1368,70 @@ void io_duplicate_stars_gparts(struct threadpool* tp,
data.sparts = sparts;
data.Ndm = Ndm;
- threadpool_map(tp, io_duplicate_hydro_sparts_mapper, sparts, Nstars,
+ threadpool_map(tp, io_duplicate_stars_gparts_mapper, sparts, Nstars,
sizeof(struct spart), 0, &data);
}
+void io_duplicate_black_holes_gparts_mapper(void* restrict data,
+ int Nblackholes,
+ void* restrict extra_data) {
+
+ struct duplication_data* temp = (struct duplication_data*)extra_data;
+ const int Ndm = temp->Ndm;
+ struct bpart* bparts = (struct bpart*)data;
+ const ptrdiff_t offset = bparts - temp->bparts;
+ struct gpart* gparts = temp->gparts + offset;
+
+ for (int i = 0; i < Nblackholes; ++i) {
+
+ /* Duplicate the crucial information */
+ gparts[i + Ndm].x[0] = bparts[i].x[0];
+ gparts[i + Ndm].x[1] = bparts[i].x[1];
+ gparts[i + Ndm].x[2] = bparts[i].x[2];
+
+ gparts[i + Ndm].v_full[0] = bparts[i].v[0];
+ gparts[i + Ndm].v_full[1] = bparts[i].v[1];
+ gparts[i + Ndm].v_full[2] = bparts[i].v[2];
+
+ gparts[i + Ndm].mass = bparts[i].mass;
+
+ /* Set gpart type */
+ gparts[i + Ndm].type = swift_type_black_hole;
+
+ /* Link the particles */
+ gparts[i + Ndm].id_or_neg_offset = -(long long)(offset + i);
+ bparts[i].gpart = &gparts[i + Ndm];
+ }
+}
+
+/**
+ * @brief Copy every #bpart into the corresponding #gpart and link them.
+ *
+ * This function assumes that the DM particles, gas particles and star particles
+ * are all at the start of the gparts array and adds the black hole particles
+ * afterwards
+ *
+ * @param tp The current #threadpool.
+ * @param bparts The array of #bpart freshly read in.
+ * @param gparts The array of #gpart freshly read in with all the DM, gas
+ * and star particles at the start.
+ * @param Nblackholes The number of blackholes particles read in.
+ * @param Ndm The number of DM, gas and star particles read in.
+ */
+void io_duplicate_black_holes_gparts(struct threadpool* tp,
+ struct bpart* const bparts,
+ struct gpart* const gparts,
+ size_t Nblackholes, size_t Ndm) {
+
+ struct duplication_data data;
+ data.gparts = gparts;
+ data.bparts = bparts;
+ data.Ndm = Ndm;
+
+ threadpool_map(tp, io_duplicate_black_holes_gparts_mapper, bparts,
+ Nblackholes, sizeof(struct bpart), 0, &data);
+}
+
/**
* @brief Copy every non-inhibited #part into the parts_written array.
*
@@ -1348,6 +1506,40 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
count, Nsparts_written);
}
+/**
+ * @brief Copy every non-inhibited #bpart into the bparts_written array.
+ *
+ * @param bparts The array of #bpart containing all particles.
+ * @param bparts_written The array of #bpart to fill with particles we want to
+ * write.
+ * @param Nbparts The total number of #part.
+ * @param Nbparts_written The total number of #part to write.
+ */
+void io_collect_bparts_to_write(const struct bpart* restrict bparts,
+ struct bpart* restrict bparts_written,
+ const size_t Nbparts,
+ const size_t Nbparts_written) {
+
+ size_t count = 0;
+
+ /* Loop over all parts */
+ for (size_t i = 0; i < Nbparts; ++i) {
+
+ /* And collect the ones that have not been removed */
+ if (bparts[i].time_bin != time_bin_inhibited &&
+ bparts[i].time_bin != time_bin_not_created) {
+
+ bparts_written[count] = bparts[i];
+ count++;
+ }
+ }
+
+ /* Check that everything is fine */
+ if (count != Nbparts_written)
+ error("Collected the wrong number of s-particles (%zu vs. %zu expected)",
+ count, Nbparts_written);
+}
+
/**
* @brief Copy every non-inhibited DM #gpart into the gparts_written array.
*
diff --git a/src/common_io.h b/src/common_io.h
index eb1ee0a804f324d897842fb2a0ca33fc07e769d6..c93414d1ffca8f7ead7a1ff29d387967f82a9cb2 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -38,6 +38,7 @@ struct part;
struct gpart;
struct velociraptor_gpart_data;
struct spart;
+struct bpart;
struct xpart;
struct io_props;
struct engine;
@@ -113,6 +114,10 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
struct spart* restrict sparts_written,
const size_t Nsparts,
const size_t Nsparts_written);
+void io_collect_bparts_to_write(const struct bpart* restrict bparts,
+ struct bpart* restrict bparts_written,
+ const size_t Nbparts,
+ const size_t Nbparts_written);
void io_collect_gparts_to_write(const struct gpart* restrict gparts,
const struct velociraptor_gpart_data* vr_data,
struct gpart* restrict gparts_written,
@@ -128,6 +133,10 @@ void io_duplicate_stars_gparts(struct threadpool* tp,
struct spart* const sparts,
struct gpart* const gparts, size_t Nstars,
size_t Ndm);
+void io_duplicate_black_holes_gparts(struct threadpool* tp,
+ struct bpart* const bparts,
+ struct gpart* const gparts, size_t Nstars,
+ size_t Ndm);
void io_check_output_fields(const struct swift_params* params,
const long long N_total[3]);
diff --git a/src/drift.h b/src/drift.h
index 7e874fe0ceabe5b091cc7c5bb53adbef2c9a3efd..faa8743a86e47ab53d9880e086d5acf454ca7c38 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -23,6 +23,7 @@
#include "../config.h"
/* Local headers. */
+#include "black_holes.h"
#include "const.h"
#include "debug.h"
#include "dimension.h"
@@ -150,4 +151,42 @@ __attribute__((always_inline)) INLINE static void drift_spart(
}
}
+/**
+ * @brief Perform the 'drift' operation on a #bpart
+ *
+ * @param bp The #bpart to drift.
+ * @param dt_drift The drift time-step.
+ * @param ti_old Integer start of time-step (for debugging checks).
+ * @param ti_current Integer end of time-step (for debugging checks).
+ */
+__attribute__((always_inline)) INLINE static void drift_bpart(
+ struct bpart *restrict bp, double dt_drift, integertime_t ti_old,
+ integertime_t ti_current) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (bp->ti_drift != ti_old)
+ error(
+ "s-particle has not been drifted to the current time "
+ "bp->ti_drift=%lld, "
+ "c->ti_old=%lld, ti_current=%lld",
+ bp->ti_drift, ti_old, ti_current);
+
+ bp->ti_drift = ti_current;
+#endif
+
+ /* Drift... */
+ bp->x[0] += bp->v[0] * dt_drift;
+ bp->x[1] += bp->v[1] * dt_drift;
+ bp->x[2] += bp->v[2] * dt_drift;
+
+ /* Predict the values of the extra fields */
+ black_holes_predict_extra(bp, dt_drift);
+
+ /* Compute offsets since last cell construction */
+ for (int k = 0; k < 3; k++) {
+ const float dx = bp->v[k] * dt_drift;
+ bp->x_diff[k] -= dx;
+ }
+}
+
#endif /* SWIFT_DRIFT_H */
diff --git a/src/engine.c b/src/engine.c
index e0747dde19f44f20b370b2d1d806f4a7e757bb3c..eb064ec0dff2fde05016faee1628b0fc58a1d7a7 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -117,6 +117,7 @@ const char *engine_policy_names[] = {"none",
"structure finding",
"star formation",
"feedback",
+ "black holes",
"time-step limiter"};
/** The rank of the engine as a global variable (for messages). */
@@ -365,6 +366,14 @@ static void ENGINE_REDISTRIBUTE_DEST_MAPPER(spart);
*/
static void ENGINE_REDISTRIBUTE_DEST_MAPPER(gpart);
+/**
+ * @brief Accumulate the counts of black holes particles per cell.
+ * Threadpool helper for accumulating the counts of particles per cell.
+ *
+ * bpart version.
+ */
+static void ENGINE_REDISTRIBUTE_DEST_MAPPER(bpart);
+
#endif /* redist_mapper_data */
#ifdef WITH_MPI /* savelink_mapper_data */
@@ -435,17 +444,29 @@ static void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(spart, 1);
static void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(spart, 0);
#endif
+/**
+ * @brief Save position of bpart-gpart links.
+ * Threadpool helper for accumulating the counts of particles per cell.
+ */
+#ifdef SWIFT_DEBUG_CHECKS
+static void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 1);
+#else
+static void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 0);
+#endif
+
#endif /* savelink_mapper_data */
#ifdef WITH_MPI /* relink_mapper_data */
-/* Support for relinking parts, gparts and sparts after moving between nodes. */
+/* Support for relinking parts, gparts, sparts and bparts after moving between
+ * nodes. */
struct relink_mapper_data {
int nodeID;
int nr_nodes;
int *counts;
int *s_counts;
int *g_counts;
+ int *b_counts;
struct space *s;
};
@@ -468,6 +489,7 @@ static void engine_redistribute_relink_mapper(void *map_data, int num_elements,
int *counts = mydata->counts;
int *g_counts = mydata->g_counts;
int *s_counts = mydata->s_counts;
+ int *b_counts = mydata->b_counts;
struct space *s = mydata->s;
for (int i = 0; i < num_elements; i++) {
@@ -478,11 +500,13 @@ static void engine_redistribute_relink_mapper(void *map_data, int num_elements,
size_t offset_parts = 0;
size_t offset_gparts = 0;
size_t offset_sparts = 0;
+ size_t offset_bparts = 0;
for (int n = 0; n < node; n++) {
int ind_recv = n * nr_nodes + nodeID;
offset_parts += counts[ind_recv];
offset_gparts += g_counts[ind_recv];
offset_sparts += s_counts[ind_recv];
+ offset_bparts += b_counts[ind_recv];
}
/* Number of gparts sent from this node. */
@@ -513,6 +537,17 @@ static void engine_redistribute_relink_mapper(void *map_data, int num_elements,
s->gparts[k].id_or_neg_offset = -partner_index;
s->sparts[partner_index].gpart = &s->gparts[k];
}
+
+ /* Does this gpart have a black hole partner ? */
+ else if (s->gparts[k].type == swift_type_black_hole) {
+
+ const ptrdiff_t partner_index =
+ offset_bparts - s->gparts[k].id_or_neg_offset;
+
+ /* Re-link */
+ s->gparts[k].id_or_neg_offset = -partner_index;
+ s->sparts[partner_index].gpart = &s->gparts[k];
+ }
}
}
}
@@ -548,11 +583,13 @@ void engine_redistribute(struct engine *e) {
struct part *parts = s->parts;
struct gpart *gparts = s->gparts;
struct spart *sparts = s->sparts;
+ struct bpart *bparts = s->bparts;
ticks tic = getticks();
size_t nr_parts = s->nr_parts;
size_t nr_gparts = s->nr_gparts;
size_t nr_sparts = s->nr_sparts;
+ size_t nr_bparts = s->nr_bparts;
/* Start by moving inhibited particles to the end of the arrays */
for (size_t k = 0; k < nr_parts; /* void */) {
@@ -599,6 +636,27 @@ void engine_redistribute(struct engine *e) {
}
}
+ /* Now move inhibited black hole particles to the end of the arrays */
+ for (size_t k = 0; k < nr_bparts; /* void */) {
+ if (bparts[k].time_bin == time_bin_inhibited ||
+ bparts[k].time_bin == time_bin_not_created) {
+ nr_bparts -= 1;
+
+ /* Swap the particle */
+ memswap(&s->bparts[k], &s->bparts[nr_bparts], sizeof(struct bpart));
+
+ /* Swap the link with the gpart */
+ if (s->bparts[k].gpart != NULL) {
+ s->bparts[k].gpart->id_or_neg_offset = -k;
+ }
+ if (s->bparts[nr_bparts].gpart != NULL) {
+ s->bparts[nr_bparts].gpart->id_or_neg_offset = -nr_bparts;
+ }
+ } else {
+ k++;
+ }
+ }
+
/* Finally do the same with the gravity particles */
for (size_t k = 0; k < nr_gparts; /* void */) {
if (gparts[k].time_bin == time_bin_inhibited ||
@@ -613,13 +671,19 @@ void engine_redistribute(struct engine *e) {
s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_stars) {
s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+ } else if (s->gparts[k].type == swift_type_black_hole) {
+ s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
+
if (s->gparts[nr_gparts].type == swift_type_gas) {
s->parts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
} else if (s->gparts[nr_gparts].type == swift_type_stars) {
s->sparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
+ } else if (s->gparts[nr_gparts].type == swift_type_black_hole) {
+ s->bparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
+ &s->gparts[nr_gparts];
}
} else {
k++;
@@ -770,6 +834,70 @@ void engine_redistribute(struct engine *e) {
}
swift_free("s_dest", s_dest);
+ /* Get destination of each b-particle */
+ int *b_counts;
+ if ((b_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
+ error("Failed to allocate b_counts temporary buffer.");
+
+ int *b_dest;
+ if ((b_dest = (int *)swift_malloc("b_dest", sizeof(int) * nr_bparts)) == NULL)
+ error("Failed to allocate b_dest temporary buffer.");
+
+ redist_data.counts = b_counts;
+ redist_data.dest = b_dest;
+ redist_data.base = (void *)bparts;
+
+ threadpool_map(&e->threadpool, engine_redistribute_dest_mapper_bpart, bparts,
+ nr_bparts, sizeof(struct bpart), 0, &redist_data);
+
+ /* Sort the particles according to their cell index. */
+ if (nr_bparts > 0)
+ space_bparts_sort(s->bparts, b_dest, &b_counts[nodeID * nr_nodes], nr_nodes,
+ 0);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Verify that the bpart have been sorted correctly. */
+ for (size_t k = 0; k < nr_bparts; k++) {
+ const struct bpart *bp = &s->bparts[k];
+
+ if (bp->time_bin == time_bin_inhibited)
+ error("Inhibited particle found after sorting!");
+
+ if (bp->time_bin == time_bin_not_created)
+ error("Inhibited particle found after sorting!");
+
+ /* New cell index */
+ const int new_cid =
+ cell_getid(s->cdim, bp->x[0] * s->iwidth[0], bp->x[1] * s->iwidth[1],
+ bp->x[2] * s->iwidth[2]);
+
+ /* New cell of this bpart */
+ const struct cell *c = &s->cells_top[new_cid];
+ const int new_node = c->nodeID;
+
+ if (s_dest[k] != new_node)
+ error("bpart's new node index not matching sorted index.");
+
+ if (bp->x[0] < c->loc[0] || bp->x[0] > c->loc[0] + c->width[0] ||
+ bp->x[1] < c->loc[1] || bp->x[1] > c->loc[1] + c->width[1] ||
+ bp->x[2] < c->loc[2] || bp->x[2] > c->loc[2] + c->width[2])
+ error("bpart not sorted into the right top-level cell!");
+ }
+#endif
+
+ /* We need to re-link the gpart partners of bparts. */
+ if (nr_bparts > 0) {
+
+ struct savelink_mapper_data savelink_data;
+ savelink_data.nr_nodes = nr_nodes;
+ savelink_data.counts = s_counts;
+ savelink_data.parts = (void *)bparts;
+ savelink_data.nodeID = nodeID;
+ threadpool_map(&e->threadpool, engine_redistribute_savelink_mapper_bpart,
+ nodes, nr_nodes, sizeof(int), 0, &savelink_data);
+ }
+ swift_free("b_dest", b_dest);
+
/* Get destination of each g-particle */
int *g_counts;
if ((g_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
@@ -788,7 +916,7 @@ void engine_redistribute(struct engine *e) {
/* Sort the gparticles according to their cell index. */
if (nr_gparts > 0)
- space_gparts_sort(s->gparts, s->parts, s->sparts, g_dest,
+ space_gparts_sort(s->gparts, s->parts, s->sparts, s->bparts, g_dest,
&g_counts[nodeID * nr_nodes], nr_nodes);
#ifdef SWIFT_DEBUG_CHECKS
@@ -829,30 +957,37 @@ void engine_redistribute(struct engine *e) {
MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to allreduce particle transfer counts.");
- /* Get all the s_counts from all the nodes. */
+ /* Get all the g_counts from all the nodes. */
if (MPI_Allreduce(MPI_IN_PLACE, g_counts, nr_nodes * nr_nodes, MPI_INT,
MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to allreduce gparticle transfer counts.");
- /* Get all the g_counts from all the nodes. */
+ /* Get all the s_counts from all the nodes. */
if (MPI_Allreduce(MPI_IN_PLACE, s_counts, nr_nodes * nr_nodes, MPI_INT,
MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to allreduce sparticle transfer counts.");
+ /* Get all the b_counts from all the nodes. */
+ if (MPI_Allreduce(MPI_IN_PLACE, b_counts, nr_nodes * nr_nodes, MPI_INT,
+ MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
+ error("Failed to allreduce bparticle transfer counts.");
+
/* Report how many particles will be moved. */
if (e->verbose) {
if (e->nodeID == 0) {
- size_t total = 0, g_total = 0, s_total = 0;
- size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0;
+ size_t total = 0, g_total = 0, s_total = 0, b_total = 0;
+ size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0;
for (int p = 0, r = 0; p < nr_nodes; p++) {
for (int n = 0; n < nr_nodes; n++) {
total += counts[r];
g_total += g_counts[r];
s_total += s_counts[r];
+ b_total += b_counts[r];
if (p == n) {
unmoved += counts[r];
g_unmoved += g_counts[r];
s_unmoved += s_counts[r];
+ b_unmoved += b_counts[r];
}
r++;
}
@@ -868,19 +1003,26 @@ void engine_redistribute(struct engine *e) {
message("%ld of %ld (%.2f%%) of s-particles moved", s_total - s_unmoved,
s_total,
100.0 * (double)(s_total - s_unmoved) / (double)s_total);
+ if (b_total > 0)
+ message("%ld of %ld (%.2f%%) of b-particles moved", b_total - b_unmoved,
+ b_total,
+ 100.0 * (double)(b_total - b_unmoved) / (double)b_total);
}
}
- /* Now each node knows how many parts, sparts and gparts will be transferred
- * to every other node.
- * Get the new numbers of particles for this node. */
- size_t nr_parts_new = 0, nr_gparts_new = 0, nr_sparts_new = 0;
+ /* Now each node knows how many parts, sparts, bparts, and gparts will be
+ * transferred to every other node. Get the new numbers of particles for this
+ * node. */
+ size_t nr_parts_new = 0, nr_gparts_new = 0, nr_sparts_new = 0,
+ nr_bparts_new = 0;
for (int k = 0; k < nr_nodes; k++)
nr_parts_new += counts[k * nr_nodes + nodeID];
for (int k = 0; k < nr_nodes; k++)
nr_gparts_new += g_counts[k * nr_nodes + nodeID];
for (int k = 0; k < nr_nodes; k++)
nr_sparts_new += s_counts[k * nr_nodes + nodeID];
+ for (int k = 0; k < nr_nodes; k++)
+ nr_bparts_new += b_counts[k * nr_nodes + nodeID];
/* Now exchange the particles, type by type to keep the memory required
* under control. */
@@ -919,6 +1061,15 @@ void engine_redistribute(struct engine *e) {
s->nr_sparts = nr_sparts_new;
s->size_sparts = engine_redistribute_alloc_margin * nr_sparts_new;
+ /* Black holes particles. */
+ new_parts = engine_do_redistribute(
+ "bparts", b_counts, (char *)s->bparts, nr_bparts_new,
+ sizeof(struct bpart), bpart_align, bpart_mpi_type, nr_nodes, nodeID);
+ swift_free("bparts", s->bparts);
+ s->bparts = (struct bpart *)new_parts;
+ s->nr_bparts = nr_bparts_new;
+ s->size_bparts = engine_redistribute_alloc_margin * nr_bparts_new;
+
/* All particles have now arrived. Time for some final operations on the
stuff we just received */
@@ -930,6 +1081,7 @@ void engine_redistribute(struct engine *e) {
relink_data.counts = counts;
relink_data.g_counts = g_counts;
relink_data.s_counts = s_counts;
+ relink_data.b_counts = b_counts;
relink_data.nodeID = nodeID;
relink_data.nr_nodes = nr_nodes;
@@ -941,6 +1093,7 @@ void engine_redistribute(struct engine *e) {
free(counts);
free(g_counts);
free(s_counts);
+ free(b_counts);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that all parts are in the right place. */
@@ -968,10 +1121,18 @@ void engine_redistribute(struct engine *e) {
error("Received s-particle (%zu) that does not belong here (nodeID=%i).",
k, cells[cid].nodeID);
}
+ for (size_t k = 0; k < nr_bparts_new; k++) {
+ const int cid = cell_getid(s->cdim, s->bparts[k].x[0] * s->iwidth[0],
+ s->bparts[k].x[1] * s->iwidth[1],
+ s->bparts[k].x[2] * s->iwidth[2]);
+ if (cells[cid].nodeID != nodeID)
+ error("Received p-particle (%zu) that does not belong here (nodeID=%i).",
+ k, cells[cid].nodeID);
+ }
/* Verify that the links are correct */
- part_verify_links(s->parts, s->gparts, s->sparts, nr_parts_new, nr_gparts_new,
- nr_sparts_new, e->verbose);
+ part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, nr_parts_new,
+ nr_gparts_new, nr_sparts_new, nr_bparts_new, e->verbose);
#endif
@@ -980,14 +1141,18 @@ void engine_redistribute(struct engine *e) {
int my_cells = 0;
for (int k = 0; k < nr_cells; k++)
if (cells[k].nodeID == nodeID) my_cells += 1;
- message("node %i now has %zu parts, %zu sparts and %zu gparts in %i cells.",
- nodeID, nr_parts_new, nr_sparts_new, nr_gparts_new, my_cells);
+ message(
+ "node %i now has %zu parts, %zu sparts, %zu bparts and %zu gparts in "
+ "%i cells.",
+ nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_gparts_new,
+ my_cells);
}
/* Flag that we do not have any extra particles any more */
s->nr_extra_parts = 0;
s->nr_extra_gparts = 0;
s->nr_extra_sparts = 0;
+ s->nr_extra_bparts = 0;
/* Flag that a redistribute has taken place */
e->step_props |= engine_step_prop_redistribute;
@@ -2186,22 +2351,25 @@ void engine_rebuild(struct engine *e, int repartitioned,
const ticks tic2 = getticks();
/* Update the global counters of particles */
- long long num_particles[3] = {
+ long long num_particles[4] = {
(long long)(e->s->nr_parts - e->s->nr_extra_parts),
(long long)(e->s->nr_gparts - e->s->nr_extra_gparts),
- (long long)(e->s->nr_sparts - e->s->nr_extra_sparts)};
+ (long long)(e->s->nr_sparts - e->s->nr_extra_sparts),
+ (long long)(e->s->nr_bparts - e->s->nr_extra_bparts)};
#ifdef WITH_MPI
- MPI_Allreduce(MPI_IN_PLACE, num_particles, 3, MPI_LONG_LONG, MPI_SUM,
+ MPI_Allreduce(MPI_IN_PLACE, num_particles, 4, MPI_LONG_LONG, MPI_SUM,
MPI_COMM_WORLD);
#endif
e->total_nr_parts = num_particles[0];
e->total_nr_gparts = num_particles[1];
e->total_nr_sparts = num_particles[2];
+ e->total_nr_bparts = num_particles[3];
/* Flag that there are no inhibited particles */
e->nr_inhibited_parts = 0;
e->nr_inhibited_gparts = 0;
e->nr_inhibited_sparts = 0;
+ e->nr_inhibited_bparts = 0;
if (e->verbose)
message("updating particle counts took %.3f %s.",
@@ -2216,8 +2384,9 @@ void engine_rebuild(struct engine *e, int repartitioned,
engine_recompute_displacement_constraint(e);
#ifdef SWIFT_DEBUG_CHECKS
- part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->nr_parts,
- e->s->nr_gparts, e->s->nr_sparts, e->verbose);
+ part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
+ e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
+ e->s->nr_bparts, e->verbose);
#endif
/* Initial cleaning up session ? */
@@ -2276,6 +2445,7 @@ void engine_rebuild(struct engine *e, int repartitioned,
e->updates_since_rebuild = 0;
e->g_updates_since_rebuild = 0;
e->s_updates_since_rebuild = 0;
+ e->b_updates_since_rebuild = 0;
/* Flag that a rebuild has taken place */
e->step_props |= engine_step_prop_rebuild;
@@ -2975,6 +3145,7 @@ void engine_first_init_particles(struct engine *e) {
space_first_init_parts(e->s, e->verbose);
space_first_init_gparts(e->s, e->verbose);
space_first_init_sparts(e->s, e->verbose);
+ space_first_init_bparts(e->s, e->verbose);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -3027,6 +3198,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
space_init_parts(s, e->verbose);
space_init_gparts(s, e->verbose);
space_init_sparts(s, e->verbose);
+ space_init_bparts(s, e->verbose);
/* Update the cooling function */
if ((e->policy & engine_policy_cooling) ||
@@ -3090,6 +3262,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
space_init_parts(e->s, e->verbose);
space_init_gparts(e->s, e->verbose);
space_init_sparts(e->s, e->verbose);
+ space_init_bparts(e->s, e->verbose);
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
@@ -3215,8 +3388,9 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
#ifdef SWIFT_DEBUG_CHECKS
space_check_timesteps(e->s);
- part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->nr_parts,
- e->s->nr_gparts, e->s->nr_sparts, e->verbose);
+ part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
+ e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
+ e->s->nr_bparts, e->verbose);
#endif
/* Ready to go */
@@ -4123,8 +4297,8 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the links are correct */
- part_verify_links(s->parts, s->gparts, s->sparts, s->nr_parts, s->nr_gparts,
- s->nr_sparts, e->verbose);
+ part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
+ s->nr_gparts, s->nr_sparts, s->nr_bparts, e->verbose);
#endif
if (e->verbose)
diff --git a/src/engine.h b/src/engine.h
index e0c52dafadeb2996600e3e3b01349ae718c985dc..c9557c189cae013e3c21c02b54ca65c8c9907ec8 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -75,9 +75,10 @@ enum engine_policy {
engine_policy_structure_finding = (1 << 16),
engine_policy_star_formation = (1 << 17),
engine_policy_feedback = (1 << 18),
- engine_policy_limiter = (1 << 19)
+ engine_policy_black_holes = (1 << 19),
+ engine_policy_limiter = (1 << 20)
};
-#define engine_maxpolicy 20
+#define engine_maxpolicy 21
extern const char *engine_policy_names[engine_maxpolicy + 1];
/**
@@ -213,18 +214,22 @@ struct engine {
integertime_t ti_beg_max;
/* Number of particles updated in the previous step */
- long long updates, g_updates, s_updates;
+ long long updates, g_updates, s_updates, b_updates;
/* Number of updates since the last rebuild */
long long updates_since_rebuild;
long long g_updates_since_rebuild;
long long s_updates_since_rebuild;
+ long long b_updates_since_rebuild;
/* Properties of the previous step */
int step_props;
/* Total numbers of particles in the system. */
- long long total_nr_parts, total_nr_gparts, total_nr_sparts;
+ long long total_nr_parts;
+ long long total_nr_gparts;
+ long long total_nr_sparts;
+ long long total_nr_bparts;
/* Total numbers of cells (top-level and sub-cells) in the system. */
long long total_nr_cells;
@@ -233,12 +238,17 @@ struct engine {
long long total_nr_tasks;
/* The total number of inhibited particles in the system. */
- long long nr_inhibited_parts, nr_inhibited_gparts, nr_inhibited_sparts;
+ long long nr_inhibited_parts;
+ long long nr_inhibited_gparts;
+ long long nr_inhibited_sparts;
+ long long nr_inhibited_bparts;
#ifdef SWIFT_DEBUG_CHECKS
/* Total number of particles removed from the system since the last rebuild */
- long long count_inhibited_parts, count_inhibited_gparts,
- count_inhibited_sparts;
+ long long count_inhibited_parts;
+ long long count_inhibited_gparts;
+ long long count_inhibited_sparts;
+ long long count_inhibited_bparts;
#endif
/* Total mass in the simulation */
diff --git a/src/engine_drift.c b/src/engine_drift.c
index 1b0711619d68da02753f307190ca3a0624feecce..fee9ed1f11de2d5e372bd70554bd053b8313ac80 100644
--- a/src/engine_drift.c
+++ b/src/engine_drift.c
@@ -172,6 +172,54 @@ void engine_do_drift_all_spart_mapper(void *map_data, int num_elements,
}
}
+/**
+ * @brief Mapper function to drift *all* the #bpart to the current time.
+ *
+ * @param map_data An array of #cell%s.
+ * @param num_elements Chunk size.
+ * @param extra_data Pointer to an #engine.
+ */
+void engine_do_drift_all_bpart_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+
+ const struct engine *e = (const struct engine *)extra_data;
+ const int restarting = e->restarting;
+ struct space *s = e->s;
+ struct cell *cells_top;
+ int *local_cells_top;
+
+ if (restarting) {
+
+ /* When restarting, we loop over all top-level cells */
+ cells_top = (struct cell *)map_data;
+ local_cells_top = NULL;
+
+ } else {
+
+ /* In any other case, we use the list of local cells with tasks */
+ cells_top = s->cells_top;
+ local_cells_top = (int *)map_data;
+ }
+
+ for (int ind = 0; ind < num_elements; ind++) {
+
+ struct cell *c;
+
+ /* When restarting, the list of local cells does not
+ yet exist. We use the raw list of top-level cells instead */
+ if (restarting)
+ c = &cells_top[ind];
+ else
+ c = &cells_top[local_cells_top[ind]];
+
+ if (c->nodeID == e->nodeID) {
+
+ /* Drift all the particles */
+ cell_drift_bpart(c, e, /* force the drift=*/1);
+ }
+ }
+}
+
/**
* @brief Mapper function to drift *all* the multipoles to the current time.
*
@@ -257,6 +305,11 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
/* default chunk */ 0, e);
}
+ if (e->s->nr_bparts > 0) {
+ threadpool_map(&e->threadpool, engine_do_drift_all_bpart_mapper,
+ e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
+ /* default chunk */ 0, e);
+ }
if (drift_mpoles && (e->policy & engine_policy_self_gravity)) {
threadpool_map(&e->threadpool, engine_do_drift_all_multipole_mapper,
e->s->local_cells_with_tasks_top,
@@ -274,6 +327,16 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
e->s->cells_top, e->s->nr_cells, sizeof(struct cell),
/* default chunk */ 0, e);
}
+ if (e->s->nr_sparts > 0) {
+ threadpool_map(&e->threadpool, engine_do_drift_all_spart_mapper,
+ e->s->cells_top, e->s->nr_cells, sizeof(struct cell),
+ /* default chunk */ 0, e);
+ }
+ if (e->s->nr_bparts > 0) {
+ threadpool_map(&e->threadpool, engine_do_drift_all_bpart_mapper,
+ e->s->cells_top, e->s->nr_cells, sizeof(struct cell),
+ /* default chunk */ 0, e);
+ }
if (e->s->nr_gparts > 0) {
threadpool_map(&e->threadpool, engine_do_drift_all_gpart_mapper,
e->s->cells_top, e->s->nr_cells, sizeof(struct cell),
@@ -294,8 +357,9 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
space_check_drift_point(
e->s, e->ti_current,
drift_mpoles && (e->policy & engine_policy_self_gravity));
- part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->nr_parts,
- e->s->nr_gparts, e->s->nr_sparts, e->verbose);
+ part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
+ e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
+ e->s->nr_bparts, e->verbose);
#endif
if (e->verbose)
diff --git a/src/io_properties.h b/src/io_properties.h
index c45edb2641e374e2cfaec6c3251aff7d18f361d6..8e5d31c9bdec0db3308006e391453ff5970c0127 100644
--- a/src/io_properties.h
+++ b/src/io_properties.h
@@ -60,6 +60,13 @@ typedef void (*conversion_func_spart_double)(const struct engine*,
typedef void (*conversion_func_spart_long_long)(const struct engine*,
const struct spart*,
long long*);
+typedef void (*conversion_func_bpart_float)(const struct engine*,
+ const struct bpart*, float*);
+typedef void (*conversion_func_bpart_double)(const struct engine*,
+ const struct bpart*, double*);
+typedef void (*conversion_func_bpart_long_long)(const struct engine*,
+ const struct bpart*,
+ long long*);
/**
* @brief The properties of a given dataset for i/o
@@ -101,6 +108,7 @@ struct io_props {
const struct xpart* xparts;
const struct gpart* gparts;
const struct spart* sparts;
+ const struct bpart* bparts;
/* Are we converting? */
int conversion;
@@ -119,6 +127,11 @@ struct io_props {
conversion_func_spart_float convert_spart_f;
conversion_func_spart_double convert_spart_d;
conversion_func_spart_long_long convert_spart_l;
+
+ /* Conversion function for bpart */
+ conversion_func_bpart_float convert_bpart_f;
+ conversion_func_bpart_double convert_bpart_d;
+ conversion_func_bpart_long_long convert_bpart_l;
};
/**
@@ -157,6 +170,7 @@ INLINE static struct io_props io_make_input_field_(
r.xparts = NULL;
r.gparts = NULL;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 0;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -167,6 +181,9 @@ INLINE static struct io_props io_make_input_field_(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -204,6 +221,7 @@ INLINE static struct io_props io_make_output_field_(
r.parts = NULL;
r.gparts = NULL;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 0;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -214,6 +232,9 @@ INLINE static struct io_props io_make_output_field_(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -258,6 +279,7 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
r.xparts = xparts;
r.gparts = NULL;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = functionPtr;
r.convert_part_d = NULL;
@@ -268,6 +290,9 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -304,6 +329,7 @@ INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
r.xparts = xparts;
r.gparts = NULL;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = functionPtr;
@@ -314,6 +340,9 @@ INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -350,6 +379,7 @@ INLINE static struct io_props io_make_output_field_convert_part_LONGLONG(
r.xparts = xparts;
r.gparts = NULL;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -360,6 +390,9 @@ INLINE static struct io_props io_make_output_field_convert_part_LONGLONG(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -402,6 +435,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
r.xparts = NULL;
r.gparts = gparts;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -412,6 +446,9 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -446,6 +483,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE(
r.xparts = NULL;
r.gparts = gparts;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -456,6 +494,9 @@ INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -490,6 +531,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_LONGLONG(
r.xparts = NULL;
r.gparts = gparts;
r.sparts = NULL;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -500,6 +542,9 @@ INLINE static struct io_props io_make_output_field_convert_gpart_LONGLONG(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -542,6 +587,7 @@ INLINE static struct io_props io_make_output_field_convert_spart_FLOAT(
r.xparts = NULL;
r.gparts = NULL;
r.sparts = sparts;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -552,6 +598,9 @@ INLINE static struct io_props io_make_output_field_convert_spart_FLOAT(
r.convert_spart_f = functionPtr;
r.convert_spart_d = NULL;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -586,6 +635,7 @@ INLINE static struct io_props io_make_output_field_convert_spart_DOUBLE(
r.xparts = NULL;
r.gparts = NULL;
r.sparts = sparts;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -596,6 +646,9 @@ INLINE static struct io_props io_make_output_field_convert_spart_DOUBLE(
r.convert_spart_f = NULL;
r.convert_spart_d = functionPtr;
r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
return r;
}
@@ -630,6 +683,7 @@ INLINE static struct io_props io_make_output_field_convert_spart_LONGLONG(
r.xparts = NULL;
r.gparts = NULL;
r.sparts = sparts;
+ r.bparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
r.convert_part_d = NULL;
@@ -640,6 +694,161 @@ INLINE static struct io_props io_make_output_field_convert_spart_LONGLONG(
r.convert_spart_f = NULL;
r.convert_spart_d = NULL;
r.convert_spart_l = functionPtr;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
+
+ return r;
+}
+
+/**
+ * @brief Constructs an #io_props (with conversion) from its parameters
+ */
+#define io_make_output_field_convert_bpart(name, type, dim, units, bpart, \
+ convert) \
+ io_make_output_field_convert_bpart_##type(name, type, dim, units, \
+ sizeof(bpart[0]), bpart, convert)
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param spartSize The size in byte of the particle
+ * @param sparts The particle array
+ * @param functionPtr The function used to convert a g-particle to a float
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_bpart_FLOAT(
+ const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+ enum unit_conversion_factor units, size_t bpartSize,
+ const struct bpart* bparts, conversion_func_bpart_float functionPtr) {
+
+ struct io_props r;
+ strcpy(r.name, name);
+ r.type = type;
+ r.dimension = dimension;
+ r.importance = UNUSED;
+ r.units = units;
+ r.field = NULL;
+ r.partSize = bpartSize;
+ r.parts = NULL;
+ r.xparts = NULL;
+ r.gparts = NULL;
+ r.sparts = NULL;
+ r.bparts = bparts;
+ r.conversion = 1;
+ r.convert_part_f = NULL;
+ r.convert_part_d = NULL;
+ r.convert_part_l = NULL;
+ r.convert_gpart_f = NULL;
+ r.convert_gpart_d = NULL;
+ r.convert_gpart_l = NULL;
+ r.convert_spart_f = NULL;
+ r.convert_spart_d = NULL;
+ r.convert_spart_l = NULL;
+ r.convert_bpart_f = functionPtr;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = NULL;
+
+ return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param spartSize The size in byte of the particle
+ * @param sparts The particle array
+ * @param functionPtr The function used to convert a s-particle to a double
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_bpart_DOUBLE(
+ const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+ enum unit_conversion_factor units, size_t bpartSize,
+ const struct bpart* bparts, conversion_func_bpart_double functionPtr) {
+
+ struct io_props r;
+ strcpy(r.name, name);
+ r.type = type;
+ r.dimension = dimension;
+ r.importance = UNUSED;
+ r.units = units;
+ r.field = NULL;
+ r.partSize = bpartSize;
+ r.parts = NULL;
+ r.xparts = NULL;
+ r.gparts = NULL;
+ r.sparts = NULL;
+ r.bparts = bparts;
+ r.conversion = 1;
+ r.convert_part_f = NULL;
+ r.convert_part_d = NULL;
+ r.convert_part_l = NULL;
+ r.convert_gpart_f = NULL;
+ r.convert_gpart_d = NULL;
+ r.convert_gpart_l = NULL;
+ r.convert_spart_f = NULL;
+ r.convert_spart_d = NULL;
+ r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = functionPtr;
+ r.convert_bpart_l = NULL;
+
+ return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param spartSize The size in byte of the particle
+ * @param sparts The particle array
+ * @param functionPtr The function used to convert a s-particle to a double
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_bpart_LONGLONG(
+ const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+ enum unit_conversion_factor units, size_t bpartSize,
+ const struct bpart* bparts, conversion_func_bpart_long_long functionPtr) {
+
+ struct io_props r;
+ strcpy(r.name, name);
+ r.type = type;
+ r.dimension = dimension;
+ r.importance = UNUSED;
+ r.units = units;
+ r.field = NULL;
+ r.partSize = bpartSize;
+ r.parts = NULL;
+ r.xparts = NULL;
+ r.gparts = NULL;
+ r.sparts = NULL;
+ r.bparts = bparts;
+ r.conversion = 1;
+ r.convert_part_f = NULL;
+ r.convert_part_d = NULL;
+ r.convert_part_l = NULL;
+ r.convert_gpart_f = NULL;
+ r.convert_gpart_d = NULL;
+ r.convert_gpart_l = NULL;
+ r.convert_spart_f = NULL;
+ r.convert_spart_d = NULL;
+ r.convert_spart_l = NULL;
+ r.convert_bpart_f = NULL;
+ r.convert_bpart_d = NULL;
+ r.convert_bpart_l = functionPtr;
return r;
}
diff --git a/src/part.c b/src/part.c
index ec3627d728f69f469cc7d75eb2beb9ae39ed107e..7fb27acae568a86e4f439fdac01dfe39374f9b54 100644
--- a/src/part.c
+++ b/src/part.c
@@ -96,6 +96,22 @@ void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
}
}
+/**
+ * @brief Re-link the #bpart%s associated with the list of #gpart%s.
+ *
+ * @param gparts The list of #gpart.
+ * @param N The number of particles to re-link;
+ * @param bparts The global #bpart array in which to find the #gpart offsets.
+ */
+void part_relink_bparts_to_gparts(struct gpart *gparts, size_t N,
+ struct bpart *bparts) {
+ for (size_t k = 0; k < N; k++) {
+ if (gparts[k].type == swift_type_black_hole) {
+ bparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+ }
+ }
+}
+
/**
* @brief Re-link both the #part%s and #spart%s associated with the list of
* #gpart%s.
@@ -104,14 +120,18 @@ void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
* @param N The number of particles to re-link;
* @param parts The global #part array in which to find the #gpart offsets.
* @param sparts The global #spart array in which to find the #gpart offsets.
+ * @param bparts The global #bpart array in which to find the #gpart offsets.
*/
void part_relink_all_parts_to_gparts(struct gpart *gparts, size_t N,
- struct part *parts, struct spart *sparts) {
+ struct part *parts, struct spart *sparts,
+ struct bpart *bparts) {
for (size_t k = 0; k < N; k++) {
if (gparts[k].type == swift_type_gas) {
parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
} else if (gparts[k].type == swift_type_stars) {
sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+ } else if (gparts[k].type == swift_type_black_hole) {
+ bparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
}
}
}
@@ -126,14 +146,17 @@ void part_relink_all_parts_to_gparts(struct gpart *gparts, size_t N,
* @param parts The #part array.
* @param gparts The #gpart array.
* @param sparts The #spart array.
+ * @param bparts The #bpart array.
* @param nr_parts The number of #part in the array.
* @param nr_gparts The number of #gpart in the array.
* @param nr_sparts The number of #spart in the array.
+ * @param nr_bparts The number of #bpart in the array.
* @param verbose Do we report verbosely in case of success ?
*/
void part_verify_links(struct part *parts, struct gpart *gparts,
- struct spart *sparts, size_t nr_parts, size_t nr_gparts,
- size_t nr_sparts, int verbose) {
+ struct spart *sparts, struct bpart *bparts,
+ size_t nr_parts, size_t nr_gparts, size_t nr_sparts,
+ size_t nr_bparts, int verbose) {
ticks tic = getticks();
@@ -202,6 +225,33 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
if (gparts[k].time_bin != spart->time_bin)
error("Linked particles are not at the same time !");
}
+
+ else if (gparts[k].type == swift_type_black_hole) {
+
+ /* Check that it is linked */
+ if (gparts[k].id_or_neg_offset > 0)
+ error("Black holes gpart not linked to anything !");
+
+ /* Find its link */
+ const struct bpart *bpart = &bparts[-gparts[k].id_or_neg_offset];
+
+ /* Check the reverse link */
+ if (bpart->gpart != &gparts[k]) error("Linking problem !");
+
+ /* Check that the particles are at the same place */
+ if (gparts[k].x[0] != bpart->x[0] || gparts[k].x[1] != bpart->x[1] ||
+ gparts[k].x[2] != bpart->x[2])
+ error(
+ "Linked particles are not at the same position !\n"
+ "gp->x=[%e %e %e] sp->x=[%e %e %e] diff=[%e %e %e]",
+ gparts[k].x[0], gparts[k].x[1], gparts[k].x[2], bpart->x[0],
+ bpart->x[1], bpart->x[2], gparts[k].x[0] - bpart->x[0],
+ gparts[k].x[1] - bpart->x[1], gparts[k].x[2] - bpart->x[2]);
+
+ /* Check that the particles are at the same time */
+ if (gparts[k].time_bin != bpart->time_bin)
+ error("Linked particles are not at the same time !");
+ }
}
/* Now check that all parts are linked */
@@ -250,6 +300,29 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
}
}
+ /* Now check that all bparts are linked */
+ for (size_t k = 0; k < nr_bparts; ++k) {
+
+ /* Ok, there is a link */
+ if (bparts[k].gpart != NULL) {
+
+ /* Check the link */
+ if (bparts[k].gpart->id_or_neg_offset != -(ptrdiff_t)k) {
+ error("Linking problem !");
+
+ /* Check that the particles are at the same place */
+ if (bparts[k].x[0] != bparts[k].gpart->x[0] ||
+ bparts[k].x[1] != bparts[k].gpart->x[1] ||
+ bparts[k].x[2] != bparts[k].gpart->x[2])
+ error("Linked particles are not at the same position !");
+
+ /* Check that the particles are at the same time */
+ if (bparts[k].time_bin != bparts[k].gpart->time_bin)
+ error("Linked particles are not at the same time !");
+ }
+ }
+ }
+
if (verbose) message("All links OK");
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -262,6 +335,7 @@ MPI_Datatype part_mpi_type;
MPI_Datatype xpart_mpi_type;
MPI_Datatype gpart_mpi_type;
MPI_Datatype spart_mpi_type;
+MPI_Datatype bpart_mpi_type;
/**
* @brief Registers MPI particle types.
@@ -294,5 +368,10 @@ void part_create_mpi_types(void) {
MPI_Type_commit(&spart_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for sparts.");
}
+ if (MPI_Type_contiguous(sizeof(struct bpart) / sizeof(unsigned char),
+ MPI_BYTE, &bpart_mpi_type) != MPI_SUCCESS ||
+ MPI_Type_commit(&bpart_mpi_type) != MPI_SUCCESS) {
+ error("Failed to create MPI type for bparts.");
+ }
}
#endif
diff --git a/src/part.h b/src/part.h
index 069a5075c1a8aa1037f37df2f1ce0168e1130a5f..3b2f01180d67a597b71086c46283dae2d7de49b9 100644
--- a/src/part.h
+++ b/src/part.h
@@ -40,6 +40,7 @@
#define xpart_align 128
#define spart_align 128
#define gpart_align 128
+#define bpart_align 128
/* Import the right hydro particle definition */
#if defined(MINIMAL_SPH)
@@ -103,6 +104,13 @@
#error "Invalid choice of star particle"
#endif
+/* Import the right black hole particle definition */
+#if defined(BLACK_HOLES_NONE)
+#include "./black_holes/Default/black_holes_part.h"
+#else
+#error "Invalid choice of black hole particle"
+#endif
+
void part_relink_gparts_to_parts(struct part *parts, size_t N,
ptrdiff_t offset);
void part_relink_gparts_to_sparts(struct spart *sparts, size_t N,
@@ -111,11 +119,15 @@ void part_relink_parts_to_gparts(struct gpart *gparts, size_t N,
struct part *parts);
void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
struct spart *sparts);
+void part_relink_bparts_to_gparts(struct gpart *gparts, size_t N,
+ struct bpart *bparts);
void part_relink_all_parts_to_gparts(struct gpart *gparts, size_t N,
- struct part *parts, struct spart *sparts);
+ struct part *parts, struct spart *sparts,
+ struct bpart *bparts);
void part_verify_links(struct part *parts, struct gpart *gparts,
- struct spart *sparts, size_t nr_parts, size_t nr_gparts,
- size_t nr_sparts, int verbose);
+ struct spart *sparts, struct bpart *bparts,
+ size_t nr_parts, size_t nr_gparts, size_t nr_sparts,
+ size_t nr_bparts, int verbose);
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
@@ -123,6 +135,7 @@ extern MPI_Datatype part_mpi_type;
extern MPI_Datatype xpart_mpi_type;
extern MPI_Datatype gpart_mpi_type;
extern MPI_Datatype spart_mpi_type;
+extern MPI_Datatype bpart_mpi_type;
void part_create_mpi_types(void);
#endif
diff --git a/src/runner.c b/src/runner.c
index 6f031ded7252537ea37f564634daa89b7eebda30..bf0c32ec85cfae56dad6812bea173c0b302f4bcb 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -2033,6 +2033,22 @@ void runner_do_drift_spart(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_drift_spart);
}
+/**
+ * @brief Drift all bpart in a cell.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_drift_bpart(struct runner *r, struct cell *c, int timer) {
+
+ TIMER_TIC;
+
+ cell_drift_bpart(c, r->e, 0);
+
+ if (timer) TIMER_TOC(timer_drift_bpart);
+}
+
/**
* @brief Perform the first half-kick on all the active particles in a cell.
*
@@ -2432,30 +2448,35 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
+ const int bcount = c->black_holes.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
+ struct bpart *restrict bparts = c->black_holes.parts;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
- !cell_is_active_stars(c, e)) {
+ !cell_is_active_stars(c, e) && !cell_is_active_black_holes(c, e)) {
c->hydro.updated = 0;
c->grav.updated = 0;
c->stars.updated = 0;
+ c->black_holes.updated = 0;
return;
}
- int updated = 0, g_updated = 0, s_updated = 0;
- int inhibited = 0, g_inhibited = 0, s_inhibited = 0;
+ int updated = 0, g_updated = 0, s_updated = 0, b_updated = 0;
+ int inhibited = 0, g_inhibited = 0, s_inhibited = 0, b_inhibited = 0;
integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_end_max = 0,
ti_hydro_beg_max = 0;
integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0,
ti_gravity_beg_max = 0;
integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_end_max = 0,
ti_stars_beg_max = 0;
+ integertime_t ti_black_holes_end_min = max_nr_timesteps, ti_black_holes_end_max = 0,
+ ti_black_holes_beg_max = 0;
/* No children? */
if (!c->split) {
@@ -2664,6 +2685,67 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
}
}
+
+ /* Loop over the star particles in this cell. */
+ for (int k = 0; k < bcount; k++) {
+
+ /* Get a handle on the part. */
+ struct bpart *restrict bp = &bparts[k];
+
+ /* need to be updated ? */
+ if (bpart_is_active(bp, e)) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Current end of time-step */
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current, bp->time_bin);
+
+ if (ti_end != ti_current)
+ error("Computing time-step of rogue particle.");
+#endif
+ /* Get new time-step */
+ const integertime_t ti_new_step = get_bpart_timestep(bp, e);
+
+ /* Update particle */
+ bp->time_bin = get_time_bin(ti_new_step);
+ bp->gpart->time_bin = get_time_bin(ti_new_step);
+
+ /* Number of updated s-particles */
+ b_updated++;
+ g_updated++;
+
+ ti_black_holes_end_min = min(ti_current + ti_new_step, ti_black_holes_end_min);
+ ti_black_holes_end_max = max(ti_current + ti_new_step, ti_black_holes_end_max);
+ ti_gravity_end_min = min(ti_current + ti_new_step, ti_gravity_end_min);
+ ti_gravity_end_max = max(ti_current + ti_new_step, ti_gravity_end_max);
+
+ /* What is the next starting point for this cell ? */
+ ti_black_holes_beg_max = max(ti_current, ti_black_holes_beg_max);
+ ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
+
+ /* star particle is inactive but not inhibited */
+ } else {
+
+ /* Count the number of inhibited particles */
+ if (bpart_is_inhibited(bp, e)) ++b_inhibited;
+
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current, bp->time_bin);
+
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, bp->time_bin);
+
+ ti_black_holes_end_min = min(ti_end, ti_black_holes_end_min);
+ ti_black_holes_end_max = max(ti_end, ti_black_holes_end_max);
+ ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
+ ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
+
+ /* What is the next starting point for this cell ? */
+ ti_black_holes_beg_max = max(ti_beg, ti_black_holes_beg_max);
+ ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+ }
+ }
+
} else {
/* Loop over the progeny. */
@@ -2678,18 +2760,28 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
updated += cp->hydro.updated;
g_updated += cp->grav.updated;
s_updated += cp->stars.updated;
+ b_updated += cp->black_holes.updated;
+
inhibited += cp->hydro.inhibited;
g_inhibited += cp->grav.inhibited;
s_inhibited += cp->stars.inhibited;
- ti_hydro_end_min = min(cp->hydro.ti_end_min, ti_hydro_end_min);
+ b_inhibited += cp->black_holes.inhibited;
+
+ ti_hydro_end_min = min(cp->hydro.ti_end_min, ti_hydro_end_min);
ti_hydro_end_max = max(cp->hydro.ti_end_max, ti_hydro_end_max);
ti_hydro_beg_max = max(cp->hydro.ti_beg_max, ti_hydro_beg_max);
- ti_gravity_end_min = min(cp->grav.ti_end_min, ti_gravity_end_min);
+
+ ti_gravity_end_min = min(cp->grav.ti_end_min, ti_gravity_end_min);
ti_gravity_end_max = max(cp->grav.ti_end_max, ti_gravity_end_max);
ti_gravity_beg_max = max(cp->grav.ti_beg_max, ti_gravity_beg_max);
- ti_stars_end_min = min(cp->stars.ti_end_min, ti_stars_end_min);
+
+ ti_stars_end_min = min(cp->stars.ti_end_min, ti_stars_end_min);
ti_stars_end_max = max(cp->grav.ti_end_max, ti_stars_end_max);
ti_stars_beg_max = max(cp->grav.ti_beg_max, ti_stars_beg_max);
+
+ ti_black_holes_end_min = min(cp->black_holes.ti_end_min, ti_black_holes_end_min);
+ ti_black_holes_end_max = max(cp->grav.ti_end_max, ti_black_holes_end_max);
+ ti_black_holes_beg_max = max(cp->grav.ti_beg_max, ti_black_holes_beg_max);
}
}
}
@@ -2698,9 +2790,13 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
c->hydro.updated = updated;
c->grav.updated = g_updated;
c->stars.updated = s_updated;
+ c->black_holes.updated = b_updated;
+
c->hydro.inhibited = inhibited;
c->grav.inhibited = g_inhibited;
c->stars.inhibited = s_inhibited;
+ c->black_holes.inhibited = b_inhibited;
+
c->hydro.ti_end_min = ti_hydro_end_min;
c->hydro.ti_end_max = ti_hydro_end_max;
c->hydro.ti_beg_max = ti_hydro_beg_max;
@@ -2710,6 +2806,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
c->stars.ti_end_min = ti_stars_end_min;
c->stars.ti_end_max = ti_stars_end_max;
c->stars.ti_beg_max = ti_stars_beg_max;
+ c->black_holes.ti_end_min = ti_black_holes_end_min;
+ c->black_holes.ti_end_max = ti_black_holes_end_max;
+ c->black_holes.ti_beg_max = ti_black_holes_beg_max;
#ifdef SWIFT_DEBUG_CHECKS
if (c->hydro.ti_end_min == e->ti_current &&
@@ -2721,6 +2820,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
if (c->stars.ti_end_min == e->ti_current &&
c->stars.ti_end_min < max_nr_timesteps)
error("End of next stars step is current time!");
+ if (c->black_holes.ti_end_min == e->ti_current &&
+ c->black_holes.ti_end_min < max_nr_timesteps)
+ error("End of next black holes step is current time!");
#endif
if (timer) TIMER_TOC(timer_timestep);
@@ -3449,6 +3551,9 @@ void *runner_main(void *data) {
case task_type_drift_spart:
runner_do_drift_spart(r, ci, 1);
break;
+ case task_type_drift_bpart:
+ runner_do_drift_bpart(r, ci, 1);
+ break;
case task_type_drift_gpart:
runner_do_drift_gpart(r, ci, 1);
break;
@@ -3484,6 +3589,9 @@ void *runner_main(void *data) {
if (t->subtype == task_subtype_tend_spart) {
free(t->buff);
}
+ if (t->subtype == task_subtype_tend_bpart) {
+ free(t->buff);
+ }
break;
case task_type_recv:
if (t->subtype == task_subtype_tend_part) {
@@ -3495,6 +3603,9 @@ void *runner_main(void *data) {
} else if (t->subtype == task_subtype_tend_spart) {
cell_unpack_end_step_stars(ci, (struct pcell_step_stars *)t->buff);
free(t->buff);
+ } else if (t->subtype == task_subtype_tend_bpart) {
+ cell_unpack_end_step_black_holes(ci, (struct pcell_step_black_holes *)t->buff);
+ free(t->buff);
} else if (t->subtype == task_subtype_xv) {
runner_do_recv_part(r, ci, 1, 1);
} else if (t->subtype == task_subtype_rho) {
diff --git a/src/scheduler.c b/src/scheduler.c
index 80ce572002b06991e1720db36e4f1b2fe90daa62..9ce053152590d21e07f09b98c05b499aa120cb00 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1540,6 +1540,8 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
const float gcount_j = (t->cj != NULL) ? t->cj->grav.count : 0.f;
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;
switch (t->type) {
case task_type_sort:
@@ -1633,6 +1635,9 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_drift_spart:
cost = wscale * scount_i;
break;
+ case task_type_drift_bpart:
+ cost = wscale * bcount_i;
+ break;
case task_type_init_grav:
cost = wscale * gcount_i;
break;
@@ -1658,13 +1663,13 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
cost = wscale * (count_i + scount_i);
break;
case task_type_kick1:
- cost = wscale * count_i + wscale * gcount_i;
+ cost = wscale * (count_i + gcount_i + scount_i + bcount_i);
break;
case task_type_kick2:
- cost = wscale * count_i + wscale * gcount_i;
+ cost = wscale * (count_i + gcount_i + scount_i + bcount_i);
break;
case task_type_timestep:
- cost = wscale * count_i + wscale * gcount_i;
+ cost = wscale * (count_i + gcount_i + scount_i + bcount_i);
break;
case task_type_send:
if (count_i < 1e5)
@@ -1877,6 +1882,13 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
t->buff, t->ci->mpi.pcell_size * sizeof(struct pcell_step_stars),
MPI_BYTE, t->ci->nodeID, t->flags, subtaskMPI_comms[t->subtype],
&t->req);
+ } else if (t->subtype == task_subtype_tend_bpart) {
+ t->buff = (struct pcell_step_black_holes *)malloc(
+ sizeof(struct pcell_step_black_holes) * t->ci->mpi.pcell_size);
+ err = MPI_Irecv(
+ t->buff, t->ci->mpi.pcell_size * sizeof(struct pcell_step_black_holes),
+ MPI_BYTE, t->ci->nodeID, t->flags, subtaskMPI_comms[t->subtype],
+ &t->req);
} else if (t->subtype == task_subtype_xv ||
t->subtype == task_subtype_rho ||
t->subtype == task_subtype_gradient) {
@@ -1965,6 +1977,25 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
MPI_BYTE, t->cj->nodeID, t->flags, subtaskMPI_comms[t->subtype],
&t->req);
}
+ } else if (t->subtype == task_subtype_tend_bpart) {
+ t->buff = (struct pcell_step_black_holes *)malloc(
+ sizeof(struct pcell_step_black_holes) * t->ci->mpi.pcell_size);
+ cell_pack_end_step_black_holes(t->ci, (struct pcell_step_black_holes *)t->buff);
+
+ if ((t->ci->mpi.pcell_size * sizeof(struct pcell_step_black_holes)) >
+ s->mpi_message_limit) {
+ err = MPI_Isend(
+ t->buff,
+ t->ci->mpi.pcell_size * sizeof(struct pcell_step_black_holes),
+ MPI_BYTE, t->cj->nodeID, t->flags, subtaskMPI_comms[t->subtype],
+ &t->req);
+ } else {
+ err = MPI_Issend(
+ t->buff,
+ t->ci->mpi.pcell_size * sizeof(struct pcell_step_black_holes),
+ MPI_BYTE, t->cj->nodeID, t->flags, subtaskMPI_comms[t->subtype],
+ &t->req);
+ }
} else if (t->subtype == task_subtype_xv ||
t->subtype == task_subtype_rho ||
t->subtype == task_subtype_gradient) {
diff --git a/src/single_io.c b/src/single_io.c
index 917dd880495300cea8f5d16df85513c6386b53f2..83802fca25e7ff8a325883bbbdaa3c1827dfc328 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -36,6 +36,7 @@
#include "single_io.h"
/* Local includes. */
+#include "black_holes_io.h"
#include "chemistry_io.h"
#include "common_io.h"
#include "cooling_io.h"
@@ -348,14 +349,17 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
* @param parts (output) Array of #part particles.
* @param gparts (output) Array of #gpart particles.
* @param sparts (output) Array of #spart particles.
+ * @param bparts (output) Array of #bpart particles.
* @param Ngas (output) number of Gas particles read.
* @param Ngparts (output) The number of #gpart read.
* @param Nstars (output) The number of #spart read.
+ * @param Nblackholes (output) The number of #bpart read.
* @param flag_entropy (output) 1 if the ICs contained Entropy in the
* InternalEnergy field
* @param with_hydro Are we reading gas particles ?
* @param with_gravity Are we reading/creating #gpart arrays ?
* @param with_stars Are we reading star particles ?
+ * @param with_black_hole Are we reading black hole particles ?
* @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
* @param cleanup_sqrt_a Are we cleaning-up the sqrt(a) factors in the Gadget
* IC velocities?
@@ -374,9 +378,10 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
struct part** parts, struct gpart** gparts,
- struct spart** sparts, size_t* Ngas, size_t* Ngparts,
- size_t* Nstars, int* flag_entropy, int with_hydro,
- int with_gravity, int with_stars, int cleanup_h,
+ struct spart** sparts, struct bpart** bparts, size_t* Ngas,
+ size_t* Ngparts, size_t* Nstars, size_t* Nblackholes,
+ int* flag_entropy, int with_hydro, int with_gravity,
+ int with_stars, int with_black_holes, int cleanup_h,
int cleanup_sqrt_a, double h, double a, int n_threads,
int dry_run) {
@@ -518,12 +523,22 @@ void read_ic_single(const char* fileName,
bzero(*sparts, *Nstars * sizeof(struct spart));
}
+ /* Allocate memory to store star particles */
+ if (with_black_holes) {
+ *Nblackholes = N[swift_type_black_hole];
+ if (swift_memalign("bparts", (void**)bparts, bpart_align,
+ *Nblackholes * sizeof(struct bpart)) != 0)
+ error("Error while allocating memory for stars particles");
+ bzero(*bparts, *Nblackholes * sizeof(struct bpart));
+ }
+
/* Allocate memory to store all gravity particles */
if (with_gravity) {
Ndm = N[swift_type_dark_matter];
*Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
N[swift_type_dark_matter] +
- (with_stars ? N[swift_type_stars] : 0);
+ (with_stars ? N[swift_type_stars] : 0) +
+ (with_black_holes ? N[swift_type_black_hole] : 0);
if (swift_memalign("gparts", (void**)gparts, gpart_align,
*Ngparts * sizeof(struct gpart)) != 0)
error("Error while allocating memory for gravity particles");
@@ -579,6 +594,13 @@ void read_ic_single(const char* fileName,
}
break;
+ case swift_type_black_hole:
+ if (with_black_holes) {
+ Nparticles = *Nblackholes;
+ black_holes_read_particles(*bparts, list, &num_fields);
+ }
+ break;
+
default:
message("Particle Type %d not yet supported. Particles ignored", ptype);
}
@@ -610,6 +632,11 @@ void read_ic_single(const char* fileName,
if (with_stars)
io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+ /* Duplicate the black hole particles into gparts */
+ if (with_black_holes)
+ io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
+ Ndm + *Ngas + *Nstars);
+
threadpool_clean(&tp);
}
@@ -648,6 +675,7 @@ void write_output_single(struct engine* e, const char* baseName,
const struct xpart* xparts = e->s->xparts;
const struct gpart* gparts = e->s->gparts;
const struct spart* sparts = e->s->sparts;
+ const struct bpart* bparts = e->s->bparts;
struct swift_params* params = e->parameter_file;
const int with_cosmology = e->policy & engine_policy_cosmology;
const int with_cooling = e->policy & engine_policy_cooling;
@@ -663,6 +691,7 @@ void write_output_single(struct engine* e, const char* baseName,
const size_t Ntot = e->s->nr_gparts;
const size_t Ngas = e->s->nr_parts;
const size_t Nstars = e->s->nr_sparts;
+ const size_t Nblackholes = e->s->nr_bparts;
// const size_t Nbaryons = Ngas + Nstars;
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
@@ -673,17 +702,17 @@ void write_output_single(struct engine* e, const char* baseName,
e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
const size_t Nstars_written =
e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
- const size_t Nbaryons_written = Ngas_written + Nstars_written;
+ const size_t Nblackholes_written =
+ e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
+ const size_t Nbaryons_written =
+ Ngas_written + Nstars_written + Nblackholes_written;
const size_t Ndm_written =
Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
/* Format things in a Gadget-friendly array */
- long long N_total[swift_type_count] = {(long long)Ngas_written,
- (long long)Ndm_written,
- 0,
- 0,
- (long long)Nstars_written,
- 0};
+ long long N_total[swift_type_count] = {
+ (long long)Ngas_written, (long long)Ndm_written, 0, 0,
+ (long long)Nstars_written, (long long)Nblackholes_written};
/* File name */
char fileName[FILENAME_BUFFER_SIZE];
@@ -906,6 +935,7 @@ void write_output_single(struct engine* e, const char* baseName,
struct gpart* gparts_written = NULL;
struct velociraptor_gpart_data* gpart_group_data_written = NULL;
struct spart* sparts_written = NULL;
+ struct bpart* bparts_written = NULL;
/* Write particle fields from the particle structure */
switch (ptype) {
@@ -1055,6 +1085,39 @@ void write_output_single(struct engine* e, const char* baseName,
}
} break;
+ case swift_type_black_hole: {
+ if (Nblackholes == Nblackholes_written) {
+
+ /* No inhibted particles: easy case */
+ N = Nblackholes;
+ black_holes_write_particles(bparts, list, &num_fields);
+ if (with_stf) {
+ num_fields += velociraptor_write_bparts(bparts, list + num_fields);
+ }
+ } else {
+
+ /* Ok, we need to fish out the particles we want */
+ N = Nblackholes_written;
+
+ /* Allocate temporary arrays */
+ if (swift_memalign("bparts_written", (void**)&bparts_written,
+ bpart_align,
+ Nblackholes_written * sizeof(struct bpart)) != 0)
+ error("Error while allocating temporart memory for bparts");
+
+ /* Collect the particles we want to write */
+ io_collect_bparts_to_write(bparts, bparts_written, Nblackholes,
+ Nblackholes_written);
+
+ /* Select the fields to write */
+ black_holes_write_particles(bparts_written, list, &num_fields);
+ if (with_stf) {
+ num_fields +=
+ velociraptor_write_bparts(bparts_written, list + num_fields);
+ }
+ }
+ } break;
+
default:
error("Particle Type %d not yet supported. Aborting", ptype);
}
@@ -1080,6 +1143,7 @@ void write_output_single(struct engine* e, const char* baseName,
if (gpart_group_data_written)
swift_free("gpart_group_written", gpart_group_data_written);
if (sparts_written) swift_free("sparts_written", sparts_written);
+ if (bparts_written) swift_free("bparts_written", bparts_written);
/* Close particle group */
H5Gclose(h_grp);
diff --git a/src/single_io.h b/src/single_io.h
index 62285c3da210243e76347f33780146604673656f..9ff04c893378d7f8c01621e7dbf387dc939d0157 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -33,9 +33,10 @@
void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
struct part** parts, struct gpart** gparts,
- struct spart** sparts, size_t* Ngas, size_t* Ndm,
- size_t* Nstars, int* flag_entropy, int with_hydro,
- int with_gravity, int with_stars, int cleanup_h,
+ struct spart** sparts, struct bpart** bparts, size_t* Ngas,
+ size_t* Ndm, size_t* Nstars, size_t* Nblackholes,
+ int* flag_entropy, int with_hydro, int with_gravity,
+ int with_stars, int with_black_holes, int cleanup_h,
int cleanup_sqrt_a, double h, double a, int nr_threads,
int dry_run);
diff --git a/src/space.c b/src/space.c
index a2a08be232bb8b8a439e23efb6d143d0c069771b..63f4489d27dd5ea738c67be92a171283eb3916fd 100644
--- a/src/space.c
+++ b/src/space.c
@@ -41,6 +41,7 @@
/* Local headers. */
#include "atomic.h"
+#include "black_holes.h"
#include "chemistry.h"
#include "const.h"
#include "cooling.h"
@@ -78,6 +79,9 @@ int space_extra_parts = space_extra_parts_default;
/*! Number of extra #spart we allocate memory for per top-level cell */
int space_extra_sparts = space_extra_sparts_default;
+/*! Number of extra #bpart we allocate memory for per top-level cell */
+int space_extra_bparts = space_extra_bparts_default;
+
/*! Number of extra #gpart we allocate memory for per top-level cell */
int space_extra_gparts = space_extra_gparts_default;
@@ -120,9 +124,11 @@ struct index_data {
size_t count_inhibited_part;
size_t count_inhibited_gpart;
size_t count_inhibited_spart;
+ size_t count_inhibited_bpart;
size_t count_extra_part;
size_t count_extra_gpart;
size_t count_extra_spart;
+ size_t count_extra_bpart;
};
/**
@@ -704,16 +710,19 @@ void space_allocate_extras(struct space *s, int verbose) {
size_t nr_parts = s->nr_parts;
size_t nr_gparts = s->nr_gparts;
size_t nr_sparts = s->nr_sparts;
+ size_t nr_bparts = s->nr_bparts;
/* The current number of actual particles */
size_t nr_actual_parts = nr_parts - s->nr_extra_parts;
size_t nr_actual_gparts = nr_gparts - s->nr_extra_gparts;
size_t nr_actual_sparts = nr_sparts - s->nr_extra_sparts;
+ size_t nr_actual_bparts = nr_bparts - s->nr_extra_bparts;
/* The number of particles we allocated memory for (MPI overhead) */
size_t size_parts = s->size_parts;
size_t size_gparts = s->size_gparts;
size_t size_sparts = s->size_sparts;
+ size_t size_bparts = s->size_bparts;
int *local_cells = (int *)malloc(sizeof(int) * s->nr_cells);
if (local_cells == NULL)
@@ -732,15 +741,18 @@ void space_allocate_extras(struct space *s, int verbose) {
const size_t expected_num_extra_parts = nr_local_cells * space_extra_parts;
const size_t expected_num_extra_gparts = nr_local_cells * space_extra_gparts;
const size_t expected_num_extra_sparts = nr_local_cells * space_extra_sparts;
+ const size_t expected_num_extra_bparts = nr_local_cells * space_extra_bparts;
if (verbose) {
- message("Currently have %zd/%zd/%zd real particles.", nr_actual_parts,
- nr_actual_gparts, nr_actual_sparts);
- message("Currently have %zd/%zd/%zd spaces for extra particles.",
- s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sparts);
- message("Requesting space for future %zd/%zd/%zd part/gpart/sparts.",
- expected_num_extra_parts, expected_num_extra_gparts,
- expected_num_extra_sparts);
+ message("Currently have %zd/%zd/%zd/%zd real particles.", nr_actual_parts,
+ nr_actual_gparts, nr_actual_sparts, nr_actual_bparts);
+ message("Currently have %zd/%zd/%zd/%zd spaces for extra particles.",
+ s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sparts,
+ s->nr_extra_bparts);
+ message(
+ "Requesting space for future %zd/%zd/%zd/%zd part/gpart/sparts/bparts.",
+ expected_num_extra_parts, expected_num_extra_gparts,
+ expected_num_extra_sparts, expected_num_extra_bparts);
}
if (expected_num_extra_parts < s->nr_extra_parts)
@@ -1018,11 +1030,94 @@ void space_allocate_extras(struct space *s, int verbose) {
s->nr_extra_sparts = expected_num_extra_sparts;
}
+ /* Do we have enough space for the extra bparts (i.e. we haven't used up any)
+ * ? */
+ if (nr_actual_bparts + expected_num_extra_bparts > nr_bparts) {
+
+ /* Ok... need to put some more in the game */
+
+ /* Do we need to reallocate? */
+ if (nr_actual_bparts + expected_num_extra_bparts > size_bparts) {
+
+ size_bparts = (nr_actual_bparts + expected_num_extra_bparts) *
+ engine_redistribute_alloc_margin;
+
+ if (verbose)
+ message("Re-allocating bparts array from %zd to %zd", s->size_bparts,
+ size_bparts);
+
+ /* Create more space for parts */
+ struct bpart *bparts_new = NULL;
+ if (swift_memalign("bparts", (void **)&bparts_new, bpart_align,
+ sizeof(struct bpart) * size_bparts) != 0)
+ error("Failed to allocate new bpart data");
+ memcpy(bparts_new, s->bparts, sizeof(struct bpart) * s->size_bparts);
+ swift_free("bparts", s->bparts);
+ s->bparts = bparts_new;
+
+ /* Update the counter */
+ s->size_bparts = size_bparts;
+ }
+
+ /* Turn some of the allocated spares into particles we can use */
+ for (size_t i = nr_bparts; i < nr_actual_bparts + expected_num_extra_bparts;
+ ++i) {
+ bzero(&s->bparts[i], sizeof(struct bpart));
+ s->bparts[i].time_bin = time_bin_not_created;
+ s->bparts[i].id = -42;
+ }
+
+ /* Put the spare particles in their correct cell */
+ int local_cell_id = 0;
+ int current_cell = local_cells[local_cell_id];
+ int count_in_cell = 0;
+ size_t count_extra_bparts = 0;
+ for (size_t i = 0; i < nr_actual_bparts + expected_num_extra_bparts; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (current_cell == s->nr_cells)
+ error("Cell counter beyond the maximal nr. cells.");
+#endif
+
+ if (s->bparts[i].time_bin == time_bin_not_created) {
+
+ /* We want the extra particles to be at the centre of their cell */
+ s->bparts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
+ s->bparts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
+ s->bparts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
+ ++count_in_cell;
+ count_extra_bparts++;
+ }
+
+ /* Once we have reached the number of extra bpart per cell, we move to the
+ * next */
+ if (count_in_cell == space_extra_bparts) {
+ ++local_cell_id;
+
+ if (local_cell_id == nr_local_cells) break;
+
+ current_cell = local_cells[local_cell_id];
+ count_in_cell = 0;
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (count_extra_bparts != expected_num_extra_bparts)
+ error("Constructed the wrong number of extra bparts (%zd vs. %zd)",
+ count_extra_bparts, expected_num_extra_bparts);
+#endif
+
+ /* Update the counters */
+ s->nr_bparts = nr_actual_bparts + expected_num_extra_bparts;
+ s->nr_extra_bparts = expected_num_extra_bparts;
+ }
+
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the links are correct */
- if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0))
- part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
- nr_sparts, verbose);
+ if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0) ||
+ (nr_gparts > 0 && nr_bparts > 0))
+ part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, nr_parts,
+ nr_gparts, nr_sparts, nr_bparts, verbose);
#endif
/* Free the list of local cells */
@@ -1063,26 +1158,31 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
size_t nr_parts = s->nr_parts;
size_t nr_gparts = s->nr_gparts;
size_t nr_sparts = s->nr_sparts;
+ size_t nr_bparts = s->nr_bparts;
/* The number of particles we allocated memory for */
size_t size_parts = s->size_parts;
size_t size_gparts = s->size_gparts;
size_t size_sparts = s->size_sparts;
+ size_t size_bparts = s->size_bparts;
/* Counter for the number of inhibited particles found on the node */
size_t count_inhibited_parts = 0;
size_t count_inhibited_gparts = 0;
size_t count_inhibited_sparts = 0;
+ size_t count_inhibited_bparts = 0;
/* Counter for the number of extra particles found on the node */
size_t count_extra_parts = 0;
size_t count_extra_gparts = 0;
size_t count_extra_sparts = 0;
+ size_t count_extra_bparts = 0;
/* Number of particles we expect to have after strays exchange */
const size_t h_index_size = size_parts + space_expected_max_nr_strays;
const size_t g_index_size = size_gparts + space_expected_max_nr_strays;
const size_t s_index_size = size_sparts + space_expected_max_nr_strays;
+ const size_t b_index_size = size_bparts + space_expected_max_nr_strays;
/* Allocate arrays to store the indices of the cells where particles
belong. We allocate extra space to allow for particles we may
@@ -1090,7 +1190,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
int *h_index = (int *)swift_malloc("h_index", sizeof(int) * h_index_size);
int *g_index = (int *)swift_malloc("g_index", sizeof(int) * g_index_size);
int *s_index = (int *)swift_malloc("s_index", sizeof(int) * s_index_size);
- if (h_index == NULL || g_index == NULL || s_index == NULL)
+ int *b_index = (int *)swift_malloc("b_index", sizeof(int) * b_index_size);
+ if (h_index == NULL || g_index == NULL || s_index == NULL || b_index == NULL)
error("Failed to allocate temporary particle indices.");
/* Allocate counters of particles that will land in each cell */
@@ -1100,8 +1201,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
(int *)swift_malloc("cell_gpart_counts", sizeof(int) * s->nr_cells);
int *cell_spart_counts =
(int *)swift_malloc("cell_spart_counts", sizeof(int) * s->nr_cells);
+ int *cell_bpart_counts =
+ (int *)swift_malloc("cell_bpart_counts", sizeof(int) * s->nr_cells);
if (cell_part_counts == NULL || cell_gpart_counts == NULL ||
- cell_spart_counts == NULL)
+ cell_spart_counts == NULL || cell_bpart_counts == NULL)
error("Failed to allocate cell particle count buffer.");
/* Initialise the counters, including buffer space for future particles */
@@ -1109,6 +1212,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
cell_part_counts[i] = 0;
cell_gpart_counts[i] = 0;
cell_spart_counts[i] = 0;
+ cell_bpart_counts[i] = 0;
}
/* Run through the particles and get their cell index. */
@@ -1124,6 +1228,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
space_sparts_get_cell_index(s, s_index, cell_spart_counts,
&count_inhibited_sparts, &count_extra_sparts,
verbose);
+ if (nr_bparts > 0)
+ space_bparts_get_cell_index(s, b_index, cell_bpart_counts,
+ &count_inhibited_bparts, &count_extra_bparts,
+ verbose);
#ifdef SWIFT_DEBUG_CHECKS
/* Some safety checks */
@@ -1133,6 +1241,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
error("We just repartitioned but still found inhibited sparts.");
if (repartitioned && count_inhibited_gparts)
error("We just repartitioned but still found inhibited gparts.");
+ if (repartitioned && count_inhibited_bparts)
+ error("We just repartitioned but still found inhibited bparts.");
if (count_extra_parts != s->nr_extra_parts)
error(
@@ -1146,6 +1256,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
error(
"Number of extra sparts in the spart array not matching the space "
"counter.");
+ if (count_extra_bparts != s->nr_extra_bparts)
+ error(
+ "Number of extra bparts in the bpart array not matching the space "
+ "counter.");
#endif
/* Move non-local parts and inhibited parts to the end of the list. */
@@ -1156,6 +1270,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
if (h_index[k] == -1 || cells_top[h_index[k]].nodeID != local_nodeID) {
/* One fewer particle */
+
nr_parts -= 1;
/* Swap the particle */
@@ -1248,6 +1363,55 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
error("Counts of inhibited s-particles do not match!");
#endif /* SWIFT_DEBUG_CHECKS */
+ /* Move non-local bparts and inhibited bparts to the end of the list. */
+ if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_bparts > 0)) {
+ for (size_t k = 0; k < nr_bparts; /* void */) {
+
+ /* Inhibited particle or foreign particle */
+ if (b_index[k] == -1 || cells_top[b_index[k]].nodeID != local_nodeID) {
+
+ /* One fewer particle */
+ nr_bparts -= 1;
+
+ /* Swap the particle */
+ memswap(&s->bparts[k], &s->bparts[nr_bparts], sizeof(struct bpart));
+
+ /* Swap the link with the gpart */
+ if (s->bparts[k].gpart != NULL) {
+ s->bparts[k].gpart->id_or_neg_offset = -k;
+ }
+ if (s->bparts[nr_bparts].gpart != NULL) {
+ s->bparts[nr_bparts].gpart->id_or_neg_offset = -nr_bparts;
+ }
+
+ /* Swap the index */
+ memswap(&b_index[k], &b_index[nr_bparts], sizeof(int));
+
+ } else {
+ /* Increment when not exchanging otherwise we need to retest "k".*/
+ k++;
+ }
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that all bparts are in the correct place. */
+ size_t check_count_inhibited_bpart = 0;
+ for (size_t k = 0; k < nr_bparts; k++) {
+ if (b_index[k] == -1 || cells_top[b_index[k]].nodeID != local_nodeID) {
+ error("Failed to move all non-local bparts to send list");
+ }
+ }
+ for (size_t k = nr_bparts; k < s->nr_bparts; k++) {
+ if (b_index[k] != -1 && cells_top[b_index[k]].nodeID == local_nodeID) {
+ error("Failed to remove local bparts from send list");
+ }
+ if (b_index[k] == -1) ++check_count_inhibited_bpart;
+ }
+ if (check_count_inhibited_bpart != count_inhibited_bparts)
+ error("Counts of inhibited b-particles do not match!");
+#endif /* SWIFT_DEBUG_CHECKS */
+
/* Move non-local gparts and inhibited parts to the end of the list. */
if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_gparts > 0)) {
for (size_t k = 0; k < nr_gparts; /* void */) {
@@ -1266,13 +1430,19 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_stars) {
s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+ } else if (s->gparts[k].type == swift_type_black_hole) {
+ s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
+
if (s->gparts[nr_gparts].type == swift_type_gas) {
s->parts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
} else if (s->gparts[nr_gparts].type == swift_type_stars) {
s->sparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
+ } else if (s->gparts[nr_gparts].type == swift_type_black_hole) {
+ s->bparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
+ &s->gparts[nr_gparts];
}
/* Swap the index */
@@ -1312,6 +1482,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
size_t nr_parts_exchanged = s->nr_parts - nr_parts;
size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts;
size_t nr_sparts_exchanged = s->nr_sparts - nr_sparts;
+ size_t nr_bparts_exchanged = s->nr_bparts - nr_bparts;
engine_exchange_strays(s->e, nr_parts, &h_index[nr_parts],
&nr_parts_exchanged, nr_gparts, &g_index[nr_gparts],
&nr_gparts_exchanged, nr_sparts, &s_index[nr_sparts],
@@ -1321,6 +1492,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
s->nr_parts = nr_parts + nr_parts_exchanged;
s->nr_gparts = nr_gparts + nr_gparts_exchanged;
s->nr_sparts = nr_sparts + nr_sparts_exchanged;
+ s->nr_bparts = nr_bparts + nr_bparts_exchanged;
} else {
#ifdef SWIFT_DEBUG_CHECKS
@@ -1339,6 +1511,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
cell_part_counts[k] = 0;
cell_spart_counts[k] = 0;
cell_gpart_counts[k] = 0;
+ cell_bpart_counts[k] = 0;
}
}
@@ -1364,6 +1537,17 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
s_index = sind_new;
}
+ /* Re-allocate the index array for the bparts if needed.. */
+ if (s->nr_bparts + 1 > s_index_size) {
+ int *bind_new;
+ if ((bind_new = (int *)swift_malloc(
+ "b_index", sizeof(int) * (s->nr_bparts + 1))) == NULL)
+ error("Failed to allocate temporary s-particle indices.");
+ memcpy(bind_new, b_index, sizeof(int) * nr_bparts);
+ swift_free("b_index", b_index);
+ b_index = bind_new;
+ }
+
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]};
@@ -1395,11 +1579,26 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
nr_sparts = s->nr_sparts;
+ /* Assign each received bpart to its cell. */
+ for (size_t k = nr_bparts; k < s->nr_bparts; k++) {
+ const struct bpart *const bp = &s->bparts[k];
+ b_index[k] =
+ cell_getid(cdim, bp->x[0] * ih[0], bp->x[1] * ih[1], bp->x[2] * ih[2]);
+ cell_bpart_counts[b_index[k]]++;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (cells_top[b_index[k]].nodeID != local_nodeID)
+ error("Received s-part that does not belong to me (nodeID=%i).",
+ cells_top[b_index[k]].nodeID);
+#endif
+ }
+ nr_bparts = s->nr_bparts;
+
#else /* WITH_MPI */
- /* Update the part and spart counters */
+ /* Update the part, spart and bpart counters */
s->nr_parts = nr_parts;
s->nr_sparts = nr_sparts;
+ s->nr_bparts = nr_bparts;
#endif /* WITH_MPI */
@@ -1464,6 +1663,36 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
#endif /* SWIFT_DEBUG_CHECKS */
+ /* Sort the bparts according to their cells. */
+ if (nr_bparts > 0)
+ space_bparts_sort(s->bparts, b_index, cell_bpart_counts, s->nr_cells, 0);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Verify that the bpart have been sorted correctly. */
+ for (size_t k = 0; k < nr_bparts; k++) {
+ const struct bpart *bp = &s->bparts[k];
+
+ if (bp->time_bin == time_bin_inhibited)
+ error("Inhibited particle sorted into a cell!");
+
+ /* New cell index */
+ const int new_bind =
+ cell_getid(s->cdim, bp->x[0] * s->iwidth[0], bp->x[1] * s->iwidth[1],
+ bp->x[2] * s->iwidth[2]);
+
+ /* New cell of this bpart */
+ const struct cell *c = &s->cells_top[new_bind];
+
+ if (b_index[k] != new_bind)
+ error("bpart's new cell index not matching sorted index.");
+
+ if (bp->x[0] < c->loc[0] || bp->x[0] > c->loc[0] + c->width[0] ||
+ bp->x[1] < c->loc[1] || bp->x[1] > c->loc[1] + c->width[1] ||
+ bp->x[2] < c->loc[2] || bp->x[2] > c->loc[2] + c->width[2])
+ error("bpart not sorted into the right top-level cell!");
+ }
+#endif /* SWIFT_DEBUG_CHECKS */
+
/* Extract the cell counts from the sorted indices. Deduct the extra
* particles. */
size_t last_index = 0;
@@ -1488,11 +1717,25 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
}
+ /* Extract the cell counts from the sorted indices. Deduct the extra
+ * particles. */
+ size_t last_bindex = 0;
+ b_index[nr_bparts] = s->nr_cells; // sentinel.
+ for (size_t k = 0; k < nr_bparts; k++) {
+ if (b_index[k] < b_index[k + 1]) {
+ cells_top[b_index[k]].black_holes.count =
+ k - last_bindex + 1 - space_extra_bparts;
+ last_bindex = k + 1;
+ }
+ }
+
/* We no longer need the indices as of here. */
swift_free("h_index", h_index);
swift_free("cell_part_counts", cell_part_counts);
swift_free("s_index", s_index);
swift_free("cell_spart_counts", cell_spart_counts);
+ swift_free("b_index", b_index);
+ swift_free("cell_bpart_counts", cell_bpart_counts);
#ifdef WITH_MPI
@@ -1532,10 +1775,11 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
s->nr_inhibited_parts = 0;
s->nr_inhibited_gparts = 0;
s->nr_inhibited_sparts = 0;
+ s->nr_inhibited_bparts = 0;
/* Sort the gparts according to their cells. */
if (nr_gparts > 0)
- space_gparts_sort(s->gparts, s->parts, s->sparts, g_index,
+ space_gparts_sort(s->gparts, s->parts, s->sparts, s->bparts, g_index,
cell_gpart_counts, s->nr_cells);
#ifdef SWIFT_DEBUG_CHECKS
@@ -1582,9 +1826,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the links are correct */
- if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0))
- part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
- nr_sparts, verbose);
+ if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0) ||
+ (nr_gparts > 0 && nr_bparts > 0))
+ part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, nr_parts,
+ nr_gparts, nr_sparts, nr_bparts, verbose);
#endif
/* Hook the cells up to the parts. Make list of local and non-empty cells */
@@ -1593,6 +1838,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
struct xpart *xfinger = s->xparts;
struct gpart *gfinger = s->gparts;
struct spart *sfinger = s->sparts;
+ struct bpart *bfinger = s->bparts;
s->nr_cells_with_particles = 0;
s->nr_local_cells_with_particles = 0;
s->nr_local_cells = 0;
@@ -1602,6 +1848,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
c->grav.ti_old_part = ti_current;
c->grav.ti_old_multipole = ti_current;
c->stars.ti_old_part = ti_current;
+ c->black_holes.ti_old_part = ti_current;
#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH)
c->cellID = -last_cell_id;
@@ -1609,23 +1856,26 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
#endif
const int is_local = (c->nodeID == engine_rank);
- const int has_particles =
- (c->hydro.count > 0) || (c->grav.count > 0) || (c->stars.count > 0);
+ const int has_particles = (c->hydro.count > 0) || (c->grav.count > 0) ||
+ (c->stars.count > 0) | (c->black_holes.count > 0);
if (is_local) {
c->hydro.parts = finger;
c->hydro.xparts = xfinger;
c->grav.parts = gfinger;
c->stars.parts = sfinger;
+ c->black_holes.parts = bfinger;
c->hydro.count_total = c->hydro.count + space_extra_parts;
c->grav.count_total = c->grav.count + space_extra_gparts;
c->stars.count_total = c->stars.count + space_extra_sparts;
+ c->black_holes.count_total = c->black_holes.count + space_extra_bparts;
finger = &finger[c->hydro.count_total];
xfinger = &xfinger[c->hydro.count_total];
gfinger = &gfinger[c->grav.count_total];
sfinger = &sfinger[c->stars.count_total];
+ bfinger = &bfinger[c->black_holes.count_total];
/* Add this cell to the list of local cells */
s->local_cells_top[s->nr_local_cells] = k;
@@ -2152,6 +2402,126 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
atomic_add_f(&s->sum_spart_vel_norm, sum_vel_norm);
}
+/**
+ * @brief #threadpool mapper function to compute the b-particle cell indices.
+ *
+ * @param map_data Pointer towards the b-particles.
+ * @param nr_bparts The number of b-particles to treat.
+ * @param extra_data Pointers to the space and index list
+ */
+void space_bparts_get_cell_index_mapper(void *map_data, int nr_bparts,
+ void *extra_data) {
+
+ /* Unpack the data */
+ struct bpart *restrict bparts = (struct bpart *)map_data;
+ struct index_data *data = (struct index_data *)extra_data;
+ struct space *s = data->s;
+ int *const ind = data->ind + (ptrdiff_t)(bparts - s->bparts);
+
+ /* Get some constants */
+ const double dim_x = s->dim[0];
+ const double dim_y = s->dim[1];
+ const double dim_z = s->dim[2];
+ const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+ const double ih_x = s->iwidth[0];
+ const double ih_y = s->iwidth[1];
+ const double ih_z = s->iwidth[2];
+
+ /* Init the local count buffer. */
+ int *cell_counts = (int *)calloc(sizeof(int), s->nr_cells);
+ if (cell_counts == NULL)
+ error("Failed to allocate temporary cell count buffer.");
+
+ /* Init the local collectors */
+ float min_mass = FLT_MAX;
+ float sum_vel_norm = 0.f;
+ size_t count_inhibited_bpart = 0;
+ size_t count_extra_bpart = 0;
+
+ for (int k = 0; k < nr_bparts; k++) {
+
+ /* Get the particle */
+ struct bpart *restrict bp = &bparts[k];
+
+ const double old_pos_x = bp->x[0];
+ const double old_pos_y = bp->x[1];
+ const double old_pos_z = bp->x[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (!s->periodic) {
+ if (old_pos_x < 0. || old_pos_x > dim_x)
+ error("Particle outside of volume along X.");
+ if (old_pos_y < 0. || old_pos_y > dim_y)
+ error("Particle outside of volume along Y.");
+ if (old_pos_z < 0. || old_pos_z > dim_z)
+ error("Particle outside of volume along Z.");
+ }
+#endif
+
+ /* Put it back into the simulation volume */
+ const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Get its cell index */
+ const int index =
+ cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (index < 0 || index >= cdim[0] * cdim[1] * cdim[2])
+ error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, cdim[0],
+ cdim[1], cdim[2], pos_x, pos_y, pos_z);
+
+ if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. ||
+ pos_y < 0. || pos_z < 0.)
+ error("Particle outside of simulation box. p->x=[%e %e %e]", pos_x, pos_y,
+ pos_z);
+#endif
+
+ /* Is this particle to be removed? */
+ if (bp->time_bin == time_bin_inhibited) {
+ ind[k] = -1;
+ ++count_inhibited_bpart;
+ } else if (bp->time_bin == time_bin_not_created) {
+ /* Is this a place-holder for on-the-fly creation? */
+ ind[k] = index;
+ cell_counts[index]++;
+ ++count_extra_bpart;
+ } else {
+ /* List its top-level cell index */
+ ind[k] = index;
+ cell_counts[index]++;
+
+ /* Compute minimal mass */
+ min_mass = min(min_mass, bp->mass);
+
+ /* Compute sum of velocity norm */
+ sum_vel_norm +=
+ bp->v[0] * bp->v[0] + bp->v[1] * bp->v[1] + bp->v[2] * bp->v[2];
+
+ /* Update the position */
+ bp->x[0] = pos_x;
+ bp->x[1] = pos_y;
+ bp->x[2] = pos_z;
+ }
+ }
+
+ /* Write the counts back to the global array. */
+ for (int k = 0; k < s->nr_cells; k++)
+ if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
+ free(cell_counts);
+
+ /* Write the count of inhibited and extra bparts */
+ if (count_inhibited_bpart)
+ atomic_add(&data->count_inhibited_bpart, count_inhibited_bpart);
+ if (count_extra_bpart)
+ atomic_add(&data->count_extra_bpart, count_extra_bpart);
+
+ /* Write back the minimal part mass and velocity sum */
+ atomic_min_f(&s->min_bpart_mass, min_mass);
+ atomic_add_f(&s->sum_bpart_vel_norm, sum_vel_norm);
+}
+
/**
* @brief Computes the cell index of all the particles.
*
@@ -2183,9 +2553,11 @@ void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
data.count_inhibited_part = 0;
data.count_inhibited_gpart = 0;
data.count_inhibited_spart = 0;
+ data.count_inhibited_bpart = 0;
data.count_extra_part = 0;
data.count_extra_gpart = 0;
data.count_extra_spart = 0;
+ data.count_extra_bpart = 0;
threadpool_map(&s->e->threadpool, space_parts_get_cell_index_mapper, s->parts,
s->nr_parts, sizeof(struct part), 0, &data);
@@ -2229,9 +2601,11 @@ void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
data.count_inhibited_part = 0;
data.count_inhibited_gpart = 0;
data.count_inhibited_spart = 0;
+ data.count_inhibited_bpart = 0;
data.count_extra_part = 0;
data.count_extra_gpart = 0;
data.count_extra_spart = 0;
+ data.count_extra_bpart = 0;
threadpool_map(&s->e->threadpool, space_gparts_get_cell_index_mapper,
s->gparts, s->nr_gparts, sizeof(struct gpart), 0, &data);
@@ -2275,9 +2649,11 @@ void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
data.count_inhibited_part = 0;
data.count_inhibited_gpart = 0;
data.count_inhibited_spart = 0;
+ data.count_inhibited_bpart = 0;
data.count_extra_part = 0;
data.count_extra_gpart = 0;
data.count_extra_spart = 0;
+ data.count_extra_bpart = 0;
threadpool_map(&s->e->threadpool, space_sparts_get_cell_index_mapper,
s->sparts, s->nr_sparts, sizeof(struct spart), 0, &data);
@@ -2290,6 +2666,54 @@ void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
clocks_getunit());
}
+/**
+ * @brief Computes the cell index of all the b-particles.
+ *
+ * Also computes the minimal mass of all #bpart.
+ *
+ * @param s The #space.
+ * @param bind The array of indices to fill.
+ * @param cell_counts The cell counters to update.
+ * @param count_inhibited_bparts (return) The number of #bpart to remove.
+ * @param count_extra_bparts (return) The number of #bpart for on-the-fly
+ * creation.
+ * @param verbose Are we talkative ?
+ */
+void space_bparts_get_cell_index(struct space *s, int *bind, int *cell_counts,
+ size_t *count_inhibited_bparts,
+ size_t *count_extra_bparts, int verbose) {
+
+ const ticks tic = getticks();
+
+ /* Re-set the counters */
+ s->min_bpart_mass = FLT_MAX;
+ s->sum_bpart_vel_norm = 0.f;
+
+ /* Pack the extra information */
+ struct index_data data;
+ data.s = s;
+ data.ind = bind;
+ data.cell_counts = cell_counts;
+ data.count_inhibited_part = 0;
+ data.count_inhibited_gpart = 0;
+ data.count_inhibited_spart = 0;
+ data.count_inhibited_bpart = 0;
+ data.count_extra_part = 0;
+ data.count_extra_gpart = 0;
+ data.count_extra_spart = 0;
+ data.count_extra_bpart = 0;
+
+ threadpool_map(&s->e->threadpool, space_bparts_get_cell_index_mapper,
+ s->bparts, s->nr_bparts, sizeof(struct bpart), 0, &data);
+
+ *count_inhibited_bparts = data.count_inhibited_bpart;
+ *count_extra_bparts = data.count_extra_bpart;
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
/**
* @brief Sort the particles and condensed particles according to the given
* indices.
@@ -2415,19 +2839,81 @@ void space_sparts_sort(struct spart *sparts, int *restrict ind,
swift_free("sparts_offsets", offsets);
}
+/**
+ * @brief Sort the b-particles according to the given indices.
+ *
+ * @param bparts The array of #bpart to sort.
+ * @param ind The indices with respect to which the #bpart are sorted.
+ * @param counts Number of particles per index.
+ * @param num_bins Total number of bins (length of counts).
+ * @param bparts_offset Offset of the #bpart array from the global #bpart.
+ * array.
+ */
+void space_bparts_sort(struct bpart *bparts, int *restrict ind,
+ int *restrict counts, int num_bins,
+ ptrdiff_t bparts_offset) {
+ /* Create the offsets array. */
+ size_t *offsets = NULL;
+ if (swift_memalign("bparts_offsets", (void **)&offsets,
+ SWIFT_STRUCT_ALIGNMENT,
+ sizeof(size_t) * (num_bins + 1)) != 0)
+ error("Failed to allocate temporary cell offsets array.");
+
+ offsets[0] = 0;
+ for (int k = 1; k <= num_bins; k++) {
+ offsets[k] = offsets[k - 1] + counts[k - 1];
+ counts[k - 1] = 0;
+ }
+
+ /* Loop over local cells. */
+ for (int cid = 0; cid < num_bins; cid++) {
+ for (size_t k = offsets[cid] + counts[cid]; k < offsets[cid + 1]; k++) {
+ counts[cid]++;
+ int target_cid = ind[k];
+ if (target_cid == cid) {
+ continue;
+ }
+ struct bpart temp_bpart = bparts[k];
+ while (target_cid != cid) {
+ size_t j = offsets[target_cid] + counts[target_cid]++;
+ while (ind[j] == target_cid) {
+ j = offsets[target_cid] + counts[target_cid]++;
+ }
+ memswap(&bparts[j], &temp_bpart, sizeof(struct bpart));
+ memswap(&ind[j], &target_cid, sizeof(int));
+ if (bparts[j].gpart)
+ bparts[j].gpart->id_or_neg_offset = -(j + bparts_offset);
+ }
+ bparts[k] = temp_bpart;
+ ind[k] = target_cid;
+ if (bparts[k].gpart)
+ bparts[k].gpart->id_or_neg_offset = -(k + bparts_offset);
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ for (int k = 0; k < num_bins; k++)
+ if (offsets[k + 1] != offsets[k] + counts[k])
+ error("Bad offsets after shuffle.");
+#endif /* SWIFT_DEBUG_CHECKS */
+
+ swift_free("bparts_offsets", offsets);
+}
+
/**
* @brief Sort the g-particles according to the given indices.
*
* @param gparts The array of #gpart to sort.
* @param parts Global #part array for re-linking.
* @param sparts Global #spart array for re-linking.
+ * @param bparts Global #bpart array for re-linking.
* @param ind The indices with respect to which the gparts are sorted.
* @param counts Number of particles per index.
* @param num_bins Total number of bins (length of counts).
*/
void space_gparts_sort(struct gpart *gparts, struct part *parts,
- struct spart *sparts, int *restrict ind,
- int *restrict counts, int num_bins) {
+ struct spart *sparts, struct bpart *bparts,
+ int *restrict ind, int *restrict counts, int num_bins) {
/* Create the offsets array. */
size_t *offsets = NULL;
if (swift_memalign("gparts_offsets", (void **)&offsets,
@@ -2461,6 +2947,8 @@ void space_gparts_sort(struct gpart *gparts, struct part *parts,
parts[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
} else if (gparts[j].type == swift_type_stars) {
sparts[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
+ } else if (gparts[j].type == swift_type_black_hole) {
+ bparts[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
}
}
gparts[k] = temp_gpart;
@@ -2469,6 +2957,8 @@ void space_gparts_sort(struct gpart *gparts, struct part *parts,
parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
} else if (gparts[k].type == swift_type_stars) {
sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+ } else if (gparts[k].type == swift_type_black_hole) {
+ bparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
}
}
}
@@ -2645,31 +3135,37 @@ void space_map_cells_pre(struct space *s, int full,
*/
void space_split_recursive(struct space *s, struct cell *c,
struct cell_buff *buff, struct cell_buff *sbuff,
- struct cell_buff *gbuff) {
+ struct cell_buff *bbuff, struct cell_buff *gbuff) {
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
+ const int bcount = c->black_holes.count;
const int with_self_gravity = s->with_self_gravity;
const int depth = c->depth;
int maxdepth = 0;
float h_max = 0.0f;
float stars_h_max = 0.f;
+ float black_holes_h_max = 0.f;
integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_end_max = 0,
ti_hydro_beg_max = 0;
integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0,
ti_gravity_beg_max = 0;
integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_end_max = 0,
ti_stars_beg_max = 0;
+ integertime_t ti_black_holes_end_min = max_nr_timesteps,
+ ti_black_holes_end_max = 0, ti_black_holes_beg_max = 0;
struct part *parts = c->hydro.parts;
struct gpart *gparts = c->grav.parts;
struct spart *sparts = c->stars.parts;
+ struct bpart *bparts = c->black_holes.parts;
struct xpart *xparts = c->hydro.xparts;
struct engine *e = s->e;
const integertime_t ti_current = e->ti_current;
/* If the buff is NULL, allocate it, and remember to free it. */
- const int allocate_buffer = (buff == NULL && gbuff == NULL && sbuff == NULL);
+ const int allocate_buffer =
+ (buff == NULL && gbuff == NULL && sbuff == NULL && bbuff == NULL);
if (allocate_buffer) {
if (count > 0) {
if (swift_memalign("tempbuff", (void **)&buff, SWIFT_STRUCT_ALIGNMENT,
@@ -2719,6 +3215,22 @@ void space_split_recursive(struct space *s, struct cell *c,
sbuff[k].x[2] = sparts[k].x[2];
}
}
+ if (bcount > 0) {
+ if (swift_memalign("tempbbuff", (void **)&bbuff, SWIFT_STRUCT_ALIGNMENT,
+ sizeof(struct cell_buff) * bcount) != 0)
+ error("Failed to allocate temporary indices.");
+ for (int k = 0; k < bcount; k++) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (bparts[k].time_bin == time_bin_inhibited)
+ error("Inhibited particle present in space_split()");
+ if (bparts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
+#endif
+ bbuff[k].x[0] = bparts[k].x[0];
+ bbuff[k].x[1] = bparts[k].x[1];
+ bbuff[k].x[2] = bparts[k].x[2];
+ }
+ }
}
/* Check the depth. */
@@ -2747,13 +3259,16 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->hydro.count = 0;
cp->grav.count = 0;
cp->stars.count = 0;
+ cp->black_holes.count = 0;
cp->hydro.count_total = 0;
cp->grav.count_total = 0;
cp->stars.count_total = 0;
+ cp->black_holes.count_total = 0;
cp->hydro.ti_old_part = c->hydro.ti_old_part;
cp->grav.ti_old_part = c->grav.ti_old_part;
cp->grav.ti_old_multipole = c->grav.ti_old_multipole;
cp->stars.ti_old_part = c->stars.ti_old_part;
+ cp->black_holes.ti_old_part = c->black_holes.ti_old_part;
cp->loc[0] = c->loc[0];
cp->loc[1] = c->loc[1];
cp->loc[2] = c->loc[2];
@@ -2772,6 +3287,8 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->stars.h_max = 0.f;
cp->stars.dx_max_part = 0.f;
cp->stars.dx_max_sort = 0.f;
+ cp->black_holes.h_max = 0.f;
+ cp->black_holes.dx_max_part = 0.f;
cp->nodeID = c->nodeID;
cp->parent = c;
cp->top = c->top;
@@ -2783,6 +3300,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->grav.do_sub_drift = 0;
cp->hydro.do_sub_drift = 0;
cp->stars.do_sub_drift = 0;
+ cp->black_holes.do_sub_drift = 0;
cp->hydro.do_sub_limiter = 0;
cp->hydro.do_limiter = 0;
#ifdef WITH_MPI
@@ -2794,12 +3312,12 @@ void space_split_recursive(struct space *s, struct cell *c,
}
/* Split the cell's partcle data. */
- cell_split(c, c->hydro.parts - s->parts, c->stars.parts - s->sparts, buff,
- sbuff, gbuff);
+ cell_split(c, c->hydro.parts - s->parts, c->stars.parts - s->sparts,
+ c->black_holes.parts - s->bparts, buff, sbuff, bbuff, gbuff);
/* Buffers for the progenitors */
struct cell_buff *progeny_buff = buff, *progeny_gbuff = gbuff,
- *progeny_sbuff = sbuff;
+ *progeny_sbuff = sbuff, *progeny_bbuff = bbuff;
for (int k = 0; k < 8; k++) {
@@ -2807,7 +3325,8 @@ void space_split_recursive(struct space *s, struct cell *c,
struct cell *cp = c->progeny[k];
/* Remove any progeny with zero particles. */
- if (cp->hydro.count == 0 && cp->grav.count == 0 && cp->stars.count == 0) {
+ if (cp->hydro.count == 0 && cp->grav.count == 0 && cp->stars.count == 0 &&
+ cp->black_holes.count == 0) {
space_recycle(s, cp);
c->progeny[k] = NULL;
@@ -2815,13 +3334,14 @@ void space_split_recursive(struct space *s, struct cell *c,
} else {
/* Recurse */
- space_split_recursive(s, cp, progeny_buff, progeny_sbuff,
+ space_split_recursive(s, cp, progeny_buff, progeny_sbuff, progeny_bbuff,
progeny_gbuff);
/* Update the pointers in the buffers */
progeny_buff += cp->hydro.count;
progeny_gbuff += cp->grav.count;
progeny_sbuff += cp->stars.count;
+ progeny_bbuff += cp->black_holes.count;
/* Update the cell-wide properties */
h_max = max(h_max, cp->hydro.h_max);
@@ -2835,6 +3355,12 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_stars_end_min = min(ti_stars_end_min, cp->stars.ti_end_min);
ti_stars_end_max = min(ti_stars_end_max, cp->stars.ti_end_max);
ti_stars_beg_max = min(ti_stars_beg_max, cp->stars.ti_beg_max);
+ ti_black_holes_end_min =
+ min(ti_black_holes_end_min, cp->black_holes.ti_end_min);
+ ti_black_holes_end_max =
+ min(ti_black_holes_end_max, cp->black_holes.ti_end_max);
+ ti_black_holes_beg_max =
+ min(ti_black_holes_beg_max, cp->black_holes.ti_beg_max);
/* Increase the depth */
if (cp->maxdepth > maxdepth) maxdepth = cp->maxdepth;
@@ -2961,7 +3487,9 @@ void space_split_recursive(struct space *s, struct cell *c,
timebin_t hydro_time_bin_min = num_time_bins, hydro_time_bin_max = 0;
timebin_t gravity_time_bin_min = num_time_bins, gravity_time_bin_max = 0;
- timebin_t stars_time_bin_min = num_time_bins;
+ timebin_t stars_time_bin_min = num_time_bins, stars_time_bin_max = 0;
+ timebin_t black_holes_time_bin_min = num_time_bins,
+ black_holes_time_bin_max = 0;
/* parts: Get dt_min/dt_max and h_max. */
for (int k = 0; k < count; k++) {
@@ -3003,10 +3531,8 @@ void space_split_recursive(struct space *s, struct cell *c,
if (sparts[k].time_bin == time_bin_inhibited)
error("Inhibited s-particle present in space_split()");
#endif
- gravity_time_bin_min = min(gravity_time_bin_min, sparts[k].time_bin);
- gravity_time_bin_max = max(gravity_time_bin_max, sparts[k].time_bin);
stars_time_bin_min = min(stars_time_bin_min, sparts[k].time_bin);
-
+ stars_time_bin_max = max(stars_time_bin_max, sparts[k].time_bin);
stars_h_max = max(stars_h_max, sparts[k].h);
/* Reset x_diff */
@@ -3015,6 +3541,26 @@ void space_split_recursive(struct space *s, struct cell *c,
sparts[k].x_diff[2] = 0.f;
}
+ /* bparts: Get dt_min/dt_max */
+ for (int k = 0; k < bcount; k++) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (bparts[k].time_bin == time_bin_not_created)
+ error("Extra s-particle present in space_split()");
+ if (bparts[k].time_bin == time_bin_inhibited)
+ error("Inhibited s-particle present in space_split()");
+#endif
+ black_holes_time_bin_min =
+ min(black_holes_time_bin_min, bparts[k].time_bin);
+ black_holes_time_bin_max =
+ max(black_holes_time_bin_max, bparts[k].time_bin);
+ black_holes_h_max = max(black_holes_h_max, bparts[k].h);
+
+ /* Reset x_diff */
+ bparts[k].x_diff[0] = 0.f;
+ bparts[k].x_diff[1] = 0.f;
+ bparts[k].x_diff[2] = 0.f;
+ }
+
/* Convert into integer times */
ti_hydro_end_min = get_integer_time_end(ti_current, hydro_time_bin_min);
ti_hydro_end_max = get_integer_time_end(ti_current, hydro_time_bin_max);
@@ -3025,6 +3571,15 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_gravity_beg_max =
get_integer_time_begin(ti_current + 1, gravity_time_bin_max);
ti_stars_end_min = get_integer_time_end(ti_current, stars_time_bin_min);
+ ti_stars_end_max = get_integer_time_end(ti_current, stars_time_bin_max);
+ ti_stars_beg_max =
+ get_integer_time_begin(ti_current + 1, stars_time_bin_max);
+ ti_black_holes_end_min =
+ get_integer_time_end(ti_current, black_holes_time_bin_min);
+ ti_black_holes_end_max =
+ get_integer_time_end(ti_current, black_holes_time_bin_max);
+ ti_black_holes_beg_max =
+ get_integer_time_begin(ti_current + 1, black_holes_time_bin_max);
/* Construct the multipole and the centre of mass*/
if (s->with_self_gravity) {
@@ -3066,6 +3621,10 @@ void space_split_recursive(struct space *s, struct cell *c,
c->stars.ti_end_max = ti_stars_end_max;
c->stars.ti_beg_max = ti_stars_beg_max;
c->stars.h_max = stars_h_max;
+ c->black_holes.ti_end_min = ti_black_holes_end_min;
+ c->black_holes.ti_end_max = ti_black_holes_end_max;
+ c->black_holes.ti_beg_max = ti_black_holes_beg_max;
+ c->black_holes.h_max = black_holes_h_max;
c->maxdepth = maxdepth;
/* Set ownership according to the start of the parts array. */
@@ -3075,6 +3634,9 @@ void space_split_recursive(struct space *s, struct cell *c,
else if (s->nr_sparts > 0)
c->owner = ((c->stars.parts - s->sparts) % s->nr_sparts) * s->nr_queues /
s->nr_sparts;
+ else if (s->nr_bparts > 0)
+ c->owner = ((c->black_holes.parts - s->bparts) % s->nr_bparts) *
+ s->nr_queues / s->nr_bparts;
else if (s->nr_gparts > 0)
c->owner = ((c->grav.parts - s->gparts) % s->nr_gparts) * s->nr_queues /
s->nr_gparts;
@@ -3086,6 +3648,7 @@ void space_split_recursive(struct space *s, struct cell *c,
if (buff != NULL) swift_free("tempbuff", buff);
if (gbuff != NULL) swift_free("tempgbuff", gbuff);
if (sbuff != NULL) swift_free("tempsbuff", sbuff);
+ if (bbuff != NULL) swift_free("tempbbuff", bbuff);
}
}
@@ -3107,7 +3670,7 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
/* Loop over the non-empty cells */
for (int ind = 0; ind < num_cells; ind++) {
struct cell *c = &cells_top[local_cells_with_particles[ind]];
- space_split_recursive(s, c, NULL, NULL, NULL);
+ space_split_recursive(s, c, NULL, NULL, NULL, NULL);
}
#ifdef SWIFT_DEBUG_CHECKS
@@ -3384,6 +3947,17 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
sp->x[1] = gp->x[1];
sp->x[2] = gp->x[2];
}
+
+ else if (gp->type == swift_type_black_hole) {
+
+ /* Get it's black hole friend */
+ struct bpart *bp = &s->bparts[-gp->id_or_neg_offset];
+
+ /* Synchronize positions */
+ bp->x[0] = gp->x[0];
+ bp->x[1] = gp->x[1];
+ bp->x[2] = gp->x[2];
+ }
}
}
@@ -3652,6 +4226,77 @@ void space_first_init_sparts(struct space *s, int verbose) {
clocks_getunit());
}
+void space_first_init_bparts_mapper(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct bpart *restrict bp = (struct bpart *)map_data;
+ const struct space *restrict s = (struct space *)extra_data;
+ const struct engine *e = s->e;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ const ptrdiff_t delta = bp - s->bparts;
+#endif
+
+ const float initial_h = s->initial_bpart_h;
+
+ const struct cosmology *cosmo = e->cosmology;
+ const float a_factor_vel = cosmo->a;
+
+ /* Convert velocities to internal units */
+ for (int k = 0; k < count; k++) {
+
+ bp[k].v[0] *= a_factor_vel;
+ bp[k].v[1] *= a_factor_vel;
+ bp[k].v[2] *= a_factor_vel;
+
+ /* Imposed smoothing length from parameter file */
+ if (initial_h != -1.f) {
+ bp[k].h = initial_h;
+ }
+
+#ifdef HYDRO_DIMENSION_2D
+ bp[k].x[2] = 0.f;
+ bp[k].v[2] = 0.f;
+#endif
+
+#ifdef HYDRO_DIMENSION_1D
+ bp[k].x[1] = bp[k].x[2] = 0.f;
+ bp[k].v[1] = bp[k].v[2] = 0.f;
+#endif
+ }
+
+ /* Initialise the rest */
+ for (int k = 0; k < count; k++) {
+
+ black_holes_first_init_bpart(&bp[k]);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (bp[k].gpart && bp[k].gpart->id_or_neg_offset != -(k + delta))
+ error("Invalid gpart -> bpart link");
+
+ /* Initialise the time-integration check variables */
+ bp[k].ti_drift = 0;
+ bp[k].ti_kick = 0;
+#endif
+ }
+}
+
+/**
+ * @brief Initialises all the b-particles by setting them into a valid state
+ *
+ * Calls stars_first_init_bpart() on all the particles
+ */
+void space_first_init_bparts(struct space *s, int verbose) {
+ const ticks tic = getticks();
+ if (s->nr_bparts > 0)
+ threadpool_map(&s->e->threadpool, space_first_init_bparts_mapper, s->bparts,
+ s->nr_bparts, sizeof(struct bpart), 0, s);
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
void space_init_parts_mapper(void *restrict map_data, int count,
void *restrict extra_data) {
@@ -3730,6 +4375,32 @@ void space_init_sparts(struct space *s, int verbose) {
clocks_getunit());
}
+void space_init_bparts_mapper(void *restrict map_data, int bcount,
+ void *restrict extra_data) {
+
+ struct bpart *restrict bparts = (struct bpart *)map_data;
+ for (int k = 0; k < bcount; k++) black_holes_init_bpart(&bparts[k]);
+}
+
+/**
+ * @brief Calls the #bpart initialisation function on all particles in the
+ * space.
+ *
+ * @param s The #space.
+ * @param verbose Are we talkative?
+ */
+void space_init_bparts(struct space *s, int verbose) {
+
+ const ticks tic = getticks();
+
+ if (s->nr_bparts > 0)
+ threadpool_map(&s->e->threadpool, space_init_bparts_mapper, s->bparts,
+ s->nr_bparts, sizeof(struct bpart), 0, NULL);
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
void space_convert_quantities_mapper(void *restrict map_data, int count,
void *restrict extra_data) {
struct space *s = (struct space *)extra_data;
@@ -3776,9 +4447,11 @@ void space_convert_quantities(struct space *s, int verbose) {
* @param parts Array of Gas particles.
* @param gparts Array of Gravity particles.
* @param sparts Array of stars particles.
+ * @param bparts Array of black hole particles.
* @param Npart The number of Gas particles in the space.
* @param Ngpart The number of Gravity particles in the space.
* @param Nspart The number of stars particles in the space.
+ * @param Nbpart The number of black hole particles in the space.
* @param periodic flag whether the domain is periodic or not.
* @param replicate How many replications along each direction do we want?
* @param generate_gas_in_ics Are we generating gas particles from the gparts?
@@ -3796,10 +4469,10 @@ void space_convert_quantities(struct space *s, int verbose) {
void space_init(struct space *s, struct swift_params *params,
const struct cosmology *cosmo, double dim[3],
struct part *parts, struct gpart *gparts, struct spart *sparts,
- size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
- int replicate, int generate_gas_in_ics, int hydro,
- int self_gravity, int star_formation, int verbose,
- int dry_run) {
+ struct bpart *bparts, size_t Npart, size_t Ngpart,
+ size_t Nspart, size_t Nbpart, int periodic, int replicate,
+ int generate_gas_in_ics, int hydro, int self_gravity,
+ int star_formation, int verbose, int dry_run) {
/* Clean-up everything */
bzero(s, sizeof(struct space));
@@ -3815,24 +4488,31 @@ void space_init(struct space *s, struct swift_params *params,
s->nr_parts = Npart;
s->nr_gparts = Ngpart;
s->nr_sparts = Nspart;
+ s->nr_bparts = Nbpart;
s->size_parts = Npart;
s->size_gparts = Ngpart;
s->size_sparts = Nspart;
+ s->size_bparts = Nbpart;
s->nr_inhibited_parts = 0;
s->nr_inhibited_gparts = 0;
s->nr_inhibited_sparts = 0;
+ s->nr_inhibited_bparts = 0;
s->nr_extra_parts = 0;
s->nr_extra_gparts = 0;
s->nr_extra_sparts = 0;
+ s->nr_extra_bparts = 0;
s->parts = parts;
s->gparts = gparts;
s->sparts = sparts;
+ s->bparts = bparts;
s->min_part_mass = FLT_MAX;
s->min_gpart_mass = FLT_MAX;
s->min_spart_mass = FLT_MAX;
+ s->min_bpart_mass = FLT_MAX;
s->sum_part_vel_norm = 0.f;
s->sum_gpart_vel_norm = 0.f;
s->sum_spart_vel_norm = 0.f;
+ s->sum_bpart_vel_norm = 0.f;
s->nr_queues = 1; /* Temporary value until engine construction */
/* Are we generating gas from the DM-only ICs? */
@@ -3844,7 +4524,8 @@ void space_init(struct space *s, struct swift_params *params,
Ngpart = s->nr_gparts;
#ifdef SWIFT_DEBUG_CHECKS
- part_verify_links(parts, gparts, sparts, Npart, Ngpart, Nspart, 1);
+ part_verify_links(parts, gparts, sparts, bparts, Npart, Ngpart, Nspart,
+ Nbpart, 1);
#endif
}
@@ -3857,12 +4538,15 @@ void space_init(struct space *s, struct swift_params *params,
parts = s->parts;
gparts = s->gparts;
sparts = s->sparts;
+ bparts = s->bparts;
Npart = s->nr_parts;
Ngpart = s->nr_gparts;
Nspart = s->nr_sparts;
+ Nbpart = s->nr_bparts;
#ifdef SWIFT_DEBUG_CHECKS
- part_verify_links(parts, gparts, sparts, Npart, Ngpart, Nspart, 1);
+ part_verify_links(parts, gparts, sparts, bparts, Npart, Ngpart, Nspart,
+ Nbpart, 1);
#endif
}
@@ -3908,6 +4592,8 @@ void space_init(struct space *s, struct swift_params *params,
params, "Scheduler:cell_extra_sparts", space_extra_sparts_default);
space_extra_gparts = parser_get_opt_param_int(
params, "Scheduler:cell_extra_gparts", space_extra_gparts_default);
+ space_extra_bparts = parser_get_opt_param_int(
+ params, "Scheduler:cell_extra_bparts", space_extra_bparts_default);
if (verbose) {
message("max_size set to %d split_size set to %d", space_maxsize,
@@ -3933,6 +4619,12 @@ void space_init(struct space *s, struct swift_params *params,
if (s->initial_spart_h != -1.f) {
message("Imposing a star smoothing length of %e", s->initial_spart_h);
}
+ /* Read in imposed star smoothing length */
+ s->initial_bpart_h = parser_get_opt_param_float(
+ params, "InitialConditions:black_holes_smoothing_length", -1.f);
+ if (s->initial_bpart_h != -1.f) {
+ message("Imposing a BH smoothing length of %e", s->initial_bpart_h);
+ }
/* Apply shift */
double shift[3] = {0.0, 0.0, 0.0};
@@ -3955,6 +4647,11 @@ void space_init(struct space *s, struct swift_params *params,
sparts[k].x[1] += shift[1];
sparts[k].x[2] += shift[2];
}
+ for (size_t k = 0; k < Nbpart; k++) {
+ bparts[k].x[0] += shift[0];
+ bparts[k].x[1] += shift[1];
+ bparts[k].x[2] += shift[2];
+ }
}
if (!dry_run) {
@@ -4000,6 +4697,20 @@ void space_init(struct space *s, struct swift_params *params,
if (sparts[k].x[j] < 0 || sparts[k].x[j] >= s->dim[j])
error("Not all s-particles are within the specified domain.");
}
+
+ /* Same for the bparts */
+ if (periodic) {
+ for (size_t k = 0; k < Nbpart; k++)
+ for (int j = 0; j < 3; j++) {
+ while (bparts[k].x[j] < 0) bparts[k].x[j] += s->dim[j];
+ while (bparts[k].x[j] >= s->dim[j]) bparts[k].x[j] -= s->dim[j];
+ }
+ } else {
+ for (size_t k = 0; k < Nbpart; k++)
+ for (int j = 0; j < 3; j++)
+ if (bparts[k].x[j] < 0 || bparts[k].x[j] >= s->dim[j])
+ error("Not all b-particles are within the specified domain.");
+ }
}
/* Allocate the extra parts array for the gas particles. */
@@ -4143,8 +4854,8 @@ void space_replicate(struct space *s, int replicate, int verbose) {
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that everything is correct */
- part_verify_links(s->parts, s->gparts, s->sparts, s->nr_parts, s->nr_gparts,
- s->nr_sparts, verbose);
+ part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
+ s->nr_gparts, s->nr_sparts, s->nr_bparts, verbose);
#endif
}
@@ -4504,6 +5215,9 @@ void space_struct_dump(struct space *s, FILE *stream) {
if (s->nr_sparts > 0)
restart_write_blocks(s->sparts, s->nr_sparts, sizeof(struct spart), stream,
"sparts", "sparts");
+ if (s->nr_bparts > 0)
+ restart_write_blocks(s->bparts, s->nr_bparts, sizeof(struct bpart), stream,
+ "bparts", "bparts");
}
/**
@@ -4575,7 +5289,16 @@ void space_struct_restore(struct space *s, FILE *stream) {
restart_read_blocks(s->sparts, s->nr_sparts, sizeof(struct spart), stream,
NULL, "sparts");
}
-
+ s->bparts = NULL;
+ if (s->nr_bparts > 0) {
+ if (swift_memalign("bparts", (void **)&s->bparts, bpart_align,
+ s->size_bparts * sizeof(struct bpart)) != 0)
+ error("Failed to allocate restore bpart array.");
+
+ restart_read_blocks(s->bparts, s->nr_bparts, sizeof(struct bpart), stream,
+ NULL, "bparts");
+ }
+
/* Need to reconnect the gravity parts to their hydro and stars particles. */
/* Re-link the parts. */
if (s->nr_parts > 0 && s->nr_gparts > 0)
@@ -4585,10 +5308,14 @@ void space_struct_restore(struct space *s, FILE *stream) {
if (s->nr_sparts > 0 && s->nr_gparts > 0)
part_relink_sparts_to_gparts(s->gparts, s->nr_gparts, s->sparts);
+ /* Re-link the bparts. */
+ if (s->nr_bparts > 0 && s->nr_gparts > 0)
+ part_relink_bparts_to_gparts(s->gparts, s->nr_gparts, s->bparts);
+
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that everything is correct */
- part_verify_links(s->parts, s->gparts, s->sparts, s->nr_parts, s->nr_gparts,
- s->nr_sparts, 1);
+ part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
+ s->nr_gparts, s->nr_sparts, s->nr_bparts, 1);
#endif
}
diff --git a/src/space.h b/src/space.h
index 272ee41c3569cf3580678c0de885a41075f0e1b4..19acb8bffb59e57a14e4d2b138cb9a76508b9708 100644
--- a/src/space.h
+++ b/src/space.h
@@ -48,6 +48,7 @@ struct cosmology;
#define space_extra_parts_default 0
#define space_extra_gparts_default 0
#define space_extra_sparts_default 100
+#define space_extra_bparts_default 0
#define space_expected_max_nr_strays_default 100
#define space_subsize_pair_hydro_default 256000000
#define space_subsize_self_hydro_default 32000
@@ -74,6 +75,7 @@ extern int space_subdepth_diff_grav;
extern int space_extra_parts;
extern int space_extra_gparts;
extern int space_extra_sparts;
+extern int space_extra_bparts;
/**
* @brief The space in which the cells and particles reside.
@@ -167,6 +169,9 @@ struct space {
/*! The total number of #spart in the space. */
size_t nr_sparts;
+ /*! The total number of #bpart in the space. */
+ size_t nr_bparts;
+
/*! The total number of #part we allocated memory for */
size_t size_parts;
@@ -176,6 +181,9 @@ struct space {
/*! The total number of #spart we allocated memory for */
size_t size_sparts;
+ /*! The total number of #bpart we allocated memory for */
+ size_t size_bparts;
+
/*! Number of inhibted gas particles in the space */
size_t nr_inhibited_parts;
@@ -185,6 +193,9 @@ struct space {
/*! Number of inhibted star particles in the space */
size_t nr_inhibited_sparts;
+ /*! Number of inhibted black hole particles in the space */
+ size_t nr_inhibited_bparts;
+
/*! Number of extra #part we allocated (for on-the-fly creation) */
size_t nr_extra_parts;
@@ -194,6 +205,9 @@ struct space {
/*! Number of extra #spart we allocated (for on-the-fly creation) */
size_t nr_extra_sparts;
+ /*! Number of extra #bpart we allocated (for on-the-fly creation) */
+ size_t nr_extra_bparts;
+
/*! The particle data (cells have pointers to this). */
struct part *parts;
@@ -206,6 +220,9 @@ struct space {
/*! The s-particle data (cells have pointers to this). */
struct spart *sparts;
+ /*! The b-particle data (cells have pointers to this). */
+ struct bpart *bparts;
+
/*! Minimal mass of all the #part */
float min_part_mass;
@@ -215,6 +232,9 @@ struct space {
/*! Minimal mass of all the #spart */
float min_spart_mass;
+ /*! Minimal mass of all the #bpart */
+ float min_bpart_mass;
+
/*! Sum of the norm of the velocity of all the #part */
float sum_part_vel_norm;
@@ -224,9 +244,15 @@ struct space {
/*! Sum of the norm of the velocity of all the #spart */
float sum_spart_vel_norm;
+ /*! Sum of the norm of the velocity of all the #bpart */
+ float sum_bpart_vel_norm;
+
/*! Initial value of the smoothing length read from the parameter file */
float initial_spart_h;
+ /*! Initial value of the smoothing length read from the parameter file */
+ float initial_bpart_h;
+
/*! General-purpose lock for this space. */
swift_lock_type lock;
@@ -261,16 +287,19 @@ void space_free_buff_sort_indices(struct space *s);
void space_parts_sort(struct part *parts, struct xpart *xparts, int *ind,
int *counts, int num_bins, ptrdiff_t parts_offset);
void space_gparts_sort(struct gpart *gparts, struct part *parts,
- struct spart *sparts, int *ind, int *counts,
- int num_bins);
+ struct spart *sparts, struct bpart *bparts, int *ind,
+ int *counts, int num_bins);
void space_sparts_sort(struct spart *sparts, int *ind, int *counts,
int num_bins, ptrdiff_t sparts_offset);
+void space_bparts_sort(struct bpart *bparts, int *ind, int *counts,
+ int num_bins, ptrdiff_t bparts_offset);
void space_getcells(struct space *s, int nr_cells, struct cell **cells);
void space_init(struct space *s, struct swift_params *params,
const struct cosmology *cosmo, double dim[3],
struct part *parts, struct gpart *gparts, struct spart *sparts,
- size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
- int replicate, int generate_gas_in_ics, int hydro, int gravity,
+ struct bpart *bparts, size_t Npart, size_t Ngpart,
+ size_t Nspart, size_t Nbpart, int periodic, int replicate,
+ int generate_gas_in_ics, int hydro, int gravity,
int star_formation, int verbose, int dry_run);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
@@ -302,16 +331,18 @@ void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
size_t *count_inhibited_sparts,
size_t *count_extra_sparts, int verbose);
+void space_bparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
+ size_t *count_inhibited_bparts,
+ size_t *count_extra_bparts, int verbose);
void space_synchronize_particle_positions(struct space *s);
-void space_do_parts_sort(void);
-void space_do_gparts_sort(void);
-void space_do_sparts_sort(void);
void space_first_init_parts(struct space *s, int verbose);
void space_first_init_gparts(struct space *s, int verbose);
void space_first_init_sparts(struct space *s, int verbose);
+void space_first_init_bparts(struct space *s, int verbose);
void space_init_parts(struct space *s, int verbose);
void space_init_gparts(struct space *s, int verbose);
void space_init_sparts(struct space *s, int verbose);
+void space_init_bparts(struct space *s, int verbose);
void space_convert_quantities(struct space *s, int verbose);
void space_link_cleanup(struct space *s);
void space_check_drift_point(struct space *s, integertime_t ti_drift,
diff --git a/src/task.c b/src/task.c
index a123809bca797b812d6acd94bce3fd8b7a05c1f5..0c19120e5a0af72ca4d39adac6a7d4e2e0bf9e14 100644
--- a/src/task.c
+++ b/src/task.c
@@ -62,6 +62,7 @@ const char *taskID_names[task_type_count] = {"none",
"extra_ghost",
"drift_part",
"drift_spart",
+ "drift_bpart",
"drift_gpart",
"drift_gpart_out",
"end_hydro_force",
@@ -100,6 +101,7 @@ const char *subtaskID_names[task_subtype_count] = {"none",
"tend_part",
"tend_gpart",
"tend_spart",
+ "tend_bpart",
"xv",
"rho",
"gpart",
diff --git a/src/task.h b/src/task.h
index 7bfd49a79904b73a87f3f5a6f0b0176e596d3707..619bca0bd2b88ff9abc68a0f61eb31287a06bbb7 100644
--- a/src/task.h
+++ b/src/task.h
@@ -53,6 +53,7 @@ enum task_types {
task_type_extra_ghost,
task_type_drift_part,
task_type_drift_spart,
+ task_type_drift_bpart,
task_type_drift_gpart,
task_type_drift_gpart_out, /* Implicit */
task_type_end_hydro_force,
@@ -96,6 +97,7 @@ enum task_subtypes {
task_subtype_tend_part,
task_subtype_tend_gpart,
task_subtype_tend_spart,
+ task_subtype_tend_bpart,
task_subtype_xv,
task_subtype_rho,
task_subtype_gpart,
diff --git a/src/timestep.h b/src/timestep.h
index b98ce06d5e69a2e2b5cb8503322c025dc69f92c7..4b252467fbd2ec08a85ac509f480ffcdac555dca 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -230,4 +230,52 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
return new_dti;
}
+/**
+ * @brief Compute the new (integer) time-step of a given #bpart
+ *
+ * @param bp The #bpart.
+ * @param e The #engine (used to get some constants).
+ */
+__attribute__((always_inline)) INLINE static integertime_t get_bpart_timestep(
+ const struct bpart *restrict bp, const struct engine *restrict e) {
+
+ /* Stellar time-step */
+ float new_dt_black_holes = black_holes_compute_timestep(bp);
+
+ /* Gravity time-step */
+ float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;
+
+ if (e->policy & engine_policy_external_gravity)
+ new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
+ e->physical_constants, bp->gpart);
+
+ const float a_hydro[3] = {0.f, 0.f, 0.f};
+ if (e->policy & engine_policy_self_gravity)
+ new_dt_self = gravity_compute_timestep_self(
+ bp->gpart, a_hydro, e->gravity_properties, e->cosmology);
+
+ /* Take the minimum of all */
+ float new_dt = min3(new_dt_black_holes, new_dt_self, new_dt_ext);
+
+ /* Apply the maximal dibslacement constraint (FLT_MAX if non-cosmological)*/
+ new_dt = min(new_dt, e->dt_max_RMS_displacement);
+
+ /* Apply cosmology correction (This is 1 if non-cosmological) */
+ new_dt *= e->cosmology->time_step_factor;
+
+ /* Limit timestep within the allowed range */
+ new_dt = min(new_dt, e->dt_max);
+ if (new_dt < e->dt_min) {
+ error("bpart (id=%lld) wants a time-step (%e) below dt_min (%e)", bp->id,
+ new_dt, e->dt_min);
+ }
+
+ /* Convert to integer time */
+ const integertime_t new_dti = make_integer_timestep(
+ new_dt, bp->time_bin, e->ti_current, e->time_base_inv);
+
+ return new_dti;
+}
+
+
#endif /* SWIFT_TIMESTEP_H */
diff --git a/src/velociraptor_io.h b/src/velociraptor_io.h
index f18398219bfbc5cd6bb58a37b103f29527fa5589..d535e54815139e243b9a3bc40ec8dd4de2af1ac1 100644
--- a/src/velociraptor_io.h
+++ b/src/velociraptor_io.h
@@ -45,6 +45,17 @@ INLINE static void velociraptor_convert_spart_groupID(const struct engine* e,
}
}
+INLINE static void velociraptor_convert_bpart_groupID(const struct engine* e,
+ const struct bpart* bp,
+ long long* ret) {
+ if (bp->gpart == NULL)
+ ret[0] = 0.f;
+ else {
+ const ptrdiff_t offset = bp->gpart - e->s->gparts;
+ *ret = (e->s->gpart_group_data + offset)->groupID;
+ }
+}
+
__attribute__((always_inline)) INLINE static int velociraptor_write_parts(
const struct part* parts, const struct xpart* xparts,
struct io_props* list) {
@@ -75,4 +86,14 @@ __attribute__((always_inline)) INLINE static int velociraptor_write_sparts(
return 1;
}
+__attribute__((always_inline)) INLINE static int velociraptor_write_bparts(
+ const struct bpart* bparts, struct io_props* list) {
+
+ list[0] = io_make_output_field_convert_bpart(
+ "GroupID", LONGLONG, 1, UNIT_CONV_NO_UNITS, bparts,
+ velociraptor_convert_bpart_groupID);
+
+ return 1;
+}
+
#endif /* SWIFT_VELOCIRAPTOR_IO_H */