Read the properties of the EAGLE AGN model from the parameter file.

examples/EAGLE_low_z/EAGLE_6/eagle_6.yml

@@ -68,7 +68,8 @@ InitialConditions:
   periodic:   1
   cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
-
+  black_holes_smoothing_length: 0.1
+  
 EAGLEChemistry: 	     # Solar abundances
   init_abundance_metal:      0.014
   init_abundance_Hydrogen:   0.70649785
@@ -151,3 +152,10 @@ EAGLEFeedback:
   SNII_yield_factor_Magnesium:      2.0             # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel.
   SNII_yield_factor_Silicon:        1.0             # (Optional) Correction factor to apply to the Silicon yield from the SNII channel.
   SNII_yield_factor_Iron:           0.5             # (Optional) Correction factor to apply to the Iron yield from the SNII channel.
+
+# EAGLE AGN model
+EAGLEAGN:
+  max_eddington_fraction:           1.         # Maximal allowed accretion rate in units of the Eddington rate.
+  radiative_efficiency:             0.1        # Fraction of the accreted mass that gets radiated.
+  coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
+  AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin. 
\ No newline at end of file

examples/main.c

@@ -98,6 +98,7 @@ int main(int argc, char *argv[]) {
   struct stars_props stars_properties;
   struct feedback_props feedback_properties;
   struct entropy_floor_properties entropy_floor;
+  struct black_holes_props black_holes_properties;
   struct part *parts = NULL;
   struct phys_const prog_const;
   struct space s;
@@ -389,6 +390,16 @@ int main(int argc, char *argv[]) {
     return 1;
   }
 
+  if (!with_hydro && with_black_holes) {
+    if (myrank == 0) {
+      argparse_usage(&argparse);
+      printf(
+          "\nError: Cannot process black holes without gas, --hydro must be "
+          "chosen.\n");
+    }
+    return 1;
+  }
+
 /* Let's pin the main thread, now we know if affinity will be used. */
 #if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
   if (with_aff &&
@@ -755,6 +766,16 @@ int main(int argc, char *argv[]) {
     } else
       bzero(&feedback_properties, sizeof(struct feedback_props));
 
+    /* Initialise the black holes properties */
+    if (with_black_holes) {
+#ifdef BLACK_HOLES_NONE
+      error("ERROR: Running with black_holes but compiled without it.");
+#endif
+      black_holes_props_init(&black_holes_properties, &prog_const, &us, params,
+                             &hydro_properties, &cosmo);
+    } else
+      bzero(&black_holes_properties, sizeof(struct black_holes_props));
+
     /* Initialise the gravity properties */
     if (with_self_gravity)
       gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology,
@@ -1005,8 +1026,9 @@ int main(int argc, char *argv[]) {
     engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
                 engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
                 &hydro_properties, &entropy_floor, &gravity_properties,
-                &stars_properties, &feedback_properties, &mesh, &potential,
-                &cooling_func, &starform, &chemistry);
+                &stars_properties, &black_holes_properties,
+                &feedback_properties, &mesh, &potential, &cooling_func,
+                &starform, &chemistry);
     engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
                   talking, restart_file);
 

src/Makefile.am

@@ -53,7 +53,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     star_formation_struct.h star_formation.h star_formation_iact.h \
     star_formation_logger.h star_formation_logger_struct.h \
     velociraptor_struct.h velociraptor_io.h random.h memuse.h black_holes.h black_holes_io.h \
-    feedback.h feedback_struct.h feedback_properties.h
+    black_holes_properties.h feedback.h feedback_struct.h feedback_properties.h
 
 # source files for EAGLE cooling
 EAGLE_COOLING_SOURCES =
@@ -197,8 +197,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.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_iact.h \
+                 black_holes/Default/black_holes_properties.h \
                  black_holes/EAGLE/black_holes.h black_holes/EAGLE/black_holes_io.h \
-		 black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h
+		 black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h \
+                 black_holes/EAGLE/black_holes_properties.h 
 
 
 # Sources and flags for regular library

src/black_holes/Default/black_holes_properties.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_BLACK_HOLES_PROPERTIES_H
+#define SWIFT_DEFAULT_BLACK_HOLES_PROPERTIES_H
+
+#include "chemistry.h"
+#include "hydro_properties.h"
+
+/**
+ * @brief Properties of the black hole scheme.
+ *
+ * In this default scheme, we only have the properties
+ * required by the neighbour search.
+ */
+struct black_holes_props {
+
+  /*! Resolution parameter */
+  float eta_neighbours;
+
+  /*! Target weightd number of neighbours (for info only)*/
+  float target_neighbours;
+
+  /*! Smoothing length tolerance */
+  float h_tolerance;
+
+  /*! Tolerance on neighbour number  (for info only)*/
+  float delta_neighbours;
+
+  /*! Maximal number of iterations to converge h */
+  int max_smoothing_iterations;
+
+  /*! Maximal change of h over one time-step */
+  float log_max_h_change;
+};
+
+/**
+ * @brief Initialise the black hole properties from the parameter file.
+ *
+ * We read the defaults from the hydro properties.
+ *
+ * @param bp The #black_holes_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ * @param hydro_props The already read-in properties of the hydro scheme.
+ * @param cosmo The cosmological model.
+ */
+void black_holes_props_init(struct black_holes_props *bp,
+                            const struct phys_const *phys_const,
+                            const struct unit_system *us,
+                            struct swift_params *params,
+                            const struct hydro_props *hydro_props,
+                            const struct cosmology *cosmo) {
+
+  /* Kernel properties */
+  bp->eta_neighbours = parser_get_opt_param_float(
+      params, "BlackHoles:resolution_eta", hydro_props->eta_neighbours);
+
+  /* Tolerance for the smoothing length Newton-Raphson scheme */
+  bp->h_tolerance = parser_get_opt_param_float(params, "BlackHoles:h_tolerance",
+                                               hydro_props->h_tolerance);
+
+  /* Get derived properties */
+  bp->target_neighbours = pow_dimension(bp->eta_neighbours) * kernel_norm;
+  const float delta_eta = bp->eta_neighbours * (1.f + bp->h_tolerance);
+  bp->delta_neighbours =
+      (pow_dimension(delta_eta) - pow_dimension(bp->eta_neighbours)) *
+      kernel_norm;
+
+  /* Number of iterations to converge h */
+  bp->max_smoothing_iterations =
+      parser_get_opt_param_int(params, "BlackHoles:max_ghost_iterations",
+                               hydro_props->max_smoothing_iterations);
+
+  /* Time integration properties */
+  const float max_volume_change =
+      parser_get_opt_param_float(params, "BlackHoles:max_volume_change", -1);
+  if (max_volume_change == -1)
+    bp->log_max_h_change = hydro_props->log_max_h_change;
+  else
+    bp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+}
+
+#endif /* SWIFT_DEFAULT_BLACK_HOLES_PROPERTIES_H */

src/black_holes/EAGLE/black_holes.h

@@ -19,13 +19,17 @@
 #ifndef SWIFT_EAGLE_BLACK_HOLES_H
 #define SWIFT_EAGLE_BLACK_HOLES_H
 
-#include <float.h>
+/* Local includes */
+#include "black_holes_properties.h"
 #include "cosmology.h"
 #include "dimension.h"
 #include "kernel_hydro.h"
 #include "minmax.h"
 #include "physical_constants.h"
 
+/* Standard includes */
+#include <float.h>
+
 /**
  * @brief Computes the gravity time-step of a given black hole particle.
  *
@@ -162,8 +166,9 @@ black_holes_bpart_has_no_neighbours(struct bpart* restrict bp,
  *
  */
 __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
-    struct bpart* restrict bp, const struct phys_const* constants,
-    const struct cosmology* cosmo, const double dt) {
+    struct bpart* restrict bp, const struct black_holes_props* props,
+    const struct phys_const* constants, const struct cosmology* cosmo,
+    const double dt) {
 
   /* Gather some physical constants (all in internal units) */
   const double G = constants->const_newton_G;
@@ -172,9 +177,9 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
   const double sigma_Thomson = constants->const_thomson_cross_section;
 
   /* Gather the parameters of the model */
-  const double f_Edd = 1.;
-  const double epsilon_r = 0.1;
-  const double epsilon_f = 0.15;
+  const double f_Edd = props->f_Edd;
+  const double epsilon_r = props->epsilon_r;
+  const double epsilon_f = props->epsilon_f;
 
   /* (Subgrid) mass of the BH (internal units) */
   const double BH_mass = bp->subgrid_mass;
@@ -207,7 +212,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
                             gas_rho_phys * denominator_inv * denominator_inv *
                             denominator_inv;
 
-  /* Compute the Eddington rate */
+  /* Compute the Eddington rate (internal units) */
   const double Eddington_rate =
       4. * M_PI * G * BH_mass * proton_mass / (epsilon_r * c * sigma_Thomson);
 

src/black_holes/EAGLE/black_holes_properties.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_EAGLE_BLACK_HOLES_PROPERTIES_H
+#define SWIFT_EAGLE_BLACK_HOLES_PROPERTIES_H
+
+#include "chemistry.h"
+#include "hydro_properties.h"
+
+struct black_holes_props {
+
+  /* ----- Basic neighbour search properties ------ */
+
+  /*! Resolution parameter */
+  float eta_neighbours;
+
+  /*! Target weightd number of neighbours (for info only)*/
+  float target_neighbours;
+
+  /*! Smoothing length tolerance */
+  float h_tolerance;
+
+  /*! Tolerance on neighbour number  (for info only)*/
+  float delta_neighbours;
+
+  /*! Maximal number of iterations to converge h */
+  int max_smoothing_iterations;
+
+  /*! Maximal change of h over one time-step */
+  float log_max_h_change;
+
+  /* ----- Properties of the accretion model ------ */
+
+  /*! Maximal fraction of the Eddington rate allowed. */
+  float f_Edd;
+
+  /*! Radiative efficiency of the black holes. */
+  float epsilon_r;
+
+  /*! Feedback coupling efficiency of the black holes. */
+  float epsilon_f;
+
+  /* ---- Properties of the feedback model ------- */
+
+  /*! Temperature increase induced by AGN feedback (Kelvin) */
+  float AGN_delta_T_desired;
+};
+
+/**
+ * @brief Initialise the black hole properties from the parameter file.
+ *
+ * For the basic black holes neighbour finding properties we use the
+ * defaults from the hydro scheme if the users did not provide specific
+ * values.
+ *
+ * @param bp The #black_holes_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ * @param hydro_props The already read-in properties of the hydro scheme.
+ * @param cosmo The cosmological model.
+ */
+INLINE static void black_holes_props_init(struct black_holes_props *bp,
+                                          const struct phys_const *phys_const,
+                                          const struct unit_system *us,
+                                          struct swift_params *params,
+                                          const struct hydro_props *hydro_props,
+                                          const struct cosmology *cosmo) {
+
+  /* Read in the basic neighbour search properties or default to the hydro
+     ones if the user did not provide any different values */
+
+  /* Kernel properties */
+  bp->eta_neighbours = parser_get_opt_param_float(
+      params, "BlackHoles:resolution_eta", hydro_props->eta_neighbours);
+
+  /* Tolerance for the smoothing length Newton-Raphson scheme */
+  bp->h_tolerance = parser_get_opt_param_float(params, "BlackHoles:h_tolerance",
+                                               hydro_props->h_tolerance);
+
+  /* Get derived properties */
+  bp->target_neighbours = pow_dimension(bp->eta_neighbours) * kernel_norm;
+  const float delta_eta = bp->eta_neighbours * (1.f + bp->h_tolerance);
+  bp->delta_neighbours =
+      (pow_dimension(delta_eta) - pow_dimension(bp->eta_neighbours)) *
+      kernel_norm;
+
+  /* Number of iterations to converge h */
+  bp->max_smoothing_iterations =
+      parser_get_opt_param_int(params, "BlackHoles:max_ghost_iterations",
+                               hydro_props->max_smoothing_iterations);
+
+  /* Time integration properties */
+  const float max_volume_change =
+      parser_get_opt_param_float(params, "BlackHoles:max_volume_change", -1);
+  if (max_volume_change == -1)
+    bp->log_max_h_change = hydro_props->log_max_h_change;
+  else
+    bp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+
+  /* Accretion parameters ---------------------------------- */
+
+  bp->f_Edd = parser_get_param_float(params, "EAGLEAGN:max_eddington_fraction");
+  bp->epsilon_r =
+      parser_get_param_float(params, "EAGLEAGN:radiative_efficiency");
+  bp->epsilon_f =
+      parser_get_param_float(params, "EAGLEAGN:coupling_efficiency");
+
+  /* Feedback parameters ---------------------------------- */
+
+  bp->AGN_delta_T_desired =
+      parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_K");
+}
+
+#endif /* SWIFT_EAGLE_BLACK_HOLES_PROPERTIES_H */

src/black_holes_properties.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_BLACK_HOLES_PROPERTIES_H
+#define SWIFT_BLACK_HOLES_PROPERTIES_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct black_holes model */
+#if defined(BLACK_HOLES_NONE)
+#include "./black_holes/Default/black_holes_properties.h"
+#elif defined(BLACK_HOLES_EAGLE)
+#include "./black_holes/EAGLE/black_holes_properties.h"
+#else
+#error "Invalid choice of black hole model"
+#endif
+
+#endif /* SWIFT_BLACK_HOLES_PROPERTIES_H */

src/engine.c

@@ -4858,6 +4858,7 @@ void engine_unpin(void) {
  * @param entropy_floor The #entropy_floor_properties for this run.
  * @param gravity The #gravity_props used for this run.
  * @param stars The #stars_props used for this run.
+ * @param black_holes The #black_holes_props used for this run.
  * @param feedback The #feedback_props used for this run.
  * @param mesh The #pm_mesh used for the long-range periodic forces.
  * @param potential The properties of the external potential.
@@ -4873,6 +4874,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
                  struct cosmology *cosmo, struct hydro_props *hydro,
                  const struct entropy_floor_properties *entropy_floor,
                  struct gravity_props *gravity, const struct stars_props *stars,
+                 const struct black_holes_props *black_holes,
                  const struct feedback_props *feedback, struct pm_mesh *mesh,
                  const struct external_potential *potential,
                  struct cooling_function_data *cooling_func,
@@ -4940,6 +4942,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->entropy_floor = entropy_floor;
   e->gravity_properties = gravity;
   e->stars_properties = stars;
+  e->black_holes_properties = black_holes;
   e->mesh = mesh;
   e->external_potential = potential;
   e->cooling_func = cooling_func;

src/engine.h

@@ -34,6 +34,7 @@
 
 /* Includes. */
 #include "barrier.h"
+#include "black_holes_properties.h"
 #include "chemistry_struct.h"
 #include "clocks.h"
 #include "collectgroup.h"
@@ -399,6 +400,9 @@ struct engine {
   /* Properties of the star model */
   const struct stars_props *stars_properties;
 
+  /* Properties of the black hole model */
+  const struct black_holes_props *black_holes_properties;
+
   /* Properties of the self-gravity scheme */
   struct gravity_props *gravity_properties;
 
@@ -473,6 +477,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
                  struct cosmology *cosmo, struct hydro_props *hydro,
                  const struct entropy_floor_properties *entropy_floor,
                  struct gravity_props *gravity, const struct stars_props *stars,
+                 const struct black_holes_props *black_holes,
                  const struct feedback_props *feedback, struct pm_mesh *mesh,
                  const struct external_potential *potential,
                  struct cooling_function_data *cooling_func,

src/runner.c

@@ -42,6 +42,7 @@
 #include "approx_math.h"
 #include "atomic.h"
 #include "black_holes.h"
+#include "black_holes_properties.h"
 #include "cell.h"
 #include "chemistry.h"
 #include "const.h"
@@ -515,9 +516,9 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
   const struct cosmology *cosmo = e->cosmology;
   const float black_holes_h_max = e->hydro_properties->h_max;
   const float black_holes_h_min = e->hydro_properties->h_min;
-  const float eps = e->hydro_properties->h_tolerance;
+  const float eps = e->black_holes_properties->h_tolerance;
   const float black_holes_eta_dim =
-      pow_dimension(e->hydro_properties->eta_neighbours);
+      pow_dimension(e->black_holes_properties->eta_neighbours);
   const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
   int redo = 0, bcount = 0;
 
@@ -633,8 +634,9 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
               ((bp->h <= black_holes_h_min) && (f > 0.f))) {
 
             /* Compute variables required for the feedback loop */
-            black_holes_prepare_feedback(bp, e->physical_constants,
-                                         e->cosmology, dt);
+            black_holes_prepare_feedback(bp, e->black_holes_properties,
+                                         e->physical_constants, e->cosmology,
+                                         dt);
 
             /* Reset quantities computed by the feedback loop */
             black_holes_reset_feedback(bp);
@@ -731,8 +733,8 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
         h_max = max(h_max, bp->h);
 
         /* Compute variables required for the feedback loop */
-        black_holes_prepare_feedback(bp, e->physical_constants, e->cosmology,
-                                     dt);
+        black_holes_prepare_feedback(bp, e->black_holes_properties,
+                                     e->physical_constants, e->cosmology, dt);
 
         /* Reset quantities computed by the feedback loop */
         black_holes_reset_feedback(bp);

src/swift.h

@@ -25,6 +25,7 @@
 /* Local headers. */
 #include "active.h"
 #include "atomic.h"
+#include "black_holes_properties.h"
 #include "cache.h"
 #include "cell.h"
 #include "chemistry.h"