Added neighbour finding loops for the black holes and infrastrcture to create the tasks.

examples/main.c

@@ -500,6 +500,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. */

src/black_holes.h

@@ -25,6 +25,7 @@
 /* Select the correct star model */
 #if defined(BLACK_HOLES_NONE)
 #include "./black_holes/Default/black_holes.h"
+#include "./black_holes/Default/black_holes_iact.h"
 #else
 #error "Invalid choice of star model"
 #endif

src/black_holes/Default/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_DEFAULT_BH_IACT_H
+#define SWIFT_DEFAULT_BH_IACT_H
+
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param bi First 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

src/cell.c

@@ -1786,6 +1786,8 @@ void cell_clean_links(struct cell *c, void *data) {
   c->grav.mm = NULL;
   c->stars.density = NULL;
   c->stars.feedback = NULL;
+  c->black_holes.density = NULL;
+  c->black_holes.feedback = NULL;
 }
 
 /**

src/cell.h

@@ -606,6 +606,15 @@ struct cell {
     /*! The drift task for bparts */
     struct task *drift;
 
+    /*! The star ghost task itself */
+    struct task *ghost;
+
+    /*! Linked list of the tasks computing this cell's star density. */
+    struct link *density;
+
+    /*! Linked list of the tasks computing this cell's star feedback. */
+    struct link *feedback;
+
     /*! Max smoothing length in this cell. */
     double h_max;
 
@@ -627,6 +636,9 @@ struct cell {
     /*! Maximum part movement in this cell since last construction. */
     float dx_max_part;
 
+    /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+    float dx_max_part_old;
+
     /*! Maximum end of (integer) time step in this cell for black tasks. */
     integertime_t ti_end_min;
 
@@ -969,6 +981,43 @@ cell_can_recurse_in_self_stars_task(const struct cell *c) {
          (kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
 }
 
+/**
+ * @brief Can a sub-pair star task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param ci The #cell with black holes.
+ * @param cj The #cell with hydro parts.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_pair_black_holes_task(const struct cell *ci,
+                                          const struct cell *cj) {
+
+  /* Is the cell split ? */
+  /* If so, is the cut-off radius plus the max distance the parts have moved */
+  /* smaller than the sub-cell sizes ? */
+  /* Note: We use the _old values as these might have been updated by a drift */
+  return ci->split && cj->split &&
+         ((kernel_gamma * ci->black_holes.h_max_old +
+           ci->black_holes.dx_max_part_old) < 0.5f * ci->dmin) &&
+         ((kernel_gamma * cj->hydro.h_max_old + cj->hydro.dx_max_part_old) <
+          0.5f * cj->dmin);
+}
+
+/**
+ * @brief Can a sub-self black_holes task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param c The #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_self_black_holes_task(const struct cell *c) {
+
+  /* Is the cell split and not smaller than the smoothing length? */
+  return c->split &&
+         (kernel_gamma * c->black_holes.h_max_old < 0.5f * c->dmin) &&
+         (kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
+}
+
 /**
  * @brief Can a pair hydro task associated with a cell be split into smaller
  * sub-tasks.

src/runner.c

@@ -41,6 +41,7 @@
 #include "active.h"
 #include "approx_math.h"
 #include "atomic.h"
+#include "black_holes.h"
 #include "cell.h"
 #include "chemistry.h"
 #include "const.h"
@@ -126,6 +127,20 @@
 #undef FUNCTION_TASK_LOOP
 #undef FUNCTION
 
+/* Import the black hole density loop functions. */
+#define FUNCTION density
+#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
+#include "runner_doiact_black_holes.h"
+#undef FUNCTION_TASK_LOOP
+#undef FUNCTION
+
+/* Import the black hole feedback loop functions. */
+#define FUNCTION feedback
+#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
+#include "runner_doiact_black_holes.h"
+#undef FUNCTION_TASK_LOOP
+#undef FUNCTION
+
 /**
  * @brief Intermediate task after the density to check that the smoothing
  * lengths are correct.
@@ -3532,6 +3547,10 @@ void *runner_main(void *data) {
             runner_doself_branch_stars_density(r, ci);
           else if (t->subtype == task_subtype_stars_feedback)
             runner_doself_branch_stars_feedback(r, ci);
+          else if (t->subtype == task_subtype_bh_density)
+            runner_doself_branch_bh_density(r, ci);
+          else if (t->subtype == task_subtype_bh_feedback)
+            runner_doself_branch_bh_feedback(r, ci);
           else
             error("Unknown/invalid task subtype (%d).", t->subtype);
           break;
@@ -3553,6 +3572,10 @@ void *runner_main(void *data) {
             runner_dopair_branch_stars_density(r, ci, cj);
           else if (t->subtype == task_subtype_stars_feedback)
             runner_dopair_branch_stars_feedback(r, ci, cj);
+          else if (t->subtype == task_subtype_bh_density)
+            runner_dopair_branch_bh_density(r, ci, cj);
+          else if (t->subtype == task_subtype_bh_feedback)
+            runner_dopair_branch_bh_feedback(r, ci, cj);
           else
             error("Unknown/invalid task subtype (%d).", t->subtype);
           break;
@@ -3572,6 +3595,10 @@ void *runner_main(void *data) {
             runner_dosub_self_stars_density(r, ci, 1);
           else if (t->subtype == task_subtype_stars_feedback)
             runner_dosub_self_stars_feedback(r, ci, 1);
+          else if (t->subtype == task_subtype_bh_density)
+            runner_dosub_self_bh_density(r, ci, 1);
+          else if (t->subtype == task_subtype_bh_feedback)
+            runner_dosub_self_bh_feedback(r, ci, 1);
           else
             error("Unknown/invalid task subtype (%d).", t->subtype);
           break;
@@ -3591,6 +3618,10 @@ void *runner_main(void *data) {
             runner_dosub_pair_stars_density(r, ci, cj, 1);
           else if (t->subtype == task_subtype_stars_feedback)
             runner_dosub_pair_stars_feedback(r, ci, cj, 1);
+          else if (t->subtype == task_subtype_bh_density)
+            runner_dosub_pair_bh_density(r, ci, cj, 1);
+          else if (t->subtype == task_subtype_bh_feedback)
+            runner_dosub_pair_bh_feedback(r, ci, cj, 1);
           else
             error("Unknown/invalid task subtype (%d).", t->subtype);
           break;

src/space.c

@@ -235,6 +235,9 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->stars.ghost = NULL;
     c->stars.density = NULL;
     c->stars.feedback = NULL;
+    c->black_holes.ghost = NULL;
+    c->black_holes.density = NULL;
+    c->black_holes.feedback = NULL;
     c->kick1 = NULL;
     c->kick2 = NULL;
     c->timestep = NULL;

src/task.c

@@ -92,11 +92,11 @@ const char *taskID_names[task_type_count] = {"none",
 
 /* Sub-task type names. */
 const char *subtaskID_names[task_subtype_count] = {
-    "none",          "density",       "gradient",      "force",
-    "limiter",       "grav",          "external_grav", "tend_part",
-    "tend_gpart",    "tend_spart",    "tend_bpart",    "xv",
-    "rho",           "gpart",         "multipole",     "spart",
-    "stars_density", "stars_feedback"};
+    "none",          "density",        "gradient",      "force",
+    "limiter",       "grav",           "external_grav", "tend_part",
+    "tend_gpart",    "tend_spart",     "tend_bpart",    "xv",
+    "rho",           "gpart",          "multipole",     "spart",
+    "stars_density", "stars_feedback", "bh_density",    "bh_feedback"};
 
 #ifdef WITH_MPI
 /* MPI communicators for the subtypes. */

src/task.h

@@ -109,6 +109,8 @@ enum task_subtypes {
   task_subtype_spart,
   task_subtype_stars_density,
   task_subtype_stars_feedback,
+  task_subtype_bh_density,
+  task_subtype_bh_feedback,
   task_subtype_count
 } __attribute__((packed));