Added file infrastructure for the EAGLE black hole model.

configure.ac

@@ -1340,7 +1340,7 @@ case "$with_subgrid" in
 	with_subgrid_stars=EAGLE
 	with_subgrid_star_formation=EAGLE
 	with_subgrid_feedback=EAGLE
-	with_subgrid_black_holes=none
+	with_subgrid_black_holes=EAGLE
    ;;
    *)
       AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
@@ -1789,6 +1789,9 @@ case "$with_black_holes" in
    none)
       AC_DEFINE([BLACK_HOLES_NONE], [1], [No black hole model])
    ;;
+   EAGLE)
+      AC_DEFINE([BLACK_HOLES_EAGLE], [1], [EAGLE black hole model])
+   ;;
    *)
       AC_MSG_ERROR([Unknown stellar model: $with_black_holes])
    ;;

src/Makefile.am

@@ -196,7 +196,9 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  feedback/EAGLE/feedback_properties.h feedback/EAGLE/imf.h feedback/EAGLE/interpolate.h \
                  feedback/EAGLE/yield_tables.h \
                  black_holes/Default/black_holes.h black_holes/Default/black_holes_io.h \
-		 black_holes/Default/black_holes_part.h 
+		 black_holes/Default/black_holes_part.h black_holes/Default/black_holes_iact.h \
+                 black_holes/EAGLE/black_holes.h black_holes/EAGLE/black_holes_io.h \
+		 black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h
 
 
 # Sources and flags for regular library

src/black_holes.h

@@ -26,6 +26,9 @@
 #if defined(BLACK_HOLES_NONE)
 #include "./black_holes/Default/black_holes.h"
 #include "./black_holes/Default/black_holes_iact.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes.h"
+#include "./black_holes/EAGLE/black_holes_iact.h"
 #else
 #error "Invalid choice of star model"
 #endif

src/black_holes/EAGLE/#black_holes_part.h#

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_EAGLE_BLACK_HOLE_PART_H
+#define SWIFT_EAGLE_BLACK_HOLE_PART_H
+
+/* 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_EAGLE_BLACK_HOLE_PART_H */

src/black_holes/EAGLE/.#black_holes_part.h

+matthieu@laptop.2863:1556535537
\ No newline at end of file

src/black_holes/EAGLE/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_EAGLE_BLACK_HOLES_H
+#define SWIFT_EAGLE_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 black hole particle.
+ *
+ * @param bp Pointer to the s-particle data.
+ */
+__attribute__((always_inline)) INLINE static float black_holes_compute_timestep(
+    const struct bpart* const bp) {
+
+  return FLT_MAX;
+}
+
+/**
+ * @brief Initialises the b-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param bp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
+    struct bpart* bp) {
+
+  bp->time_bin = 0;
+}
+
+/**
+ * @brief Prepares a b-particle for its interactions
+ *
+ * @param bp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_init_bpart(
+    struct bpart* bp) {
+
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+    bp->ids_ngbs_density[i] = -1;
+  bp->num_ngb_density = 0;
+#endif
+
+  bp->density.wcount = 0.f;
+  bp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param bp The particle
+ * @param dt_drift The drift time-step for positions.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_predict_extra(
+    struct bpart* restrict bp, float dt_drift) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param bp The particle.
+ */
+__attribute__((always_inline)) INLINE static void
+black_holes_reset_predicted_values(struct bpart* restrict bp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param bp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void black_holes_kick_extra(
+    struct bpart* bp, float dt) {}
+
+/**
+ * @brief Finishes the calculation of density on black holes
+ *
+ * @param bp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_end_density(
+    struct bpart* bp, const struct cosmology* cosmo) {
+
+  /* Some smoothing length multiples. */
+  const float h = bp->h;
+  const float h_inv = 1.0f / h;                       /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
+  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+  /* Finish the calculation by inserting the missing h-factors */
+  bp->density.wcount *= h_inv_dim;
+  bp->density.wcount_dh *= h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #spart has 0
+ * ngbs.
+ *
+ * @param bp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+black_holes_bpart_has_no_neighbours(struct bpart* restrict bp,
+                                    const struct cosmology* cosmo) {
+
+  /* Some smoothing length multiples. */
+  const float h = bp->h;
+  const float h_inv = 1.0f / h;                 /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+  /* Re-set problematic values */
+  bp->density.wcount = kernel_root * h_inv_dim;
+  bp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is the equivalent of hydro_reset_acceleration.
+ * We do not compute the acceleration on black hole, therefore no need to use
+ * it.
+ *
+ * @param bp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void black_holes_reset_feedback(
+    struct bpart* restrict bp) {
+
+#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
+  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+    bp->ids_ngbs_force[i] = -1;
+  bp->num_ngb_force = 0;
+#endif
+}
+
+#endif /* SWIFT_EAGLE_BLACK_HOLES_H */

src/black_holes/EAGLE/black_holes_iact.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_EAGLE_BH_IACT_H
+#define SWIFT_EAGLE_BH_IACT_H
+
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param bi First bparticle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
+    const float r2, const float *dx, const float hi, const float hj,
+    struct bpart *restrict bi, const struct part *restrict pj, const float a,
+    const float H) {
+
+  float wi, wi_dx;
+
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
+
+  /* Compute the kernel function */
+  const float hi_inv = 1.0f / hi;
+  const float ui = r * hi_inv;
+  kernel_deval(ui, &wi, &wi_dx);
+
+  /* Compute contribution to the number of neighbours */
+  bi->density.wcount += wi;
+  bi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+#ifdef DEBUG_INTERACTIONS_BH
+  /* Update ngb counters */
+  if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_BH)
+    bi->ids_ngbs_density[si->num_ngb_density] = pj->id;
+
+  /* Update ngb counters */
+  ++si->num_ngb_density;
+#endif
+}
+
+/**
+ * @brief Feedback interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param bi First bparticle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_bh_feedback(const float r2, const float *dx, const float hi,
+                               const float hj, struct bpart *restrict bi,
+                               struct part *restrict pj, const float a,
+                               const float H) {
+#ifdef DEBUG_INTERACTIONS_BH
+  /* Update ngb counters */
+  if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
+    bi->ids_ngbs_force[si->num_ngb_force] = pj->id;
+
+  /* Update ngb counters */
+  ++si->num_ngb_force;
+#endif
+}
+
+#endif /* SWIFT_EAGLE_BH_IACT_H */

src/black_holes/EAGLE/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_EAGLE_BLACK_HOLES_IO_H
+#define SWIFT_EAGLE_BLACK_HOLES_IO_H
+
+#include "black_holes_part.h"
+#include "io_properties.h"
+
+/**
+ * @brief Specifies which b-particle fields to read from a dataset
+ *
+ * @param bparts The b-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void black_holes_read_particles(struct bpart *bparts,
+                                              struct io_props *list,
+                                              int *num_fields) {
+
+  /* Say how much we want to read */
+  *num_fields = 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_EAGLE_BLACK_HOLES_IO_H */

src/black_holes/EAGLE/black_holes_part.h

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

src/black_holes_io.h

@@ -24,6 +24,8 @@
 /* Load the correct star type */
 #if defined(BLACK_HOLES_NONE)
 #include "./black_holes/Default/black_holes_io.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes_io.h"
 #else
 #error "Invalid choice of star model"
 #endif

src/part.h

@@ -109,6 +109,8 @@
 /* Import the right black hole particle definition */
 #if defined(BLACK_HOLES_NONE)
 #include "./black_holes/Default/black_holes_part.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes_part.h"
 #else
 #error "Invalid choice of black hole particle"
 #endif