/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2020 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_COLIBRE_BLACK_HOLES_PROPERTIES_H #define SWIFT_COLIBRE_BLACK_HOLES_PROPERTIES_H #include "chemistry.h" #include "hydro_properties.h" #include "string.h" enum AGN_feedback_models { AGN_isotropic_model, /*< Isotropic model of AGN feedback */ AGN_minimum_distance_model /*< Minimum-distance model of AGN feedback */ }; enum BH_merger_thresholds { BH_mergers_circular_velocity, /*< v_circ at separation, as in EAGLE */ BH_mergers_escape_velocity, /*< v_esc at separation */ BH_mergers_dynamical_escape_velocity /*< combined v_esc_dyn at separation */ }; /** * @brief Properties of black holes and AGN feedback in the EAGEL model. */ 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; /* ----- Initialisation properties ------ */ /*! Mass of a BH seed at creation time */ float subgrid_seed_mass; /*! Should we use the subgrid mass specified in ICs? */ int use_subgrid_mass_from_ics; /*! Should we enforce positive subgrid masses initially? */ int with_subgrid_mass_check; /* ----- Properties of the accretion model ------ */ /*! Calculate Bondi accretion rate based on subgrid properties? */ int use_subgrid_gas_properties; /*! Calculate Bondi accretion rate for individual neighbours? */ int use_multi_phase_bondi; /*! Switch between Bondi [0] or Krumholz [1] accretion rates */ int use_krumholz; /*! In Krumholz mode, should we include the vorticity term? */ int with_krumholz_vorticity; /*! Are we applying the angular-momentum-based multiplicative term from * Rosas-Guevara et al. (2015)? */ int with_angmom_limiter; /*! Normalisation of the viscuous angular momentum accretion reduction */ float alpha_visc; /*! Radiative efficiency of the black holes. */ float epsilon_r; /*! Maximal fraction of the Eddington rate allowed. */ float f_Edd; /*! Eddington fraction threshold for recording */ float f_Edd_recording; /*! Switch for the Booth, Schaye 2009 model */ int with_boost_factor; /*! Lowest value of the boost of the Booth, Schaye 2009 model */ float boost_alpha; /*! Power law slope for the boost of the Booth, Schaye 2009 model */ float boost_beta; /*! Normalisation density (internal units) for the boost of the Booth, Schaye * 2009 model */ double boost_n_h_star; /*! Switch for nibbling mode */ int use_nibbling; /*! Minimum gas particle mass in nibbling mode */ float min_gas_mass_for_nibbling; /* ---- Properties of the feedback model ------- */ /*! AGN feedback model: isotropic or minimum distance */ enum AGN_feedback_models feedback_model; /*! Is the AGN feedback model deterministic or stochastic? */ int AGN_deterministic; /*! Feedback coupling efficiency of the black holes. */ float epsilon_f; /*! Pivot (and also maximum possible) temperature in the variable AGN dT * feedback (Kelvin) */ float AGN_delta_T_pivot; /*! BH mass at which the pivot temperature is reached. The scaling with * the BH mass is linear, dT_AGN \propto dTAGN_pivot * (M_BH/M_dTAGN) */ float AGN_delta_T_Mass; /*! Minimum heating temperature in the variable AGN dT model (Kelvin) */ float AGN_delta_T_floor; /*! Number of gas neighbours to heat in a feedback event */ float num_ngbs_to_heat; /* ---- Properties of the repositioning model --- */ /*! Maximal mass of BH to reposition */ float max_reposition_mass; /*! Maximal distance to reposition, in units of softening length */ float max_reposition_distance_ratio; /*! Switch to enable a relative velocity limit for particles to which the * black holes can reposition */ int with_reposition_velocity_threshold; /*! Maximal velocity offset of particles to which the black hole can * reposition, in units of the ambient sound speed of the black hole */ float max_reposition_velocity_ratio; /*! Minimum value of the velocity repositioning threshold */ float min_reposition_velocity_threshold; /*! Switch to enable repositioning at fixed (maximum) speed */ int set_reposition_speed; /*! Normalisation factor for repositioning velocity */ float reposition_coefficient_upsilon; /*! Reference black hole mass for repositioning scaling */ float reposition_reference_mass; /*! Repositioning velocity scaling with black hole mass */ float reposition_exponent_mass; /*! Reference gas density for repositioning scaling */ float reposition_reference_n_gas; /*! Repositioning velocity scaling with gas density */ float reposition_exponent_n_gas; /*! Correct potential of BH? */ int correct_bh_potential_for_repositioning; /* ---- Properties of the merger model ---------- */ /*! Mass ratio above which a merger is considered 'minor' */ float minor_merger_threshold; /*! Mass ratio above which a merger is considered 'major' */ float major_merger_threshold; /*! Type of merger threshold */ enum BH_merger_thresholds merger_threshold_type; /*! Maximal distance over which BHs merge, in units of softening length */ float max_merging_distance_ratio; /* ---- Common conversion factors --------------- */ /*! Conversion factor from temperature to internal energy */ float temp_to_u_factor; /*! Conversion factor from physical density to n_H [cgs] */ float rho_to_n_cgs; /* ---- Black hole time-step properties ---------- */ /*! -- Minimum allowed time-step of BH in internal units */ float time_step_min; }; /** * @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)); /* Initialisation properties ---------------------------- */ bp->subgrid_seed_mass = parser_get_param_float(params, "COLIBREAGN:subgrid_seed_mass_Msun"); /* Convert to internal units */ bp->subgrid_seed_mass *= phys_const->const_solar_mass; bp->use_subgrid_mass_from_ics = parser_get_opt_param_int( params, "COLIBREAGN:use_subgrid_mass_from_ics", 1); if (bp->use_subgrid_mass_from_ics) bp->with_subgrid_mass_check = parser_get_opt_param_int( params, "COLIBREAGN:with_subgrid_mass_check", 1); /* Accretion parameters ---------------------------------- */ bp->use_subgrid_gas_properties = parser_get_param_int(params, "COLIBREAGN:use_subgrid_gas_properties"); bp->use_multi_phase_bondi = parser_get_param_int(params, "COLIBREAGN:use_multi_phase_bondi"); if (!bp->use_multi_phase_bondi) { bp->use_krumholz = parser_get_param_int(params, "COLIBREAGN:use_krumholz"); bp->with_krumholz_vorticity = parser_get_param_int(params, "COLIBREAGN:with_krumholz_vorticity"); } bp->with_angmom_limiter = parser_get_param_int(params, "COLIBREAGN:with_angmom_limiter"); if (bp->with_angmom_limiter) bp->alpha_visc = parser_get_param_float(params, "COLIBREAGN:viscous_alpha"); bp->epsilon_r = parser_get_param_float(params, "COLIBREAGN:radiative_efficiency"); if (bp->epsilon_r > 1) error("COLIBREAGN:radiative_efficiency must be <= 1, not %f.", bp->epsilon_r); bp->f_Edd = parser_get_param_float(params, "COLIBREAGN:max_eddington_fraction"); bp->f_Edd_recording = parser_get_param_float( params, "COLIBREAGN:eddington_fraction_for_recording"); /* Booth Schaye (2009) Parameters */ bp->with_boost_factor = parser_get_param_int(params, "COLIBREAGN:with_boost_factor"); if (bp->with_boost_factor) { bp->boost_alpha = parser_get_param_float(params, "COLIBREAGN:boost_alpha"); bp->boost_beta = parser_get_param_float(params, "COLIBREAGN:boost_beta"); /* Load the density in cgs and convert to internal units */ bp->boost_n_h_star = parser_get_param_float(params, "COLIBREAGN:boost_n_h_star_cm3") / units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY); } bp->use_nibbling = parser_get_param_int(params, "COLIBREAGN:use_nibbling"); if (bp->use_nibbling) { bp->min_gas_mass_for_nibbling = parser_get_param_float( params, "COLIBREAGN:min_gas_mass_for_nibbling_Msun"); bp->min_gas_mass_for_nibbling *= phys_const->const_solar_mass; } if ((bp->min_gas_mass_for_nibbling < 1e-5 * bp->subgrid_seed_mass) || (bp->min_gas_mass_for_nibbling > 1e5 * bp->subgrid_seed_mass)) { error( "The BH seeding mass and minimal gas mass for nibbling differ by more " "than 10^5. That is probably indicating a typo in the parameter file."); } /* Feedback parameters ---------------------------------- */ char temp[40]; parser_get_param_string(params, "COLIBREAGN:AGN_feedback_model", temp); if (strcmp(temp, "Isotropic") == 0) bp->feedback_model = AGN_isotropic_model; else if (strcmp(temp, "MinimumDistance") == 0) bp->feedback_model = AGN_minimum_distance_model; else error( "The AGN feedback model must be either MinimumDistance or Isotropic," " not %s", temp); bp->AGN_deterministic = parser_get_opt_param_int( params, "COLIBREAGN:AGN_use_deterministic_feedback", 1); bp->epsilon_f = parser_get_param_float(params, "COLIBREAGN:coupling_efficiency"); bp->AGN_delta_T_pivot = parser_get_param_float(params, "COLIBREAGN:AGN_delta_T_K"); /* Check that it makes sense. */ if (bp->AGN_delta_T_pivot <= 0.f) { error("The AGN heating temperature delta T must be > 0 K, not %.5e K.", bp->AGN_delta_T_pivot); } /* Get the floor temperature in the variable AGN dT model */ bp->AGN_delta_T_floor = parser_get_opt_param_float( params, "COLIBREAGN:AGN_delta_T_floor_K", bp->AGN_delta_T_pivot); bp->AGN_delta_T_pivot /= units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE); bp->AGN_delta_T_floor /= units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE); bp->num_ngbs_to_heat = parser_get_param_float(params, "COLIBREAGN:AGN_num_ngb_to_heat"); bp->AGN_delta_T_Mass = parser_get_param_float(params, "COLIBREAGN:delta_T_Mass_Msun"); /* Convert to internal units */ bp->AGN_delta_T_Mass *= phys_const->const_solar_mass; /* Reposition parameters --------------------------------- */ bp->max_reposition_mass = parser_get_param_float(params, "COLIBREAGN:max_reposition_mass_Msun"); /* Convert to internal units */ bp->max_reposition_mass *= phys_const->const_solar_mass; bp->max_reposition_distance_ratio = parser_get_param_float( params, "COLIBREAGN:max_reposition_distance_ratio"); bp->with_reposition_velocity_threshold = parser_get_param_int( params, "COLIBREAGN:with_reposition_velocity_threshold"); if (bp->with_reposition_velocity_threshold) { bp->max_reposition_velocity_ratio = parser_get_param_float( params, "COLIBREAGN:max_reposition_velocity_ratio"); /* Prevent nonsensical input */ if (bp->max_reposition_velocity_ratio <= 0) error("max_reposition_velocity_ratio must be positive, not %f.", bp->max_reposition_velocity_ratio); bp->min_reposition_velocity_threshold = parser_get_param_float( params, "COLIBREAGN:min_reposition_velocity_threshold_km_p_s"); /* Convert from km/s to internal units */ bp->min_reposition_velocity_threshold *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); } bp->set_reposition_speed = parser_get_param_int(params, "COLIBREAGN:set_reposition_speed"); if (bp->set_reposition_speed) { bp->reposition_coefficient_upsilon = parser_get_param_float( params, "COLIBREAGN:reposition_coefficient_upsilon"); /* Prevent the user from making silly wishes */ if (bp->reposition_coefficient_upsilon <= 0) error( "reposition_coefficient_upsilon must be positive, not %f " "km/s/M_sun.", bp->reposition_coefficient_upsilon); /* Convert from km/s to internal units */ bp->reposition_coefficient_upsilon *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); /* Scaling parameters with BH mass and gas density */ bp->reposition_reference_mass = parser_get_param_float(params, "COLIBREAGN:reposition_reference_mass_Msun") * phys_const->const_solar_mass; bp->reposition_exponent_mass = parser_get_param_float(params, "COLIBREAGN:reposition_exponent_mass"); bp->reposition_reference_n_gas = parser_get_param_float( params, "COLIBREAGN:reposition_reference_density_H_p_cm3"); bp->reposition_exponent_n_gas = parser_get_param_float( params, "COLIBREAGN:reposition_exponent_density"); } bp->correct_bh_potential_for_repositioning = parser_get_param_int(params, "COLIBREAGN:with_potential_correction"); /* Merger parameters ------------------------------------- */ bp->minor_merger_threshold = parser_get_param_float(params, "COLIBREAGN:threshold_minor_merger"); bp->major_merger_threshold = parser_get_param_float(params, "COLIBREAGN:threshold_major_merger"); char temp2[40]; parser_get_param_string(params, "COLIBREAGN:merger_threshold_type", temp2); if (strcmp(temp2, "CircularVelocity") == 0) bp->merger_threshold_type = BH_mergers_circular_velocity; else if (strcmp(temp2, "EscapeVelocity") == 0) bp->merger_threshold_type = BH_mergers_escape_velocity; else if (strcmp(temp2, "DynamicalEscapeVelocity") == 0) bp->merger_threshold_type = BH_mergers_dynamical_escape_velocity; else error( "The BH merger model must be either CircularVelocity, EscapeVelocity, " "or DynamicalEscapeVelocity, not %s", temp2); bp->max_merging_distance_ratio = parser_get_param_float(params, "COLIBREAGN:merger_max_distance_ratio"); /* Common conversion factors ----------------------------- */ /* Calculate temperature to internal energy conversion factor (all internal * units) */ const double k_B = phys_const->const_boltzmann_k; const double m_p = phys_const->const_proton_mass; const double mu = hydro_props->mu_ionised; bp->temp_to_u_factor = k_B / (mu * hydro_gamma_minus_one * m_p); /* Calculate conversion factor from rho to n_H. * Note this assumes primoridal abundance */ const double X_H = hydro_props->hydrogen_mass_fraction; bp->rho_to_n_cgs = (X_H / m_p) * units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY); /* ---- Black hole time-step properties ------------------ */ const double yr_in_cgs = 365.25 * 24. * 3600.; bp->time_step_min = parser_get_param_float(params, "COLIBREAGN:minimum_timestep_yr") * yr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME); } /** * @brief Write a black_holes_props struct to the given FILE as a stream of * bytes. * * @param props the black hole properties struct * @param stream the file stream */ INLINE static void black_holes_struct_dump( const struct black_holes_props *props, FILE *stream) { restart_write_blocks((void *)props, sizeof(struct black_holes_props), 1, stream, "black_holes props", "black holes props"); } /** * @brief Restore a black_holes_props struct from the given FILE as a stream of * bytes. * * @param props the black hole properties struct * @param stream the file stream */ INLINE static void black_holes_struct_restore( const struct black_holes_props *props, FILE *stream) { restart_read_blocks((void *)props, sizeof(struct black_holes_props), 1, stream, NULL, "black holes props"); } #endif /* SWIFT_COLIBRE_BLACK_HOLES_PROPERTIES_H */