diff --git a/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst b/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst index c276229ffce96bc222c6403a1e810eb3bddd2822..eef4286a1a52d0914cc5432873787f54ca8039ad 100644 --- a/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst +++ b/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst @@ -46,37 +46,39 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box. AGN_use_deterministic_feedback: 1 # Deterministic (1) or stochastic (0) AGN feedback model AGN_feedback_model: Isotropic # AGN feedback model (Isotropic or MinimumDistance) minimum_timestep_yr: 1000.0 # Minimum time-step of black-hole particles - max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate. + max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate. include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another. - alpha_acc: 0.1 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between thin and thick disk, as dot(m) = 0.2 * alpha^2. - seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1. - AGN_jet_velocity_model: BlackHoleMass # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported. - v_jet_km_p_s: 10000. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s. + alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between thin and thick disk, as dot(m) = 0.2 * alpha^2. + mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. + seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1. + AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported. + v_jet_km_p_s: 3160. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s. v_jet_cs_ratio: 10. # This sets the jet velocity to v_jet_cs_ratio times the sound speed of the hot gas of the parent halo the black hole is in. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass'. v_jet_BH_mass_scaling_reference_mass_Msun: 3.4e3 # The reference mass used in the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'. v_jet_BH_mass_scaling_slope: 0.65 # The slope of the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'. - v_jet_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is MassLoading. - v_jet_xi: 0.707 # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local' or 'SoundSpeed'. The appropriate values (to exactly obtain the formulas as derived) are 0.63 and 0.707 for the two, respectively. - v_jet_min_km_p_s: 500 # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading', 'Local' or 'SoundSpeed'. + v_jet_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is 'MassLoading'. + v_jet_xi: 1. # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local'. If a value of 1 is used, the formulas are used as derived, i.e. they are not rescaled. + v_jet_min_km_p_s: 500 # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'. + v_jet_max_km_p_s: 1e5 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'. + temperature_hot_gas_min_K: 1e4 # The minimum temperature (in Kelvin) of hot gas to include in the calculation of the hot gas sound speed around the BH. This is only used if 'AGN_jet_velocity_model' is 'Local'. opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests. N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector). AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity' eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1. fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests. jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. + fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector. + accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0]. radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH. include_ADIOS_suppression: 0 # Whether to suppress the accretion rate in the fully thick disc regime [1] (Eddington rate below 0.2alpha^2) by the amount expected to be taken away by isotropic kinetic disk winds. - ADIOS_R_in: 30. # If include_ADIOS_accr_suppression is set to 1, this parameter controls the inner radius within which winds are not important. - ADIOS_s: 0.4 # Slope of the accretion rate - radius relationship if include_ADIOS_accr_suppression is set to 1. turn_off_secondary_feedback: 1 # If set to 1, there will be only radiative (thermal) feedback in the thin disk mode, and only jets in the thick disk mode. jet_h_r_slope: 1. # The slope of the dependence of jet efficiency on aspect ratio of the subgrid accretion disk, H/R. Default value is 1, and another reasonable value is 0 (same jet efficiency for all disks). Reality could be anything in between. This parameter is only used if turn_off_secondary_feedback is set to 0. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0. - include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. - TD_SD_eps_r_threshold: 0.75 # Parameter controlling the transition from thin to slim disk. Accretion disk will be slim if radiative efficiency satisfies eps_slim < TD_SD_eps_r_threshold * eps_thin. This parameter is only used if include_slim_disk is set to 1. + include_slim_disk: 1 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. Most of these parameters should work well generally, and should not be changed except for tests. We will discuss only some of the more important ones. You can choose whether to have only the thick and thin disk (low and high BH accretion rates, respectively), or you can also include the slim disk at super-Eddington rates with ``include_slim_disk``. You can control what type of feedback you (do not) want with ``include_jets`` and ``turn_off_radiative_feedback``. If you choose to turn off jets, everything will be modeled as a thin disk (regardless of accretion rate), since jets go hand-in-hand with the thick and the slim disk. Similarly, if you turn off radiation, everything will be treated as a thick disk. diff --git a/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml b/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml index 501961c070c5e00f22e339c6d2961430942c1a8a..1e4987e51809807150234f11fdd58ecb1b35ea6b 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml @@ -233,4 +233,69 @@ EAGLEAGN: merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. - +# Spin and jet AGN model (Husko et al. 2022) +SPINJETAGN: + subgrid_seed_mass_Msun: 1.0e4 # Black hole subgrid mass at creation time in solar masses. + use_multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle? + use_subgrid_bondi: 0 # Compute Bondi rates using the subgrid extrapolation of the gas properties around the BH? + with_angmom_limiter: 0 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term? + viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term + with_boost_factor: 0 # Are we using the model from Booth & Schaye (2009)? + boost_alpha_only: 0 # If using the boost factor, are we using a constant boost only? + boost_alpha: 1. # Lowest value for the accretion effeciency for the Booth & Schaye 2009 accretion model. + boost_beta: 2. # Slope of the power law for the Booth & Schaye 2009 model, set beta to zero for constant alpha models. + boost_n_h_star_H_p_cm3: 0.1 # Normalization of the power law for the Booth & Schaye 2009 model in cgs (cm^-3). + with_fixed_T_near_EoS: 0 # Are we using a fixed temperature to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term? + fixed_T_above_EoS_dex: 0.3 # Distance above the entropy floor for which we use a fixed sound-speed + fixed_T_near_EoS_K: 8000 # Fixed temperature assumed to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term + radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated. + use_nibbling: 1 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]? + min_gas_mass_for_nibbling_Msun: 9e5 # Minimum mass for a gas particle to be nibbled from [M_Sun]. Only used if use_nibbling is 1. + max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate. + eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold. + coupling_efficiency: 0.1 # Fraction of the radiated energy that couples to the gas in feedback events. + AGN_feedback_model: MinimumDistance # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity + AGN_use_deterministic_feedback: 1 # Deterministic (reservoir) [1] or stochastic [0] AGN feedback? + use_variable_delta_T: 1 # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0]. + AGN_with_locally_adaptive_delta_T: 1 # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_norm: 3e8 # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_reference: 1e8 # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_exponent: 0.666667 # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1). + AGN_delta_T_crit_factor: 1.0 # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_background_factor: 0.0 # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_min: 1e7 # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_max: 3e9 # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event [K] (used if use_variable_delta_T is 0 or AGN_use_nheat_with_fixed_dT is 1 AND to initialise the BHs). + AGN_use_nheat_with_fixed_dT: 0 # Switch to use the constant AGN dT, rather than the adaptive one, for calculating the energy reservoir threshold. + AGN_use_adaptive_energy_reservoir_threshold: 0 # Switch to calculate an adaptive AGN energy reservoir threshold. + AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event (only used if AGN_use_adaptive_energy_reservoir_threshold is 0). + max_reposition_mass: 1e20 # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition). + max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length. + with_reposition_velocity_threshold: 0 # Should we only reposition to particles that move slowly w.r.t. the black hole? + max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_threshold is 1. + min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_threshold is 1. + set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum? + with_potential_correction: 1 # Should the BH's own contribution to the potential be removed from the neighbour's potentials when looking for repositioning targets. + threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major' + threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor' + merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity'). + merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. + minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. + include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. + turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. + alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. + TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. + N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. + jet_velocity_model: Constant # Which jet velocity model to use. + v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'. + opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests. + eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1. + fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests. + jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. + AGN_jet_feedback_model: MinimunDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity' + include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets. + turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates. + fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. + fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. + include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. \ No newline at end of file diff --git a/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml b/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml index 0ebd235b2572ba805b436069da715e951a5feb00..c0a8d256ed011b222ad848ad42ec97d636598e2f 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml @@ -233,4 +233,69 @@ EAGLEAGN: merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. - +# Spin and jet AGN model (Husko et al. 2022) +SPINJETAGN: + subgrid_seed_mass_Msun: 1.0e4 # Black hole subgrid mass at creation time in solar masses. + use_multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle? + use_subgrid_bondi: 0 # Compute Bondi rates using the subgrid extrapolation of the gas properties around the BH? + with_angmom_limiter: 0 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term? + viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term + with_boost_factor: 0 # Are we using the model from Booth & Schaye (2009)? + boost_alpha_only: 0 # If using the boost factor, are we using a constant boost only? + boost_alpha: 1. # Lowest value for the accretion effeciency for the Booth & Schaye 2009 accretion model. + boost_beta: 2. # Slope of the power law for the Booth & Schaye 2009 model, set beta to zero for constant alpha models. + boost_n_h_star_H_p_cm3: 0.1 # Normalization of the power law for the Booth & Schaye 2009 model in cgs (cm^-3). + with_fixed_T_near_EoS: 0 # Are we using a fixed temperature to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term? + fixed_T_above_EoS_dex: 0.3 # Distance above the entropy floor for which we use a fixed sound-speed + fixed_T_near_EoS_K: 8000 # Fixed temperature assumed to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term + radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated. + use_nibbling: 1 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]? + min_gas_mass_for_nibbling_Msun: 9e5 # Minimum mass for a gas particle to be nibbled from [M_Sun]. Only used if use_nibbling is 1. + max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate. + eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold. + coupling_efficiency: 0.1 # Fraction of the radiated energy that couples to the gas in feedback events. + AGN_feedback_model: MinimumDistance # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity + AGN_use_deterministic_feedback: 1 # Deterministic (reservoir) [1] or stochastic [0] AGN feedback? + use_variable_delta_T: 1 # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0]. + AGN_with_locally_adaptive_delta_T: 1 # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_norm: 3e8 # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_reference: 1e8 # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_exponent: 0.666667 # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1). + AGN_delta_T_crit_factor: 1.0 # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_background_factor: 0.0 # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_min: 1e7 # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_max: 3e9 # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event [K] (used if use_variable_delta_T is 0 or AGN_use_nheat_with_fixed_dT is 1 AND to initialise the BHs). + AGN_use_nheat_with_fixed_dT: 0 # Switch to use the constant AGN dT, rather than the adaptive one, for calculating the energy reservoir threshold. + AGN_use_adaptive_energy_reservoir_threshold: 0 # Switch to calculate an adaptive AGN energy reservoir threshold. + AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event (only used if AGN_use_adaptive_energy_reservoir_threshold is 0). + max_reposition_mass: 1e20 # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition). + max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length. + with_reposition_velocity_threshold: 0 # Should we only reposition to particles that move slowly w.r.t. the black hole? + max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_threshold is 1. + min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_threshold is 1. + set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum? + with_potential_correction: 1 # Should the BH's own contribution to the potential be removed from the neighbour's potentials when looking for repositioning targets. + threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major' + threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor' + merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity'). + merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. + minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. + include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. + turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. + alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. + TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. + N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. + jet_velocity_model: Constant # Which jet velocity model to use. + v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'. + opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests. + eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1. + fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests. + jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. + AGN_jet_feedback_model: MinimunDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity' + include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets. + turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates. + fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. + fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. + include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. diff --git a/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml b/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml index 742a9e50a4f70cd35fc9d54084c0501120962782..b29df144d0c3f5e98c918321dd64bfc7ef5a63f6 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml @@ -233,4 +233,69 @@ EAGLEAGN: merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. - +# Spin and jet AGN model (Husko et al. 2022) +SPINJETAGN: + subgrid_seed_mass_Msun: 1.0e4 # Black hole subgrid mass at creation time in solar masses. + use_multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle? + use_subgrid_bondi: 0 # Compute Bondi rates using the subgrid extrapolation of the gas properties around the BH? + with_angmom_limiter: 0 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term? + viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term + with_boost_factor: 0 # Are we using the model from Booth & Schaye (2009)? + boost_alpha_only: 0 # If using the boost factor, are we using a constant boost only? + boost_alpha: 1. # Lowest value for the accretion effeciency for the Booth & Schaye 2009 accretion model. + boost_beta: 2. # Slope of the power law for the Booth & Schaye 2009 model, set beta to zero for constant alpha models. + boost_n_h_star_H_p_cm3: 0.1 # Normalization of the power law for the Booth & Schaye 2009 model in cgs (cm^-3). + with_fixed_T_near_EoS: 0 # Are we using a fixed temperature to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term? + fixed_T_above_EoS_dex: 0.3 # Distance above the entropy floor for which we use a fixed sound-speed + fixed_T_near_EoS_K: 8000 # Fixed temperature assumed to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term + radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated. + use_nibbling: 1 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]? + min_gas_mass_for_nibbling_Msun: 9e5 # Minimum mass for a gas particle to be nibbled from [M_Sun]. Only used if use_nibbling is 1. + max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate. + eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold. + coupling_efficiency: 0.1 # Fraction of the radiated energy that couples to the gas in feedback events. + AGN_feedback_model: MinimumDistance # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity + AGN_use_deterministic_feedback: 1 # Deterministic (reservoir) [1] or stochastic [0] AGN feedback? + use_variable_delta_T: 1 # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0]. + AGN_with_locally_adaptive_delta_T: 1 # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_norm: 3e8 # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_reference: 1e8 # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_exponent: 0.666667 # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1). + AGN_delta_T_crit_factor: 1.0 # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_background_factor: 0.0 # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_min: 1e7 # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_max: 3e9 # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event [K] (used if use_variable_delta_T is 0 or AGN_use_nheat_with_fixed_dT is 1 AND to initialise the BHs). + AGN_use_nheat_with_fixed_dT: 0 # Switch to use the constant AGN dT, rather than the adaptive one, for calculating the energy reservoir threshold. + AGN_use_adaptive_energy_reservoir_threshold: 0 # Switch to calculate an adaptive AGN energy reservoir threshold. + AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event (only used if AGN_use_adaptive_energy_reservoir_threshold is 0). + max_reposition_mass: 1e20 # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition). + max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length. + with_reposition_velocity_threshold: 0 # Should we only reposition to particles that move slowly w.r.t. the black hole? + max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_threshold is 1. + min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_threshold is 1. + set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum? + with_potential_correction: 1 # Should the BH's own contribution to the potential be removed from the neighbour's potentials when looking for repositioning targets. + threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major' + threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor' + merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity'). + merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. + minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. + include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. + turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. + alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. + TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. + N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. + jet_velocity_model: Constant # Which jet velocity model to use. + v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'. + opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests. + eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1. + fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests. + jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. + AGN_jet_feedback_model: MinimunDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity' + include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets. + turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates. + fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. + fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. + include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. diff --git a/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml b/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml index b5b69a1f81151f79e338f6f9ea93a14d248c06e7..f029c7c0194d025f233a176c9edc3dd05d30cdcb 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml @@ -232,3 +232,70 @@ EAGLEAGN: merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity'). merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. + +# Spin and jet AGN model (Husko et al. 2022) +SPINJETAGN: + subgrid_seed_mass_Msun: 1.0e4 # Black hole subgrid mass at creation time in solar masses. + use_multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle? + use_subgrid_bondi: 0 # Compute Bondi rates using the subgrid extrapolation of the gas properties around the BH? + with_angmom_limiter: 0 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term? + viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term + with_boost_factor: 0 # Are we using the model from Booth & Schaye (2009)? + boost_alpha_only: 0 # If using the boost factor, are we using a constant boost only? + boost_alpha: 1. # Lowest value for the accretion effeciency for the Booth & Schaye 2009 accretion model. + boost_beta: 2. # Slope of the power law for the Booth & Schaye 2009 model, set beta to zero for constant alpha models. + boost_n_h_star_H_p_cm3: 0.1 # Normalization of the power law for the Booth & Schaye 2009 model in cgs (cm^-3). + with_fixed_T_near_EoS: 0 # Are we using a fixed temperature to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term? + fixed_T_above_EoS_dex: 0.3 # Distance above the entropy floor for which we use a fixed sound-speed + fixed_T_near_EoS_K: 8000 # Fixed temperature assumed to compute the sound-speed of gas on the entropy floor in the Bondy-Hoyle accretion term + radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated. + use_nibbling: 1 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]? + min_gas_mass_for_nibbling_Msun: 9e5 # Minimum mass for a gas particle to be nibbled from [M_Sun]. Only used if use_nibbling is 1. + max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate. + eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold. + coupling_efficiency: 0.1 # Fraction of the radiated energy that couples to the gas in feedback events. + AGN_feedback_model: MinimumDistance # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity + AGN_use_deterministic_feedback: 1 # Deterministic (reservoir) [1] or stochastic [0] AGN feedback? + use_variable_delta_T: 1 # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0]. + AGN_with_locally_adaptive_delta_T: 1 # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_norm: 3e8 # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_reference: 1e8 # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1). + AGN_delta_T_mass_exponent: 0.666667 # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1). + AGN_delta_T_crit_factor: 1.0 # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_background_factor: 0.0 # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1). + AGN_delta_T_min: 1e7 # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_max: 3e9 # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1). + AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event [K] (used if use_variable_delta_T is 0 or AGN_use_nheat_with_fixed_dT is 1 AND to initialise the BHs). + AGN_use_nheat_with_fixed_dT: 0 # Switch to use the constant AGN dT, rather than the adaptive one, for calculating the energy reservoir threshold. + AGN_use_adaptive_energy_reservoir_threshold: 0 # Switch to calculate an adaptive AGN energy reservoir threshold. + AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event (only used if AGN_use_adaptive_energy_reservoir_threshold is 0). + max_reposition_mass: 1e20 # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition). + max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length. + with_reposition_velocity_threshold: 0 # Should we only reposition to particles that move slowly w.r.t. the black hole? + max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_threshold is 1. + min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_threshold is 1. + set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum? + with_potential_correction: 1 # Should the BH's own contribution to the potential be removed from the neighbour's potentials when looking for repositioning targets. + threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major' + threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor' + merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity'). + merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length. + minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length. + include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. + turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. + alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. + TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. + N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. + jet_velocity_model: Constant # Which jet velocity model to use. + v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'. + opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests. + eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1. + fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests. + jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. + AGN_jet_feedback_model: MinimunDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity' + include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets. + turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates. + fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. + fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. + include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. \ No newline at end of file diff --git a/examples/IdealisedCluster/README b/examples/IdealisedCluster/README index 18689194da63f0ccba8bd627bae226b41d82a39c..17ebf3bd43a8654c0ca72d1038934b94a0e7cfc5 100644 --- a/examples/IdealisedCluster/README +++ b/examples/IdealisedCluster/README @@ -15,11 +15,14 @@ different resolution and different central temperatures. The code should be configured using a NFW profile and the desired subgrid models: -./configure --with-ext-potential=nfw --enable-fixed-boundary-particles=2 +./configure --with-subgrid=x --with-ext-potential=nfw +--enable-fixed-boundary-particles=2 +where x is EAGLE if the desired subgrid model is similar to that used in +Nobels et al. (2022) and SPIN_JET_EAGLE for the model used in Husko et al. +(2022), which has jet AGN feedback instead of thermal isotropic feedback. -The first argument sets the external potential to an NFW profile and the -second argument pins the supermassive BH of the simulation to the centre -of the halo. +The remaining arguments set the external potential to an NFW profile and pin +the supermassive BH of the simulation to the centre of the halo. These ICs are ideal for testing galaxy formation subgrid model implementation. diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index f79bb757ccd0409a3b4321a10f1f2a09bb4f04cb..520d3e2a62f1a06f9b8543b84df9ece49d2d1f86 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -742,34 +742,35 @@ SPINJETAGN: minimum_timestep_yr: 1000.0 # Minimum time-step of black-hole particles include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another. - alpha_acc: 0.1 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2. - seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1. - AGN_jet_velocity_model: BlackHoleMass # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported. - v_jet_km_p_s: 10000. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s. + alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2. + mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. + seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1. + AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported. + v_jet_km_p_s: 3160. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s. v_jet_cs_ratio: 10. # This sets the jet velocity to v_jet_cs_ratio times the sound speed of the hot gas of the parent halo the black hole is in. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass'. v_jet_BH_mass_scaling_reference_mass_Msun: 3.4e3 # The reference mass used in the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'. v_jet_BH_mass_scaling_slope: 0.65 # The slope of the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'. - v_jet_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is MassLoading. - v_jet_xi: 0.707 # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local' or 'SoundSpeed'. The appropriate values (to exactly obtain the formulas as derived) are 0.63 and 0.707 for the two, respectively. - v_jet_min_km_p_s: 500 # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading', 'Local' or 'SoundSpeed'. + v_jet_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is 'MassLoading'. + v_jet_xi: 1. # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local'. If a value of 1 is used, the formulas are used as derived, i.e. they are not rescaled. + v_jet_min_km_p_s: 500 # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'. + v_jet_max_km_p_s: 1e5 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'. + temperature_hot_gas_min_K: 1e4 # The minimum temperature (in Kelvin) of hot gas to include in the calculation of the hot gas sound speed around the BH. This is only used if 'AGN_jet_velocity_model' is 'Local'. opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests. N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector). AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity' eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1. - fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests. + fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. + fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector. + accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0]. radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH. - include_ADIOS_suppression: 0 # Whether to suppress the accretion rate in the fully thick disc regime [1] (Eddington rate below 0.2alpha^2) by the amount expected to be taken away by isotropic kinetic disk winds. - ADIOS_R_in: 30. # If include_ADIOS_accr_suppression is set to 1, this parameter controls the inner radius within which winds are not important. - ADIOS_s: 0.4 # Slope of the accretion rate - radius relationship if include_ADIOS_accr_suppression is set to 1. turn_off_secondary_feedback: 1 # If set to 1, there will be only radiative (thermal) feedback in the thin disk mode, and only jets in the thick disk mode. jet_h_r_slope: 1. # The slope of the dependence of jet efficiency on aspect ratio of the subgrid accretion disk, H/R. Default value is 1, and another reasonable value is 0 (same jet efficiency for all disks). Reality could be anything in between. This parameter is only used if turn_off_secondary_feedback is set to 0. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0. - include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. - TD_SD_eps_r_threshold: 0.75 # Parameter controlling the transition from thin to slim disk. Accretion disk will be slim if radiative efficiency satisfies eps_slim < TD_SD_eps_r_threshold * eps_thin. This parameter is only used if include_slim_disk is set to 1. + include_slim_disk: 1 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. # Parameters related to the neutrinos -------------------------------------------- diff --git a/src/black_holes/SPIN_JET/black_holes.h b/src/black_holes/SPIN_JET/black_holes.h index 143958a469fb6409464c922c31a1cc309175c442..b107b31fbbe35d1bcafd53e8aebc2d10177c0dc1 100644 --- a/src/black_holes/SPIN_JET/black_holes.h +++ b/src/black_holes/SPIN_JET/black_holes.h @@ -170,6 +170,9 @@ __attribute__((always_inline)) INLINE static void black_holes_first_init_bpart( bp->thermal_energy_by_mode[i] = 0.f; for (int i = 0; i < BH_accretion_modes_count; ++i) bp->jet_energy_by_mode[i] = 0.f; + bp->jet_direction[0] = bp->angular_momentum_direction[0]; + bp->jet_direction[1] = bp->angular_momentum_direction[1]; + bp->jet_direction[2] = bp->angular_momentum_direction[2]; black_holes_mark_bpart_as_not_swallowed(&bp->merger_data); } @@ -212,6 +215,8 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart( bp->accretion_rate = 0.f; /* Optionally accumulated ngb-by-ngb */ bp->f_visc = FLT_MAX; bp->mass_at_start_of_step = bp->mass; /* bp->mass may grow in nibbling mode */ + bp->rho_gas_hot = 0.f; + bp->sound_speed_gas_hot = 0.f; /* Reset the rays carried by this BH */ ray_init(bp->rays, spinjet_blackhole_number_of_rays); @@ -312,10 +317,13 @@ __attribute__((always_inline)) INLINE static void black_holes_end_density( bp->density.wcount *= h_inv_dim; bp->density.wcount_dh *= h_inv_dim_plus_one; bp->rho_gas *= h_inv_dim; + bp->rho_gas_hot *= h_inv_dim; const float rho_inv = 1.f / bp->rho_gas; /* For the following, we also have to undo the mass smoothing * (N.B.: bp->velocity_gas is in BH frame, in internal units). */ + if (bp->rho_gas_hot > 0.) + bp->sound_speed_gas_hot *= h_inv_dim / bp->rho_gas_hot; bp->sound_speed_gas *= h_inv_dim * rho_inv; bp->velocity_gas[0] *= h_inv_dim * rho_inv; bp->velocity_gas[1] *= h_inv_dim * rho_inv; @@ -825,22 +833,18 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback( bp->accretion_rate = accr_rate; bp->eddington_fraction = accr_rate / Eddington_rate; + /* If we are in the thick disc mode, suppress the accretion rate by the + * accretion efficiency. */ + if (bp->accretion_mode == BH_thick_disc) { + accr_rate *= props->accretion_efficiency; + bp->accretion_rate *= props->accretion_efficiency; + bp->eddington_fraction *= props->accretion_efficiency; + } + /* Define feedback-related quantities that we will update and need later on */ double luminosity = 0.; double jet_power = 0.; - /* Check whether we are including ADIOS winds in the thick disk */ - if ((bp->accretion_mode == BH_thick_disc) && - (props->include_ADIOS_suppression)) { - const double Bondi_R = G * BH_mass / gas_c_phys2; - const float ADIOS_suppression = powf( - Bondi_R / (props->ADIOS_R_in * R_gravitational(BH_mass, constants)), - props->ADIOS_s); - accr_rate = accr_rate / ADIOS_suppression; - bp->accretion_rate = accr_rate; - bp->eddington_fraction = accr_rate / Eddington_rate; - } - /* How much mass will be consumed over this time step? */ double delta_m_0 = bp->accretion_rate * dt; @@ -920,23 +924,12 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback( } /* Check if we are fixing the jet along the z-axis. */ - if (props->fix_jet_efficiency) { + if (props->fix_jet_direction) { bp->angular_momentum_direction[0] = 0.; bp->angular_momentum_direction[1] = 0.; bp->angular_momentum_direction[2] = 1; } - /* The amount of mass the BH is actually swallowing, including the - effects of efficiencies. */ - const double delta_m_real = - delta_m_0 * (1. - rad_efficiency(bp, props) - jet_efficiency(bp, props)); - - /* Increase the reservoir*/ - bp->jet_reservoir += - delta_m_0 * c * c * props->eps_f_jet * jet_efficiency(bp, props); - bp->energy_reservoir += - delta_m_0 * c * c * props->epsilon_f * rad_efficiency(bp, props); - float spin_final = -1.; /* Calculate the change in the BH spin */ if (bp->subgrid_mass > 0.) { @@ -967,8 +960,35 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback( spin_final = 0.01; } - /* Update the spin and mass. */ + /* Update the spin */ bp->spin = spin_final; + + /* Update other quantities */ + bp->accretion_disk_angle = dot_product; + bp->aspect_ratio = aspect_ratio(bp, constants, props); + + /* Update efficiencies given the new spin */ + bp->radiative_efficiency = rad_efficiency(bp, props); + bp->jet_efficiency = jet_efficiency(bp, props); + + /* Final jet power at the end of the step */ + jet_power = bp->jet_efficiency * accr_rate * c * c; + + /* Final luminosity at the end of the step */ + luminosity = bp->radiative_efficiency * accr_rate * c * c; + + /* The amount of mass the BH is actually swallowing, including the + effects of the updated efficiencies. */ + const double delta_m_real = + delta_m_0 * (1. - rad_efficiency(bp, props) - jet_efficiency(bp, props)); + + /* Increase the reservoir */ + bp->jet_reservoir += + delta_m_0 * c * c * props->eps_f_jet * jet_efficiency(bp, props); + bp->energy_reservoir += + delta_m_0 * c * c * props->epsilon_f * rad_efficiency(bp, props); + + /* Update the masses */ bp->subgrid_mass += delta_m_real; bp->total_accreted_mass += delta_m_real; @@ -983,26 +1003,6 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback( bp->subgrid_mass = props->subgrid_seed_mass; } - /* Update other quantities. */ - bp->accretion_disk_angle = dot_product; - bp->aspect_ratio = aspect_ratio(bp, constants, props); - if (props->fix_radiative_efficiency) { - bp->radiative_efficiency = props->radiative_efficiency; - } else { - bp->radiative_efficiency = rad_efficiency(bp, props); - } - if (props->fix_jet_efficiency) { - bp->jet_efficiency = props->jet_efficiency; - } else { - bp->jet_efficiency = jet_efficiency(bp, props); - } - - /* Final jet power at the end of the step */ - jet_power = bp->jet_efficiency * accr_rate * c * c; - - /* Final luminosity at the end of the step */ - luminosity = bp->radiative_efficiency * accr_rate * c * c; - /* Increase the subgrid angular momentum according to what we accreted * (already in physical units, a factors from velocity and radius cancel) */ bp->accreted_angular_momentum[0] += @@ -1214,33 +1214,40 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback( */ int N_unsuccessful_jet_injections = 0; - for (int i = 0; i < number_of_jet_injections; i++) { + /* Store all of this in the black hole for delivery onto the gas. */ + bp->to_distribute.AGN_delta_u_jet = delta_u_jet; + bp->to_distribute.AGN_number_of_jet_injections = number_of_jet_injections; + + /* Check for situations where the gas has been swallowed (not used in the + * nibbling scheme, or for situations where one of the BH hemispheres is + * empty. In this case we do not perform any jet kicks at all. */ + for (int i = 0; i < number_of_jet_injections / 2; i++) { /* If the gas particle that the BH wants to heat has just been swallowed - * by the same BH, increment the counter of unsuccessful injections. If - * the particle has not been swallowed by the BH, increase the energy that - * will later be subtracted from the BH's energy reservoir. - * Loop over jet both rays.*/ - - if (i % 2 == 0) { - if (bp->rays_jet[i].id_min_length != -1) { - Jet_energy_deposited += delta_u_jet * bp->rays_jet[i].mass; - } else { - N_unsuccessful_jet_injections++; - } + * by the same BH (or if there are no particles in either of the jet + * jet rays, i.e. on either side of the BH hemisphere), increment the + * counter of unsuccessful injections. If the particle has not been + * swallowed by the BH and both hemispheres are populated, the energy + * to be deposited will be subracted from the jet reservoir. */ + + if ((bp->rays_jet[i].id_min_length != -1) && + (bp->rays_jet_pos[i].id_min_length != -1)) { + + /* Neither hemisphere of the BH is empty, nor have any particles been + * swallowed, so we can safely deposity the jet energy. */ + Jet_energy_deposited += delta_u_jet * bp->rays_jet[i].mass; + Jet_energy_deposited += delta_u_jet * bp->rays_jet_pos[i].mass; + } else { - if (bp->rays_jet_pos[i].id_min_length != -1) { - Jet_energy_deposited += delta_u_jet * bp->rays_jet_pos[i].mass; - } else { - N_unsuccessful_jet_injections++; - } + + /* Track how many unsuccessful jet injections happened. */ + N_unsuccessful_jet_injections += 2; + + /* Reduce the number of jet injections to be handed out. */ + bp->to_distribute.AGN_number_of_jet_injections -= 2; } } - /* Store all of this in the black hole for delivery onto the gas. */ - bp->to_distribute.AGN_delta_u_jet = delta_u_jet; - bp->to_distribute.AGN_number_of_jet_injections = number_of_jet_injections; - /* Subtract the deposited energy from the BH jet energy reservoir. Note * that in the stochastic case, the resulting value might be negative. * This happens when (due to the probabilistic nature of the model) the @@ -1254,6 +1261,15 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback( const int N_successful_jet_injections = number_of_jet_injections - N_unsuccessful_jet_injections; + /* If we will be kicking something, choose a new direction for the jets. + * This is done by randomly generating a unit vector within a cone of a + * given opening angle around the BH spin vector. */ + if (N_successful_jet_injections > 0) { + random_direction_in_cone( + bp->id, ti_begin, random_number_BH_kick, props->opening_angle, + bp->angular_momentum_direction, bp->jet_direction); + } + /* Increase the number of jet energy injections the black hole has kicked * so far. */ bp->AGN_number_of_jet_injections += N_successful_jet_injections; @@ -1290,13 +1306,57 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback( * fraction */ decide_mode(bp, props); - /* Calculate a BH angular momentum evolution time step. Ths timestep is - chosen so that the BH spin changes by around 10% relative to the current - spin value */ + /* The accretion/feedback mode is now possibly different, so the feedback + efficiencies need to be updated. This is important for computing the + next time-step of the BH. */ + bp->radiative_efficiency = rad_efficiency(bp, props); + bp->jet_efficiency = jet_efficiency(bp, props); + + /* Calculate a BH angular momentum evolution time step. Two conditions are + used, one ensures that the BH spin changes by a small amount over the + next time-step, while the other ensures that the spin is not wildly + redirected between two time-steps */ if ((fabsf(bp->spin) > 0.01) && (bp->accretion_rate > 0.)) { - const float dt_ang_mom = 0.1 * fabsf(bp->spin) / - fabsf(da_dln_mbh_0(bp, constants, props)) * - bp->subgrid_mass / bp->accretion_rate; + + /* How much do we want to allow the spin to change over the next time- + step? If the spin is large (above 0.1 in magnitude), we allow it + to only change by 1% of its current value. If it is small, we instead + use a fixed increment of 0.01. This is to prevent the code from + becoming very slow. */ + const float spin_magnitude = fabsf(bp->spin); + float epsilon_spin = 0.01; + if (spin_magnitude > 0.1) { + epsilon_spin = 0.01 * spin_magnitude; + } + float dt_ang_mom = epsilon_spin / + fabsf(da_dln_mbh_0(bp, constants, props)) * + bp->subgrid_mass / bp->accretion_rate; + + /* We now compute the angular-momentum (direction) time-step. We allow + the anticipated warp angular momentum along the direction perpendicular + to the current spin vector to be 10% of the current BH angular momentum, + which should correspond to a change of direction of around 5 degrees. We + also apply this only if the spin is not low. */ + if (spin_magnitude > 0.1) { + + /* Angular momentum direction of gas in kernel along the direction of + the current spin vector */ + const float cosine = (bp->spec_angular_momentum_gas[0] * + bp->angular_momentum_direction[0] + + bp->spec_angular_momentum_gas[1] * + bp->angular_momentum_direction[1] + + bp->spec_angular_momentum_gas[2] * + bp->angular_momentum_direction[2]) / + spec_ang_mom_norm; + /* Compute sine, i.e. the componenent perpendicular to that. */ + const float sine = fmaxf(0., sqrtf(1. - cosine * cosine)); + + const float dt_redirection = 0.1 * m_warp(bp, constants, props) / + bp->accretion_rate * j_BH(bp, constants) / + j_warp(bp, constants, props) / sine; + dt_ang_mom = min(dt_ang_mom, dt_redirection); + } + bp->dt_ang_mom = dt_ang_mom; } else { bp->dt_ang_mom = FLT_MAX; @@ -1487,9 +1547,10 @@ INLINE static void black_holes_create_from_gas( /* Initial seed mass */ bp->subgrid_mass = props->subgrid_seed_mass; - /* Small initial spin in random direction*/ + /* Small initial spin magnitude */ bp->spin = props->seed_spin; + /* Generate a random unit vector for the spin direction */ const float rand_cos_theta = 2. * (0.5 - random_unit_interval(bp->id, ti_current, random_number_BH_spin)); @@ -1503,6 +1564,11 @@ INLINE static void black_holes_create_from_gas( bp->angular_momentum_direction[1] = rand_sin_theta * sin(rand_phi); bp->angular_momentum_direction[2] = rand_cos_theta; + /* Point the jets in the same direction to begin with */ + bp->jet_direction[0] = bp->angular_momentum_direction[0]; + bp->jet_direction[1] = bp->angular_momentum_direction[1]; + bp->jet_direction[2] = bp->angular_momentum_direction[2]; + bp->aspect_ratio = 0.01f; bp->jet_efficiency = 0.1f; bp->radiative_efficiency = 0.1f; diff --git a/src/black_holes/SPIN_JET/black_holes_iact.h b/src/black_holes/SPIN_JET/black_holes_iact.h index f252f9d0350e7ee282af4f171547febb295f6938..1bc1860d20ca684415320d19ac1d73d28edb3ca6 100644 --- a/src/black_holes/SPIN_JET/black_holes_iact.h +++ b/src/black_holes/SPIN_JET/black_holes_iact.h @@ -102,6 +102,11 @@ runner_iact_nonsym_bh_gas_density( /* Contribution to the smoothed sound speed */ bi->sound_speed_gas += mj * cj * wi; + if (cj * cosmo->a_factor_sound_speed > bh_props->sound_speed_hot_gas_min) { + bi->sound_speed_gas_hot += mj * cj * wi; + bi->rho_gas_hot += mj * wi; + } + /* Neighbour's (drifted) velocity in the frame of the black hole * (we do include a Hubble term) */ const float dv[3] = {pj->v[0] - bi->v[0], pj->v[1] - bi->v[1], @@ -293,7 +298,7 @@ runner_iact_nonsym_bh_gas_density( /* Check if the relative velocity is a significant fraction of the jet launching velocity. If it is, set the ray correction variable to some arbitrarily large value. */ - if (relative_velocity > 0.3 * bi->v_jet) { + if (relative_velocity > 0.8 * bi->v_jet) { ray_jet_correction = 1e11 * bi->h; } @@ -301,9 +306,7 @@ runner_iact_nonsym_bh_gas_density( the order closest --> farthest. */ if (cosine_theta < 0) { quantity_to_minimize = r + ray_jet_correction; - quantity_to_minimize_pos = 1e8 * r + ray_jet_correction; } else { - quantity_to_minimize = 1e8 * r + ray_jet_correction; quantity_to_minimize_pos = r + ray_jet_correction; } break; @@ -313,7 +316,7 @@ runner_iact_nonsym_bh_gas_density( /* Check if the relative velocity is a significant fraction of the jet launching velocity. If it is, set the ray correction variable to some arbitrarily large value. */ - if (relative_velocity > 0.3 * bi->v_jet) { + if (relative_velocity > 0.8 * bi->v_jet) { ray_jet_correction = 1e13 * 1. / bi->h; } @@ -321,9 +324,7 @@ runner_iact_nonsym_bh_gas_density( the order farthest --> closest */ if (cosine_theta < 0) { quantity_to_minimize = r_inv + ray_jet_correction; - quantity_to_minimize_pos = 1e8 * r_inv + ray_jet_correction; } else { - quantity_to_minimize = 1e8 * r_inv + ray_jet_correction; quantity_to_minimize_pos = r_inv + ray_jet_correction; } break; @@ -333,7 +334,7 @@ runner_iact_nonsym_bh_gas_density( /* Check if the relative velocity is a significant fraction of the jet launching velocity. If it is, set the ray correction variable to some arbitrarily large value. */ - if (relative_velocity > 0.3 * bi->v_jet) { + if (relative_velocity > 0.8 * bi->v_jet) { ray_jet_correction = 1e3; } @@ -341,8 +342,11 @@ runner_iact_nonsym_bh_gas_density( vector of the particle (relative to the BH) and the spin vector of the BH. I.e. we launch particles along the spin axis, regardless of the distances from the BH. */ - quantity_to_minimize = cosine_theta + ray_jet_correction; - quantity_to_minimize_pos = -cosine_theta + ray_jet_correction; + if (cosine_theta < 0) { + quantity_to_minimize = cosine_theta + ray_jet_correction; + } else { + quantity_to_minimize_pos = -cosine_theta + ray_jet_correction; + } break; } case AGN_jet_minimum_density_model: { @@ -350,7 +354,7 @@ runner_iact_nonsym_bh_gas_density( /* Check if the relative velocity is a significant fraction of the jet launching velocity. If it is, set the ray correction variable to some arbitrarily large value. */ - if (relative_velocity > 0.3 * bi->v_jet) { + if (relative_velocity > 0.8 * bi->v_jet) { ray_jet_correction = 1e15 * pj->rho; } @@ -358,19 +362,18 @@ runner_iact_nonsym_bh_gas_density( low-density gas. */ if (cosine_theta < 0) { quantity_to_minimize = pj->rho + ray_jet_correction; - quantity_to_minimize_pos = 1e8 * pj->rho + ray_jet_correction; } else { - quantity_to_minimize = 1e8 * pj->rho + ray_jet_correction; quantity_to_minimize_pos = pj->rho + ray_jet_correction; } break; } } - /* Loop over rays and do the actual minimization */ - for (int i = 0; i < spinjet_blackhole_number_of_rays; i++) { + /* Do the actual minimization. */ + if (cosine_theta < 0) { ray_minimise_distance(quantity_to_minimize, bi->rays_jet, arr_size_jet, gas_id, pj->mass); + } else { ray_minimise_distance(quantity_to_minimize_pos, bi->rays_jet_pos, arr_size_jet, gas_id, pj->mass); } @@ -1001,19 +1004,15 @@ runner_iact_nonsym_bh_gas_feedback( /* Get the (physical) kick velocity, and convert to code units */ const float vel_kick = sqrtf(2. * delta_u_jet) * cosmo->a; - /* Compute velocity kick direction using a function for generating a - * random unit vector within a cone around the spin vector.*/ + /* Compute velocity kick direction using the previously generated + * jet direction.*/ float vel_kick_direction[3]; - random_direction_in_cone(bi->id, pj->id, ti_current, - random_number_BH_kick, bh_props->opening_angle, - bi->angular_momentum_direction, - vel_kick_direction); /* Include the -1./1. factor (direction) which accounts for kicks in the * opposite direction of the spin vector */ - vel_kick_direction[0] = direction * vel_kick_direction[0]; - vel_kick_direction[1] = direction * vel_kick_direction[1]; - vel_kick_direction[2] = direction * vel_kick_direction[2]; + vel_kick_direction[0] = direction * bi->jet_direction[0]; + vel_kick_direction[1] = direction * bi->jet_direction[1]; + vel_kick_direction[2] = direction * bi->jet_direction[2]; /* Get the initial velocity */ const float v_init[3] = {xpj->v_full[0], xpj->v_full[1], xpj->v_full[2]}; @@ -1064,7 +1063,7 @@ runner_iact_nonsym_bh_gas_feedback( /* Store the jet energy */ const double delta_energy_jet = delta_u_jet * hydro_get_mass(pj); tracers_after_jet_feedback(pj, xpj, with_cosmology, cosmo->a, time, - delta_energy_jet); + delta_energy_jet, vel_kick); /* Impose maximal viscosity */ hydro_diffusive_feedback_reset(pj); diff --git a/src/black_holes/SPIN_JET/black_holes_part.h b/src/black_holes/SPIN_JET/black_holes_part.h index 712c70332230f77ec5747554bacb37e1cffc866e..967449be2b062e278f511d36144ec6bf7cffa3e5 100644 --- a/src/black_holes/SPIN_JET/black_holes_part.h +++ b/src/black_holes/SPIN_JET/black_holes_part.h @@ -206,6 +206,9 @@ struct bpart { /*! The normalized spin/angular momentum vector of the BH */ float angular_momentum_direction[3]; + /* The current jet direction. */ + float jet_direction[3]; + /*! The jet efficiency */ float jet_efficiency; @@ -252,6 +255,11 @@ struct bpart { /*! Halo mass the black hole is assigned to */ float group_mass; + /*! The density and sound speed of hot gas used to calculate variable + * jet velocities. */ + float rho_gas_hot; + float sound_speed_gas_hot; + union { /*! Last time a jet event occurred */ diff --git a/src/black_holes/SPIN_JET/black_holes_properties.h b/src/black_holes/SPIN_JET/black_holes_properties.h index fd01037c5d92d3ea2d72f984ccfc2b3db67c0b58..9a044d297abd22daf8b9390453af56541de0e38d 100644 --- a/src/black_holes/SPIN_JET/black_holes_properties.h +++ b/src/black_holes/SPIN_JET/black_holes_properties.h @@ -42,12 +42,11 @@ enum AGN_jet_velocity_models { AGN_jet_velocity_BH_mass, /*< Scale the jet velocity with BH mass */ AGN_jet_velocity_mass_loading, /*< Assume constant mass loading */ AGN_jet_velocity_local, /*< Scale the jet velocity such that particles clear - the kernel exactly when a new one is launched */ - AGN_jet_velocity_sound_speed, /*< Scale the jet velocity such that the - BH kernel is replenished on the same - time-scale as it is evacuated by jets */ - AGN_jet_velocity_halo_mass /*< Scale the jet velocity with halo mass - (NOT WORKING) */ + the kernel exactly when a new one is launched, + as well as making sure that the kernel never + runs out of particles */ + AGN_jet_velocity_halo_mass /*< Scale the jet velocity with halo mass + (NOT WORKING) */ }; enum BH_merger_thresholds { @@ -297,11 +296,6 @@ struct black_holes_props { free-free absorption dominates the opacity */ enum thin_disc_regions TD_region; - /*! Parameter controlling when thin disk transitions into slim disk. This - occurs when the radiative efficiency of the slim disk falls below - TD_SD_eps_r_threshold times the radiative efficiency of the thin disk. */ - float TD_SD_eps_r_threshold; - /* ---- Jet feedback - related parameters ---------- */ /*! Global switch for whether to include jets [1] or not [0]. */ @@ -322,18 +316,6 @@ struct black_holes_props { /* Whether to use GRMHD fits for the spindown rate due to jets */ int include_GRMHD_spindown; - /* Whether to include the expected suppression of the accretion rate due to - * ADIOS winds in the thick disk regime */ - int include_ADIOS_suppression; - - /* The inner radius in the accretion rate - R relation, if ADIOS suppression - * is included */ - float ADIOS_R_in; - - /* The slope of the accretion rate - R relation, if ADIOS suppression - * is included */ - float ADIOS_s; - /*! Whether to fix the radiative efficiency to some value [1] or not [0]. */ int fix_radiative_efficiency; @@ -372,6 +354,13 @@ struct black_holes_props { /*! The minimal jet velocity to use in the variable-velocity models */ float v_jet_min; + /*! The maximal jet velocity to use in the variable-velocity models */ + float v_jet_max; + + /*! The minimum sound speed of hot gas to count when calculating the + * smoothed sound speed of gas in the kernel */ + float sound_speed_hot_gas_min; + /*! The effective (half-)opening angle of the jet. */ float opening_angle; @@ -383,13 +372,19 @@ struct black_holes_props { /*! The coupling efficiency for jet feedback. */ float eps_f_jet; - /*! Whether to fix the jet efficiency to some value [1] or not [0]. If yes, - the jets will be pointed along the z-axis. */ + /*! Whether to fix the jet efficiency to some value [1] or not [0]. */ int fix_jet_efficiency; /*! The jet efficiency to use if fix_jet_efficiency is 1. */ float jet_efficiency; + /* Whether to fix the jet directions to be along the z-axis. */ + int fix_jet_direction; + + /*! The accretion efficiency (suppression of accretion rate) to use in + * the thick disc regime (at low Eddington ratios). */ + float accretion_efficiency; + /*! The jet launching scheme to use: minimum distance, maximum distance, closest to spin axis or minimum density. */ enum AGN_jet_feedback_models jet_feedback_model; @@ -785,45 +780,6 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp, bp->include_GRMHD_spindown); } - bp->include_ADIOS_suppression = - parser_get_param_int(params, "SPINJETAGN:include_ADIOS_suppression"); - - if ((bp->include_ADIOS_suppression != 0) && - (bp->include_ADIOS_suppression != 1)) { - error( - "The include_ADIOS_suppression parameter must be either 0 or 1, " - "not %d", - bp->include_ADIOS_suppression); - } - - bp->ADIOS_R_in = parser_get_param_float(params, "SPINJETAGN:ADIOS_R_in"); - - if (bp->ADIOS_R_in <= 1.) { - error( - "The ADIOS_R_in parameter must be > 1, " - "not %f", - bp->ADIOS_R_in); - } - - bp->ADIOS_s = parser_get_param_float(params, "SPINJETAGN:ADIOS_s"); - - if ((bp->ADIOS_s < 0.) || (bp->ADIOS_s > 1.)) { - error( - "The ADIOS_s parameter must be between 0 and 1, " - "not %f", - bp->ADIOS_s); - } - - bp->TD_SD_eps_r_threshold = - parser_get_param_float(params, "SPINJETAGN:TD_SD_eps_r_threshold"); - - if ((bp->TD_SD_eps_r_threshold <= 0.) || (bp->TD_SD_eps_r_threshold >= 1.)) { - error( - "The TD_SD_eps_r_threshold parameter governing the transition between" - "thin and slim disk must be between 0. and 1., not %f", - bp->TD_SD_eps_r_threshold); - } - bp->N_jet = parser_get_param_float(params, "SPINJETAGN:N_jet"); if (bp->N_jet % 2 != 0) { @@ -859,6 +815,10 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp, parser_get_param_float(params, "SPINJETAGN:v_jet_min_km_p_s"); bp->v_jet_min *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); + bp->v_jet_max = + parser_get_param_float(params, "SPINJETAGN:v_jet_max_km_p_s"); + bp->v_jet_max *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); + } else if (strcmp(temp4, "MassLoading") == 0) { bp->AGN_jet_velocity_model = AGN_jet_velocity_mass_loading; @@ -869,6 +829,10 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp, parser_get_param_float(params, "SPINJETAGN:v_jet_min_km_p_s"); bp->v_jet_min *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); + bp->v_jet_max = + parser_get_param_float(params, "SPINJETAGN:v_jet_max_km_p_s"); + bp->v_jet_max *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); + } else if (strcmp(temp4, "Local") == 0) { bp->AGN_jet_velocity_model = AGN_jet_velocity_local; @@ -878,14 +842,18 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp, parser_get_param_float(params, "SPINJETAGN:v_jet_min_km_p_s"); bp->v_jet_min *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); - } else if (strcmp(temp4, "SoundSpeed") == 0) { - bp->AGN_jet_velocity_model = AGN_jet_velocity_sound_speed; + bp->v_jet_max = + parser_get_param_float(params, "SPINJETAGN:v_jet_max_km_p_s"); + bp->v_jet_max *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); - bp->v_jet_xi = parser_get_param_float(params, "SPINJETAGN:v_jet_xi"); + float temp_hot_gas_min = + parser_get_param_float(params, "SPINJETAGN:temperature_hot_gas_min_K"); - bp->v_jet_min = - parser_get_param_float(params, "SPINJETAGN:v_jet_min_km_p_s"); - bp->v_jet_min *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs)); + temp_hot_gas_min /= units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE); + + bp->sound_speed_hot_gas_min = + sqrtf(hydro_gamma * hydro_gamma_minus_one * temp_hot_gas_min * + bp->temp_to_u_factor); } else if (strcmp(temp4, "HaloMass") == 0) { error( @@ -928,10 +896,26 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp, parser_get_param_float(params, "SPINJETAGN:jet_efficiency"); if (bp->jet_efficiency <= 0.) { + error("The constant jet efficiency parameter must be >0., not %f", + bp->jet_efficiency); + } + + bp->fix_jet_direction = + parser_get_param_int(params, "SPINJETAGN:fix_jet_direction"); + + if ((bp->fix_jet_direction != 0) && (bp->fix_jet_direction != 1)) { error( - "The jet_efficiency corresponding to the jet efficiency " - "must be larger than 0., not %f", - bp->jet_efficiency); + "The fix_jet_direction parameter must be either 0 or 1, " + "not %d", + bp->fix_jet_direction); + } + + bp->accretion_efficiency = + parser_get_param_float(params, "SPINJETAGN:accretion_efficiency"); + + if (bp->accretion_efficiency <= 0.) { + error("The accretion efficiency must be larger than 0., not %f", + bp->accretion_efficiency); } bp->fix_radiative_efficiency = diff --git a/src/black_holes/SPIN_JET/black_holes_spin.h b/src/black_holes/SPIN_JET/black_holes_spin.h index d7207ef770f31a826946713aed7ae39483fd71e4..d7e611b0cb500467bb8a87d5a0b443338f181fe5 100644 --- a/src/black_holes/SPIN_JET/black_holes_spin.h +++ b/src/black_holes/SPIN_JET/black_holes_spin.h @@ -95,8 +95,8 @@ __attribute__((always_inline)) INLINE static float j_BH( struct bpart* bp, const struct phys_const* constants) { const float J_BH = - fabs(bp->subgrid_mass * bp->subgrid_mass * bp->spin * - constants->const_newton_G / constants->const_speed_light_c); + fabsf(bp->subgrid_mass * bp->subgrid_mass * bp->spin * + constants->const_newton_G / constants->const_speed_light_c); #ifdef SWIFT_DEBUG_CHECKS if (J_BH <= 0.) { @@ -491,11 +491,28 @@ __attribute__((always_inline)) INLINE static float eps_SD(float a, float mdot) { */ __attribute__((always_inline)) INLINE static void decide_mode( struct bpart* bp, const struct black_holes_props* props) { - if (bp->eddington_fraction < props->mdot_crit_ADAF) { + + /* For deciding the accretion mode, we want to use the Eddington fraction + * calculated using the raw, unsuppressed accretion rate. This means that + * if the disc is currently thick, its current Eddington fraction, which is + * already suppressed, needs to be unsuppressed (increased) to retrieve the + * raw Bondi-based Eddington ratio. */ + float eddington_fraction_Bondi = bp->eddington_fraction; + if (bp->accretion_mode == BH_thick_disc) + eddington_fraction_Bondi *= 1. / props->accretion_efficiency; + + if (eddington_fraction_Bondi < props->mdot_crit_ADAF) { bp->accretion_mode = BH_thick_disc; } else { - if ((eps_SD(bp->spin, bp->eddington_fraction) < - props->TD_SD_eps_r_threshold * eps_NT(bp->spin)) && + + /* The disc is assumed to be slim (super-Eddington) if the Eddington + * fraction is above 1. However, the Eddington fraction is typically + * (also in this code and model) defined using an Eddington luminosity + * with an assumed radiative efficiency of eps_r = 0.1. For this purpose + * it is more sensible to define the Eddington luminosity using the spin- + * dependent eps_NT(a). We thus assume that the disc becomes super- + * Eddington at a critical f_Edd(a) = f_Edd * (eps_NT(a)/0.1) = 1. */ + if ((eddington_fraction_Bondi * eps_NT(bp->spin) / 0.1 > 1.) && (props->include_slim_disk)) { bp->accretion_mode = BH_slim_disc; } else { @@ -538,14 +555,12 @@ __attribute__((always_inline)) INLINE static float aspect_ratio( /* Define placeholder variable for the result */ float h_0 = -1.; - /* Start branching depending on which accretion mode the BH is in */ + /* Start branching depending on which accretion mode the BH is in. + * We assume the thick and slim disc to be very similar in geometrical + * terms. */ if ((bp->accretion_mode == BH_thick_disc) || (bp->accretion_mode == BH_slim_disc)) { - if (bp->accretion_mode == BH_thick_disc) { - h_0 = props->h_0_ADAF; - } else { - h_0 = 0.5 * props->gamma_SD_inv; - } + h_0 = props->h_0_ADAF; } else { /* Start branching depending on which region of the thin disk we wish to @@ -871,8 +886,10 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet( sqrtf(halo_mass * constants->const_newton_G / virial_radius); const float sound_speed = sqrtf(5. / 3. * 0.5) * virial_velocity; - /* Return the jet velocity as some factor times the sound speed */ + /* Return the jet velocity as some factor times the sound speed, apply + * floor and ceiling values. */ v_jet = fmaxf(props->v_jet_min, props->v_jet_cs_ratio * sound_speed); + v_jet = fminf(props->v_jet_max, v_jet); } else if (props->AGN_jet_velocity_model == AGN_jet_velocity_constant) { v_jet = props->v_jet; @@ -883,33 +900,56 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet( apply a floor value*/ v_jet = sqrtf(2.f * bp->jet_efficiency / props->v_jet_mass_loading) * constants->const_speed_light_c; + + /* Apply floor and ceiling values */ v_jet = fmaxf(props->v_jet_min, v_jet); + v_jet = fminf(props->v_jet_max, v_jet); } else if (props->AGN_jet_velocity_model == AGN_jet_velocity_local) { /* Calculate jet power */ - const float jet_power = bp->jet_efficiency * bp->accretion_rate * - constants->const_speed_light_c * - constants->const_speed_light_c; - - /* Calculate jet velocity from the power, smoothing length (proper, not - comoving) and neighbour mass. Then apply floor. */ - v_jet = props->v_jet_xi * cbrt(jet_power * bp->h * cosmo->a / - (bp->ngb_mass / ((double)bp->num_ngbs))); - v_jet = fmaxf(v_jet, props->v_jet_min); - - } else if (props->AGN_jet_velocity_model == AGN_jet_velocity_sound_speed) { - - /* Calculate jet power */ - const float jet_power = bp->jet_efficiency * bp->accretion_rate * - constants->const_speed_light_c * - constants->const_speed_light_c; - - /* Calculate jet velocity from the power, smoothing length (proper, not - comoving), neighbour sound speed and neighbour mass. Apply floor. */ - v_jet = props->v_jet_xi * sqrtf(jet_power * bp->h * cosmo->a / - (bp->ngb_mass * bp->sound_speed_gas)); + const double jet_power = bp->jet_efficiency * bp->accretion_rate * + constants->const_speed_light_c * + constants->const_speed_light_c; + + /* Get the sound speed of the hot gas in the kernel. Make sure the actual + * value that is used is at least the value specified in the parameter + * file. */ + float sound_speed_hot_gas = + bp->sound_speed_gas_hot * cosmo->a_factor_sound_speed; + sound_speed_hot_gas = + max(sound_speed_hot_gas, props->sound_speed_hot_gas_min); + + /* Take the maximum of the sound speed of the hot gas and the gas velocity + * dispersion. Calculate the replenishment time-scale by assuming that it + * will replenish under the influence of whichever of those two values is + * larger. */ + const float gas_dispersion = bp->velocity_dispersion_gas * cosmo->a_inv; + const double replenishment_time_scale = + bp->h * cosmo->a / max(sound_speed_hot_gas, gas_dispersion); + + /* Calculate jet velocity from the replenishment condition, taking the + * power, smoothing length (proper, not comoving), neighbour sound speed + * and (total) neighbour mass. */ + const float v_jet_repl = + sqrtf(jet_power * replenishment_time_scale / (2. * bp->ngb_mass)); + + /* Calculate jet velocity from the crossing condition, i.e. set the + * velocity such that a new particle pair will be launched roughly when + * the previous one crosses (exits) the BH kernel. This also depends on + * power, smoothing length and neighbour mass (per particle, not total). */ + const float v_jet_cross = + cbrtf(bp->h * cosmo->a * jet_power / + (4. * bp->ngb_mass / ((double)bp->num_ngbs))); + + /* Find whichever of these two is the minimum, and multiply it by an + * arbitrary scaling factor (whose fiducial value is 1, i.e. no + * rescaling. */ + v_jet = props->v_jet_xi * fmaxf(v_jet_repl, v_jet_cross); + + /* Apply floor and ceiling values */ v_jet = fmaxf(v_jet, props->v_jet_min); + v_jet = fminf(v_jet, props->v_jet_max); } else { error( @@ -1059,7 +1099,6 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve( const float mass_ratio = m2 / m1; const float sym_mass_ratio = mass_ratio / ((mass_ratio + 1.f) * (mass_ratio + 1.f)); - const float reduced_mass = m1 * m2 / (m1 + m2); /* The absolute values of the spins are also needed */ const float spin1 = fabsf(bpi->spin); @@ -1084,21 +1123,54 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve( /* We want to compute the direction of the orbital angular momentum of the two BHs, which is used in the fits. Start by defining the coordinates in - the frame of one of the BHs (it doesn't matter which one, the total - angular momentum is the same). */ - const float relative_coordinates[3] = { - bpj->x[0] - bpi->x[0], bpj->x[1] - bpi->x[1], bpj->x[2] - bpi->x[2]}; - const float relative_velocities[3] = { - bpj->v[0] - bpi->v[0], bpj->v[1] - bpi->v[1], bpj->v[2] - bpi->v[2]}; - - /* Calculate the orbital angular momentum itself. */ + the frame of the centre of mass. For this we first need to compute the + centre of mass coodinates and velocity. */ + const float centre_of_mass[3] = { + (m1 * bpi->x[0] + m2 * bpj->x[0]) / (m1 + m2), + (m1 * bpi->x[1] + m2 * bpj->x[1]) / (m1 + m2), + (m1 * bpi->x[2] + m2 * bpj->x[2]) / (m1 + m2)}; + const float centre_of_mass_vel[3] = { + (m1 * bpi->v[0] + m2 * bpj->v[0]) / (m1 + m2), + (m1 * bpi->v[1] + m2 * bpj->v[1]) / (m1 + m2), + (m1 * bpi->v[2] + m2 * bpj->v[2]) / (m1 + m2)}; + + /* Coordinates of each of the BHs in the frame of the centre of mass. */ + const float relative_coordinates_1[3] = {bpi->x[0] - centre_of_mass[0], + bpi->x[1] - centre_of_mass[1], + bpi->x[2] - centre_of_mass[2]}; + const float relative_coordinates_2[3] = {bpj->x[0] - centre_of_mass[0], + bpj->x[1] - centre_of_mass[1], + bpj->x[2] - centre_of_mass[2]}; + + /* The velocities of each BH in the centre of mass frame. */ + const float relative_velocities_1[3] = {bpi->v[0] - centre_of_mass_vel[0], + bpi->v[1] - centre_of_mass_vel[1], + bpi->v[2] - centre_of_mass_vel[2]}; + const float relative_velocities_2[3] = {bpj->v[0] - centre_of_mass_vel[0], + bpj->v[1] - centre_of_mass_vel[1], + bpj->v[2] - centre_of_mass_vel[2]}; + + /* The angular momentum of each BH in the centre of mass frame. */ + const float angular_momentum_1[3] = { + m1 * (relative_coordinates_1[1] * relative_velocities_1[2] - + relative_coordinates_1[2] * relative_velocities_1[1]), + m1 * (relative_coordinates_1[2] * relative_velocities_1[0] - + relative_coordinates_1[0] * relative_velocities_1[2]), + m1 * (relative_coordinates_1[0] * relative_velocities_1[1] - + relative_coordinates_1[1] * relative_velocities_1[0])}; + const float angular_momentum_2[3] = { + m2 * (relative_coordinates_2[1] * relative_velocities_2[2] - + relative_coordinates_2[2] * relative_velocities_2[1]), + m2 * (relative_coordinates_2[2] * relative_velocities_2[0] - + relative_coordinates_2[0] * relative_velocities_2[2]), + m2 * (relative_coordinates_2[0] * relative_velocities_2[1] - + relative_coordinates_2[1] * relative_velocities_2[0])}; + + /* The total, final orbital angular momentum. */ const float orbital_angular_momentum[3] = { - reduced_mass * (relative_coordinates[1] * relative_velocities[2] - - relative_coordinates[2] * relative_velocities[1]), - reduced_mass * (relative_coordinates[2] * relative_velocities[0] - - relative_coordinates[0] * relative_velocities[2]), - reduced_mass * (relative_coordinates[0] * relative_velocities[1] - - relative_coordinates[1] * relative_velocities[0])}; + angular_momentum_1[0] + angular_momentum_2[0], + angular_momentum_1[1] + angular_momentum_2[1], + angular_momentum_1[2] + angular_momentum_2[2]}; /* Calculate the magnitude of the orbital angular momentum. */ const float orbital_angular_momentum_magnitude = @@ -1121,12 +1193,18 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve( system, i.e. including the orbital angular momentum and the spins. This is needed since the final spin is assumed to be along the direction of this total angular momentum. Hence here we compute the direction. */ + const float j_BH_1 = + fabsf(bpi->subgrid_mass * bpi->subgrid_mass * bpi->spin * + constants->const_newton_G / constants->const_speed_light_c); + const float j_BH_2 = + fabsf(bpj->subgrid_mass * bpj->subgrid_mass * bpj->spin * + constants->const_newton_G / constants->const_speed_light_c); float total_angular_momentum_direction[3] = { - m1 * spin1 * spin_vec1[0] + m2 * spin2 * spin_vec2[0] + + j_BH_1 * spin_vec1[0] + j_BH_2 * spin_vec2[0] + orbital_angular_momentum[0], - m1 * spin1 * spin_vec1[1] + m2 * spin2 * spin_vec2[1] + + j_BH_1 * spin_vec1[1] + j_BH_2 * spin_vec2[1] + orbital_angular_momentum[1], - m1 * spin1 * spin_vec1[2] + m2 * spin2 * spin_vec2[2] + + j_BH_1 * spin_vec1[2] + j_BH_2 * spin_vec2[2] + orbital_angular_momentum[2]}; /* The above is actually the total angular momentum, so we need to normalize diff --git a/src/black_holes_properties.h b/src/black_holes_properties.h index f0f7ae660ebbb56cb7949e517e9f334a137e53c3..c54db15a9317fcd13ac4efa6e6da50344258317f 100644 --- a/src/black_holes_properties.h +++ b/src/black_holes_properties.h @@ -33,4 +33,11 @@ #error "Invalid choice of black hole model" #endif +/*! Define a flag to be used to output jet tracer data (or not) */ +#if defined(BLACK_HOLES_SPIN_JET) +#define with_jets 1 +#else +#define with_jets 0 +#endif + #endif /* SWIFT_BLACK_HOLES_PROPERTIES_H */ diff --git a/src/random.h b/src/random.h index 530233b68fb4b32e2961e6112c90f50c778f7e3b..9b5d75ddc1338e3b78bdcddb8ad53f4b194f14f5 100644 --- a/src/random.h +++ b/src/random.h @@ -300,7 +300,6 @@ INLINE static int random_poisson(const int64_t id, const double lambda, * opening_angle radians. * * @param id_bh The ID of the (BH) particle which is doing the jet feedback. - * @param id_gas The ID of the gas particle being kicked by jet feedback. * @param ti_current The time (on the time-line) for which to generate the * random kick direction. * @param type The #random_number_type to generate. @@ -308,10 +307,12 @@ INLINE static int random_poisson(const int64_t id, const double lambda, * @param a Reference direction that defines the cone. * @param rand_cone_direction Return value. */ -INLINE static void random_direction_in_cone( - const int64_t id_bh, const int64_t id_gas, const integertime_t ti_current, - const enum random_number_type type, const float opening_angle, - const float a[3], float rand_cone_direction[3]) { +INLINE static void random_direction_in_cone(const int64_t id_bh, + const integertime_t ti_current, + const enum random_number_type type, + const float opening_angle, + const float a[3], + float rand_cone_direction[3]) { /* We want to draw a random unit vector from a cone around the unit * vector a = (a0, a1, a2). We do this in a frame x'y'z', where the z' @@ -367,17 +368,15 @@ INLINE static void random_direction_in_cone( /* Draw a random cosine confined to the range [cos(opening_angle), 1] */ const float rand_cos_theta = 1.f - (1.f - cosf(opening_angle)) * - random_unit_interval(id_bh + id_gas, ti_current, type); + random_unit_interval(id_bh, ti_current, type); /* Get the corresponding sine */ const float rand_sin_theta = sqrtf(max(0.f, (1.f - rand_cos_theta) * (1.f + rand_cos_theta))); /* Get a random equitorial angle from [0, 180] deg */ - const float rand_phi = - ((float)(2. * M_PI)) * - random_unit_interval((id_bh + id_gas) * (id_bh + id_gas), ti_current, - type); + const float rand_phi = ((float)(2. * M_PI)) * + random_unit_interval(id_bh * id_bh, ti_current, type); float rand_cos_phi, rand_sin_phi; sincosf(rand_phi, &rand_sin_phi, &rand_cos_phi); diff --git a/src/tracers/EAGLE/tracers.h b/src/tracers/EAGLE/tracers.h index 6f135882acdcc1077c645228d154e9085ffcec74..aa74338bd9aaee5264daa0457deaba4b2cfb1f17 100644 --- a/src/tracers/EAGLE/tracers.h +++ b/src/tracers/EAGLE/tracers.h @@ -214,6 +214,12 @@ static INLINE void tracers_first_init_xpart( xp->tracers_data.last_AGN_injection_scale_factor = -1.f; xp->tracers_data.density_at_last_AGN_feedback_event = -1.f; + + xp->tracers_data.hit_by_jet_feedback = 0; + xp->tracers_data.jet_feedback_energy = 0.f; + xp->tracers_data.last_AGN_jet_feedback_scale_factor = 0.f; + xp->tracers_data.last_AGN_jet_feedback_time = 0.f; + xp->tracers_data.last_jet_kick_velocity = 0.f; } /** @@ -317,7 +323,8 @@ static INLINE void tracers_after_black_holes_feedback( static INLINE void tracers_after_jet_feedback( const struct part *p, struct xpart *xp, const int with_cosmology, - const float scale_factor, const double time, const double delta_energy) { + const float scale_factor, const double time, const double delta_energy, + const float vel_kick) { if (with_cosmology) xp->tracers_data.last_AGN_jet_feedback_scale_factor = scale_factor; @@ -325,6 +332,7 @@ static INLINE void tracers_after_jet_feedback( xp->tracers_data.last_AGN_jet_feedback_time = time; xp->tracers_data.hit_by_jet_feedback++; xp->tracers_data.jet_feedback_energy += delta_energy; + xp->tracers_data.last_jet_kick_velocity = vel_kick; } /** diff --git a/src/tracers/EAGLE/tracers_io.h b/src/tracers/EAGLE/tracers_io.h index 9997449151cba56ffd542c2c2d9f28a5e1586e60..3e1ab3948e8278974107dc81e5df2ff66c80ac94 100644 --- a/src/tracers/EAGLE/tracers_io.h +++ b/src/tracers/EAGLE/tracers_io.h @@ -24,6 +24,7 @@ #include <config.h> /* Local includes */ +#include "black_holes_properties.h" #include "io_properties.h" #include "tracers.h" @@ -187,7 +188,50 @@ __attribute__((always_inline)) INLINE static int tracers_write_particles( "Star formation rates of the particles averaged over the period set by " "the first two snapshot triggers"); - return 11; + if (with_jets) { + list[11] = io_make_output_field( + "KickedByJetFeedback", CHAR, 1, UNIT_CONV_NO_UNITS, 0.f, xparts, + tracers_data.hit_by_jet_feedback, + "Flags the particles that have been directly kicked by" + "an AGN jet feedback event at some point in the past. " + "If > 0, contains the number of individual events."); + + list[12] = io_make_output_field("EnergiesReceivedFromJetFeedback", FLOAT, 1, + UNIT_CONV_ENERGY, 0.f, xparts, + tracers_data.jet_feedback_energy, + "Total amount of kinetic energy from AGN " + "jet feedback events received by the " + "particles while they were still gas " + "particles."); + + if (with_cosmology) { + + list[13] = io_make_output_field( + "LastAGNJetFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, + xparts, tracers_data.last_AGN_jet_feedback_scale_factor, + "Scale-factors at which the particles were last hit by jet " + "feedback while they were still gas particles. " + "-1 if a particle has never been hit by feedback"); + + } else { + + list[13] = io_make_output_field( + "LastAGNJetFeedbackTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f, xparts, + tracers_data.last_AGN_jet_feedback_time, + "Times at which the particles were last hit by jet" + "feedback while they were still gas particles. " + "-1 if a particle has never been hit by feedback"); + } + + list[14] = io_make_output_field("LastAGNJetKickVelocities", FLOAT, 1, + UNIT_CONV_VELOCITY, 0.f, xparts, + tracers_data.last_jet_kick_velocity, + "Kick velocity at last AGN jet event."); + + return 15; + } else { + return 11; + } } __attribute__((always_inline)) INLINE static int tracers_write_sparticles( @@ -292,7 +336,50 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles( "the first two snapshot triggers when the particle was still a gas " "particle."); - return 11; + if (with_jets) { + list[14] = io_make_output_field( + "KickedByJetFeedback", CHAR, 1, UNIT_CONV_NO_UNITS, 0.f, sparts, + tracers_data.hit_by_jet_feedback, + "Flags the particles that have been directly kicked by" + "an AGN jet feedback event at some point in the past. " + "If > 0, contains the number of individual events."); + + list[15] = io_make_output_field("EnergiesReceivedFromJetFeedback", FLOAT, 1, + UNIT_CONV_ENERGY, 0.f, sparts, + tracers_data.jet_feedback_energy, + "Total amount of kinetic energy from AGN " + "jet feedback events received by the " + "particles while they were still gas " + "particles."); + + if (with_cosmology) { + + list[16] = io_make_output_field( + "LastAGNJetFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, + sparts, tracers_data.last_AGN_jet_feedback_scale_factor, + "Scale-factors at which the particles were last hit by jet " + "feedback while they were still gas particles. " + "-1 if a particle has never been hit by feedback"); + + } else { + + list[16] = io_make_output_field( + "LastAGNJetFeedbackTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f, sparts, + tracers_data.last_AGN_jet_feedback_time, + "Times at which the particles were last hit by jet" + "feedback while they were still gas particles. " + "-1 if a particle has never been hit by feedback"); + } + + list[17] = io_make_output_field("LastAGNJetKickVelocities", FLOAT, 1, + UNIT_CONV_VELOCITY, 0.f, sparts, + tracers_data.last_jet_kick_velocity, + "Kick velocity at last AGN jet event."); + + return 18; + } else { + return 15; + } } __attribute__((always_inline)) INLINE static int tracers_write_bparticles( diff --git a/src/tracers/EAGLE/tracers_struct.h b/src/tracers/EAGLE/tracers_struct.h index d1203ddd9fc54dd8d050c966da119873b64ab649..5e8c7df73249f4ac0149c7dcfc989a54a894acc9 100644 --- a/src/tracers/EAGLE/tracers_struct.h +++ b/src/tracers/EAGLE/tracers_struct.h @@ -94,6 +94,9 @@ struct tracers_xpart_data { /*! Has this particle been hit by AGN feedback? */ char hit_by_AGN_feedback; + + /* Kick velocity at last AGN jet event */ + float last_jet_kick_velocity; }; /** diff --git a/tests/testRandomCone.c b/tests/testRandomCone.c index 2393362bad991a2fa93aecda15a9d326b9ed82de..08c953c00acf2adb6997be010b44ed045229b5ab 100644 --- a/tests/testRandomCone.c +++ b/tests/testRandomCone.c @@ -67,11 +67,8 @@ void test_cone(int64_t id_bh, const integertime_t ti_current, for (int k = 0; k < N_test; ++k) { - /* Generate random ids. */ - const long long id_p = rand() * (1LL << 31) + rand(); - /* Generate a random unit vector within a cone around unit_vector */ - random_direction_in_cone(id_p, id_bh, ti_current, type, opening_angle, + random_direction_in_cone(id_bh, ti_current, type, opening_angle, unit_vector, rand_vector); /* Check that this vector is actually within the cone we want */