Added basic support for black hole particles. Particles can be read in,...

Added basic support for black hole particles. Particles can be read in, attached to gpart, integrated forward in time, drifted and written to snapshots.

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.
     

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

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)...",

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 */

src/black_holes.h

+/*******************************************************************************
+ * 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

src/black_holes/Default/black_holes.h

+/*******************************************************************************
+ * 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 */

src/black_holes/Default/black_holes_io.h

+/*******************************************************************************
+ * 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,