From 03c8026b099827cfc76ef36450f83dece068cfa7 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Thu, 2 May 2019 10:39:36 +0200 Subject: [PATCH] Read the properties of the EAGLE AGN model from the parameter file. --- examples/EAGLE_low_z/EAGLE_6/eagle_6.yml | 10 +- examples/main.c | 26 +++- src/Makefile.am | 6 +- .../Default/black_holes_properties.h | 100 ++++++++++++++ src/black_holes/EAGLE/black_holes.h | 19 ++- .../EAGLE/black_holes_properties.h | 130 ++++++++++++++++++ src/black_holes_properties.h | 34 +++++ src/engine.c | 3 + src/engine.h | 5 + src/runner.c | 14 +- src/swift.h | 1 + 11 files changed, 330 insertions(+), 18 deletions(-) create mode 100644 src/black_holes/Default/black_holes_properties.h create mode 100644 src/black_holes/EAGLE/black_holes_properties.h create mode 100644 src/black_holes_properties.h diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml index 3632fcd82..e560edbd7 100644 --- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml +++ b/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 diff --git a/examples/main.c b/examples/main.c index e3002cc7c..deb6a8544 100644 --- a/examples/main.c +++ b/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); diff --git a/src/Makefile.am b/src/Makefile.am index 8d07cc04c..595b88a38 100644 --- a/src/Makefile.am +++ b/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 diff --git a/src/black_holes/Default/black_holes_properties.h b/src/black_holes/Default/black_holes_properties.h new file mode 100644 index 000000000..8213bd214 --- /dev/null +++ b/src/black_holes/Default/black_holes_properties.h @@ -0,0 +1,100 @@ +/******************************************************************************* + * 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 . + * + ******************************************************************************/ +#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 */ diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h index e3e2d02cd..aaeed5783 100644 --- a/src/black_holes/EAGLE/black_holes.h +++ b/src/black_holes/EAGLE/black_holes.h @@ -19,13 +19,17 @@ #ifndef SWIFT_EAGLE_BLACK_HOLES_H #define SWIFT_EAGLE_BLACK_HOLES_H -#include +/* 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 + /** * @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); diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h new file mode 100644 index 000000000..a46756c56 --- /dev/null +++ b/src/black_holes/EAGLE/black_holes_properties.h @@ -0,0 +1,130 @@ +/******************************************************************************* + * 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 . + * + ******************************************************************************/ +#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 */ diff --git a/src/black_holes_properties.h b/src/black_holes_properties.h new file mode 100644 index 000000000..3226b5391 --- /dev/null +++ b/src/black_holes_properties.h @@ -0,0 +1,34 @@ +/******************************************************************************* + * 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 . + * + ******************************************************************************/ +#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 */ diff --git a/src/engine.c b/src/engine.c index 7f2a53205..d1c591c95 100644 --- a/src/engine.c +++ b/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; diff --git a/src/engine.h b/src/engine.h index 3afa8c221..1a653fdfd 100644 --- a/src/engine.h +++ b/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, diff --git a/src/runner.c b/src/runner.c index 7385f0c99..c0939bb0b 100644 --- a/src/runner.c +++ b/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); diff --git a/src/swift.h b/src/swift.h index cd7159e63..562a1fd5f 100644 --- a/src/swift.h +++ b/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" -- 2.26.2