diff --git a/examples/EAGLE_low_z/EAGLE_6/imf_example_6.yml b/examples/EAGLE_low_z/EAGLE_6/imf_example_6.yml new file mode 100644 index 0000000000..75a737c535 --- /dev/null +++ b/examples/EAGLE_low_z/EAGLE_6/imf_example_6.yml @@ -0,0 +1,271 @@ +MetaData: + run_name: EAGLE-L0006N0094-Ref + +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98841e43 # 10^10 M_sun in grams + UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters + UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Structure finding options +StructureFinding: + config_file_name: stf_input.cfg # Name of the STF config file. + basename: ./stf # Common part of the name of output files. + scale_factor_first: 0.92 # Scale-factor of the first snaphot (cosmological run) + time_first: 0.01 # Time of the first structure finding output (in internal units). + delta_time: 1.10 # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals. + +# Cosmological parameters +Cosmology: + h: 0.6777 # Reduced Hubble constant + a_begin: 0.9090909 # Initial scale-factor of the simulation + a_end: 1.0 # Final scale factor of the simulation + Omega_cdm: 0.2587481 # Cold Dark Matter density parameter + Omega_lambda: 0.693 # Dark-energy density parameter + Omega_b: 0.0482519 # Baryon density parameter + +Scheduler: + max_top_level_cells: 8 + links_per_tasks: 500 + +# Parameters governing the time integration +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 1e-2 # The end time of the simulation (in internal units). + dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). + dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: eagle # Common part of the name of output files + scale_factor_first: 0.91 # Scale-factor of the first snaphot (cosmological run) + time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units) + delta_time: 1.01 # Time difference between consecutive outputs (in internal units) + compression: 1 + distributed: 1 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + +# Parameters governing the conserved quantities statistics +Statistics: + scale_factor_first: 0.92 # Scale-factor of the first stat dump (cosmological run) + time_first: 0.01 # Time of the first stat dump (non-cosmological run) (in internal units) + delta_time: 1.05 # Time between statistics output + +# Parameters for the self-gravity scheme +Gravity: + eta: 0.025 # Constant dimensionless multiplier for time integration. + MAC: adaptive + epsilon_fmm: 0.001 + theta_cr: 0.7 # Opening angle (Multipole acceptance criterion) + mesh_side_length: 64 + comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units). + max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units). + comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units). + max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units). + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + h_min_ratio: 0.1 # Minimal smoothing length in units of softening. + h_max: 0.5 # Maximal smoothing length in co-moving internal units. + CFL_condition: 0.2 # Courant-Friedrich-Levy condition for time integration. + minimal_temperature: 100 # (internal units) + particle_splitting: 1 + particle_splitting_mass_threshold: 7e-4 # Internal units (i.e. 7e6 Msun ~ 4 times the initial gas mass) + +# Parameters of the stars neighbour search +Stars: + resolution_eta: 1.1642 # Target smoothing length in units of the mean inter-particle separation + h_tolerance: 7e-3 + overwrite_birth_time: 1 + birth_time: 0.33333 # Pretend all the stars were born at z = 2 + luminosity_filename: ./photometry + +# Parameters for the Friends-Of-Friends algorithm +FOF: + basename: fof_output # Filename for the FOF outputs. + min_group_size: 256 # The minimum no. of particles required for a group. + linking_length_ratio: 0.2 # Linking length in units of the main inter-particle separation. + seed_black_holes_enabled: 1 # Enable seeding of black holes in FoF groups + black_hole_seed_halo_mass_Msun: 1.5e10 # Minimal halo mass in which to seed a black hole (in solar masses). + scale_factor_first: 0.91 # Scale-factor of first FoF black hole seeding calls. + delta_time: 1.005 # Scale-factor ratio between consecutive FoF black hole seeding calls. + linking_types: [0, 1, 0, 0, 0, 0, 0] + attaching_types: [1, 0, 0, 0, 1, 1, 0] + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./EAGLE_ICs_6.hdf5 # The file to read + periodic: 1 + cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget + cleanup_velocity_factors: 1 # Remove the sqrt(a) factor in the velocities inherited from Gadget + remap_ids: 1 + +EAGLEChemistry: # Solar abundances + init_abundance_metal: 0.014 + init_abundance_Hydrogen: 0.70649785 + init_abundance_Helium: 0.28055534 + init_abundance_Carbon: 2.0665436e-3 + init_abundance_Nitrogen: 8.3562563e-4 + init_abundance_Oxygen: 5.4926244e-3 + init_abundance_Neon: 1.4144605e-3 + init_abundance_Magnesium: 5.907064e-4 + init_abundance_Silicon: 6.825874e-4 + init_abundance_Iron: 1.1032152e-3 + +# EAGLE cooling parameters +EAGLECooling: + dir_name: ./coolingtables/ + H_reion_z: 11.5 + H_reion_eV_p_H: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_eV_p_H: 2.0 + +# PS2020 cooling parameters +PS2020Cooling: + dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 + H_reion_z: 11.5 + H_reion_eV_p_H: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_eV_p_H: 2.0 + delta_logTEOS_subgrid_properties: 0.3 + rapid_cooling_threshold: 0.333333 + +# EAGLE star formation parameters +EAGLEStarFormation: + SF_threshold: Zdep # Zdep (Schaye 2004) or Subgrid + SF_model: PressureLaw # PressureLaw (Schaye et al. 2008) or SchmidtLaw + KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr. + KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law. + min_over_density: 100.0 # The over-density above which star-formation is allowed. + KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3. + KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold. + EOS_entropy_margin_dex: 0.3 # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form. + threshold_norm_H_p_cm3: 0.1 # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3. + threshold_Z0: 0.002 # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation. + threshold_slope: -0.64 # When using Z-based SF threshold, slope of the metal-dependant star formation threshold + threshold_max_density_H_p_cm3: 10.0 # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3. + threshold_temperature1_K: 1000 # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming. + threshold_temperature2_K: 31622 # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit. + threshold_number_density_H_p_cm3: 10 # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit. + +# Parameters for the EAGLE "equation of state" +EAGLEEntropyFloor: + Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in. + Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin. + Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor + Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in. + Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. + Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor + +# EAGLE feedback model +EAGLEFeedback: + use_SNII_feedback: 1 # Global switch for SNII thermal (stochastic) feedback. + use_SNIa_feedback: 1 # Global switch for SNIa thermal (continuous) feedback. + use_AGB_enrichment: 1 # Global switch for enrichement from AGB stars. + use_SNII_enrichment: 1 # Global switch for enrichement from SNII stars. + use_SNIa_enrichment: 1 # Global switch for enrichement from SNIa stars. + filename: ./yieldtables/ # Path to the directory containing the EAGLE yield tables. + IMF_min_mass_Msun: 0.1 # Minimal stellar mass considered for the Chabrier IMF in solar masses. + IMF_max_mass_Msun: 100.0 # Maximal stellar mass considered for the Chabrier IMF in solar masses. + # IMF configuration block (optional). If omitted, defaults to Chabrier 2003. + IMF_Model: "Chabrier" # Allowed: Chabrier | Kroupa | Salpeter | Custom + # Optional overrides for Chabrier or to customize: + IMF_HighMassSlope: 2.3 # alpha for m > PivotMass + IMF_PivotMass: 1.0 # Msun; break between low/high-mass regimes + IMF_ChabrierMc: 0.079 # Msun; characteristic mass of lognormal + IMF_ChabrierSigmaLog10: 0.69 # dispersion in base-10 log mass + # Kroupa example: + # IMF_Model: "Kroupa" + # IMF_LowMassSlope: 1.3 # alpha below PivotMass + # IMF_HighMassSlope: 2.3 # alpha above PivotMass + # IMF_PivotMass: 0.5 # Msun + # Salpeter example: + # IMF_Model: "Salpeter" + # IMF_HighMassSlope: 2.35 # single power-law slope + # Custom broken power-law example: + # IMF_Model: "Custom" + # IMF_LowMassSlope: 1.0 + # IMF_HighMassSlope: 2.0 + # IMF_PivotMass: 0.8 + SNII_min_mass_Msun: 6.0 # Minimal mass considered for SNII stars in solar masses. + SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII stars in solar masses. + SNII_feedback_model: Random # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity + SNII_sampled_delay: 0 # Sample the SNII lifetimes to do feedback. + SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event when not sampling. + SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin. + SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs. + SNII_energy_fraction_function: EAGLE # Type of functional form to use for scaling the energy fraction with density and metallicity ('EAGLE', 'Separable', or 'Independent'). + SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event. + SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event. + SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction). + SNII_energy_fraction_n_0_H_p_cm3: 1.4588 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3. + SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction. + SNII_energy_fraction_n_n: 0.8686 # Power-law for the birth density dependance of the SNII energy fraction. + SNII_energy_fraction_use_birth_density: 1 # Are we using the density at birth to compute f_E or at feedback time? + SNII_energy_fraction_use_birth_metallicity: 1 # Are we using the metallicity at birth to compuote f_E or at feedback time? + SNIa_DTD: PowerLaw # Functional form of the SNIa delay time distribution. + SNIa_DTD_delay_Gyr: 0.04 # Stellar age after which SNIa start in Gyr (40 Myr corresponds to stars ~ 8 Msun). + SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1). + SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs. + AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s. + stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled. + stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold. + SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel. + SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel. + SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel. + SNII_yield_factor_Nitrogen: 1.0 # (Optional) Correction factor to apply to the Nitrogen yield from the SNII channel. + SNII_yield_factor_Oxygen: 1.0 # (Optional) Correction factor to apply to the Oxygen yield from the SNII channel. + SNII_yield_factor_Neon: 1.0 # (Optional) Correction factor to apply to the Neon yield from the SNII channel. + SNII_yield_factor_Magnesium: 2.0 # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel. + SNII_yield_factor_Silicon: 1.0 # (Optional) Correction factor to apply to the Silicon yield from the SNII channel. + SNII_yield_factor_Iron: 0.5 # (Optional) Correction factor to apply to the Iron yield from the SNII channel. + +# EAGLE AGN model +EAGLEAGN: + subgrid_seed_mass_Msun: 1.5e5 # 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: 1 # 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 + 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. + radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated. + 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.15 # Fraction of the radiated energy that couples to the gas in feedback events. + AGN_feedback_model: Random # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity + AGN_use_deterministic_feedback: 0 # Deterministic (reservoir) [1] or stochastic [0] AGN feedback? + use_variable_delta_T: 0 # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0]. + AGN_with_locally_adaptive_delta_T: 0 # 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_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: 2e8 # Maximal BH mass considered for BH repositioning in solar masses. + max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length. + with_reposition_velocity_threshold: 1 # 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_ratio 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_ratio 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. + use_subgrid_mass_from_ics: 0 # Use the dynamical mass as the subgrid mass since we don't have subgrid masses in the ICs. diff --git a/src/feedback/EAGLE/enrichment.h b/src/feedback/EAGLE/enrichment.h index c680784439..c7e8bfc13b 100644 --- a/src/feedback/EAGLE/enrichment.h +++ b/src/feedback/EAGLE/enrichment.h @@ -689,7 +689,21 @@ void zero_yield_table_pointers(struct yield_table* table) { */ void feedback_restore_tables(struct feedback_props* fp) { - init_imf(fp); + /* Rebuild IMF according to stored options */ + { + struct eagle_imf_options opts; + opts.model = (fp->imf_model_code == 1) ? eagle_imf_model_kroupa + : (fp->imf_model_code == 2) ? eagle_imf_model_salpeter + : (fp->imf_model_code == 3) ? eagle_imf_model_custom + : eagle_imf_model_chabrier; + opts.high_mass_slope = fp->imf_high_mass_slope; + opts.low_mass_slope = fp->imf_low_mass_slope; + opts.pivot_mass_msun = fp->imf_pivot_mass_msun; + opts.chabrier_m_c_msun = fp->imf_chabrier_m_c_msun; + opts.chabrier_sigma_log10 = fp->imf_chabrier_sigma_log10; + + init_imf_from_options(fp, &opts); + } /* Allocate yield tables */ allocate_yield_tables(fp); diff --git a/src/feedback/EAGLE/imf.h b/src/feedback/EAGLE/imf.h index 564734bbfd..222f6a26b0 100644 --- a/src/feedback/EAGLE/imf.h +++ b/src/feedback/EAGLE/imf.h @@ -191,73 +191,257 @@ INLINE static double integrate_imf( } /** - * @brief Allocate space for IMF table and compute values to populate this - * table. - * - * @param feedback_props #feedback_props data structure + * @brief Supported IMF model choices when reading from YAML or options. */ -INLINE static void init_imf(struct feedback_props *feedback_props) { +enum eagle_imf_model { + eagle_imf_model_chabrier = 0, + eagle_imf_model_kroupa = 1, + eagle_imf_model_salpeter = 2, + eagle_imf_model_custom = 3 +} __attribute__((packed)); - /* Compute size of mass bins in log10 space */ - const double imf_log10_mass_bin_size = - (feedback_props->log10_imf_max_mass_msun - - feedback_props->log10_imf_min_mass_msun) / - (double)(eagle_feedback_N_imf_bins - 1); +/** + * @brief Options describing an IMF. Any value <= 0 is treated as "unset" and + * falls back to the model defaults. All units are in Msun and dimensionless + * slopes (phi ~ m^-alpha). + */ +struct eagle_imf_options { + enum eagle_imf_model model; /* pivot */ + double low_mass_slope; /*log10_imf_max_mass_msun - + feedback_props->log10_imf_min_mass_msun) / + (double)(eagle_feedback_N_imf_bins - 1); - /* Allocate IMF array */ if (swift_memalign("imf-tables", (void **)&feedback_props->imf, SWIFT_STRUCT_ALIGNMENT, eagle_feedback_N_imf_bins * sizeof(double)) != 0) error("Failed to allocate IMF bins table"); - /* Allocate array to store IMF mass bins */ if (swift_memalign("imf-tables", (void **)&feedback_props->imf_mass_bin, SWIFT_STRUCT_ALIGNMENT, eagle_feedback_N_imf_bins * sizeof(double)) != 0) error("Failed to allocate IMF bins table"); - /* Allocate array to store IMF mass bins in log10 space */ if (swift_memalign("imf-tables", (void **)&feedback_props->imf_mass_bin_log10, SWIFT_STRUCT_ALIGNMENT, eagle_feedback_N_imf_bins * sizeof(double)) != 0) error("Failed to allocate IMF bins table"); - /* Set IMF from Chabrier 2003, PASP, 115, 763 - * Eq. 17 with values from table 1 */ - for (int i = 0; i < eagle_feedback_N_imf_bins; i++) { + *imf_log10_mass_bin_size = dlog10; +} - /* Logarithmically-spaced bins in units of solar masses */ - const double log10_mass_msun = - feedback_props->log10_imf_min_mass_msun + i * imf_log10_mass_bin_size; +/** + * @brief Internal helper: normalize IMF so that ∫ m φ(m) dm = 1 across [m_min, + * m_max]. + */ +INLINE static void eagle_imf_normalize(struct feedback_props *feedback_props) { + const float norm = integrate_imf(feedback_props->log10_imf_min_mass_msun, + feedback_props->log10_imf_max_mass_msun, + eagle_imf_integration_mass_weight, + /* yields */ NULL, feedback_props); + for (int i = 0; i < eagle_feedback_N_imf_bins; i++) + feedback_props->imf[i] /= norm; +} - const double mass_msun = exp10(log10_mass_msun); +/** + * @brief Initialize a Chabrier (2003) IMF with optional overrides. + * - Low-mass: lognormal with m_c and sigma in log10. + * - High-mass: power-law with slope alpha_high, matched continuously at pivot. + */ +INLINE static void init_imf_chabrier(struct feedback_props *feedback_props, + double alpha_high_opt, + double pivot_mass_opt, double m_c_opt, + double sigma_log10_opt) { + const double alpha_high = + (alpha_high_opt > 0.0) ? alpha_high_opt : 2.3; /* default Chabrier */ + const double pivot_mass = (pivot_mass_opt > 0.0) ? pivot_mass_opt : 1.0; + const double m_c = (m_c_opt > 0.0) ? m_c_opt : 0.079; + const double sig_log10 = (sigma_log10_opt > 0.0) ? sigma_log10_opt : 0.69; + + double dlog10; + eagle_imf_allocate_arrays(feedback_props, &dlog10); + + /* Low-mass lognormal normalization constant used historically in SWIFT */ + const double log10_mc = log10(m_c); + + /* Value of the lognormal at the pivot to match continuity */ + const double log10_pivot = log10(pivot_mass); + const double phi_low_at_pivot = + 0.852464 * + exp((log10_pivot - log10_mc) * (log10_pivot - log10_mc) / + (-2.0 * sig_log10 * sig_log10)) / + pivot_mass; + + /* High-mass normalization to ensure continuity at pivot */ + const double A_high = phi_low_at_pivot * pow(pivot_mass, alpha_high); - feedback_props->imf_mass_bin[i] = mass_msun; - feedback_props->imf_mass_bin_log10[i] = log10_mass_msun; + for (int i = 0; i < eagle_feedback_N_imf_bins; i++) { + const double log10_m = feedback_props->log10_imf_min_mass_msun + i * dlog10; + const double m = exp10(log10_m); - if (mass_msun > 1.0) { + feedback_props->imf_mass_bin[i] = m; + feedback_props->imf_mass_bin_log10[i] = log10_m; - /* High-mass end */ - feedback_props->imf[i] = 0.237912 * pow(mass_msun, -2.3); + if (m > pivot_mass) { + feedback_props->imf[i] = A_high * pow(m, -alpha_high); } else { + feedback_props->imf[i] = 0.852464 * + exp((log10_m - log10_mc) * (log10_m - log10_mc) / + (-2.0 * sig_log10 * sig_log10)) / + m; + } + } + + eagle_imf_normalize(feedback_props); +} - /* Low-mass end */ - feedback_props->imf[i] = - 0.852464 * - exp((log10_mass_msun - log10(0.079)) * - (log10_mass_msun - log10(0.079)) / (-2.0 * 0.69 * 0.69)) / - mass_msun; +/** + * @brief Initialize a Kroupa (2001) broken power-law IMF. + * Defaults: alpha_low=1.3 below pivot=0.5 Msun, alpha_high=2.3 above. + */ +INLINE static void init_imf_kroupa(struct feedback_props *feedback_props, + double alpha_low_opt, double alpha_high_opt, + double pivot_mass_opt) { + const double alpha_low = (alpha_low_opt > 0.0) ? alpha_low_opt : 1.3; + const double alpha_high = (alpha_high_opt > 0.0) ? alpha_high_opt : 2.3; + const double pivot_mass = (pivot_mass_opt > 0.0) ? pivot_mass_opt : 0.5; + + double dlog10; + eagle_imf_allocate_arrays(feedback_props, &dlog10); + + /* Choose A_low arbitrarily; continuity determines A_high; mass-normalization + * comes later. */ + const double A_low = 1.0; + const double A_high = A_low * pow(pivot_mass, alpha_high - alpha_low); + + for (int i = 0; i < eagle_feedback_N_imf_bins; i++) { + const double log10_m = feedback_props->log10_imf_min_mass_msun + i * dlog10; + const double m = exp10(log10_m); + + feedback_props->imf_mass_bin[i] = m; + feedback_props->imf_mass_bin_log10[i] = log10_m; + + if (m > pivot_mass) { + feedback_props->imf[i] = A_high * pow(m, -alpha_high); + } else { + feedback_props->imf[i] = A_low * pow(m, -alpha_low); } } - /* Normalize the IMF */ - const float norm = integrate_imf(feedback_props->log10_imf_min_mass_msun, - feedback_props->log10_imf_max_mass_msun, - eagle_imf_integration_mass_weight, - /*(stellar_yields=)*/ NULL, feedback_props); + eagle_imf_normalize(feedback_props); +} - for (int i = 0; i < eagle_feedback_N_imf_bins; i++) - feedback_props->imf[i] /= norm; +/** + * @brief Initialize a Salpeter (1955) single power-law IMF. + * Default slope alpha=2.35. + */ +INLINE static void init_imf_salpeter(struct feedback_props *feedback_props, + double alpha_opt) { + const double alpha = (alpha_opt > 0.0) ? alpha_opt : 2.35; + + double dlog10; + eagle_imf_allocate_arrays(feedback_props, &dlog10); + + const double A = 1.0; /* arbitrary; normalize by mass afterwards */ + + for (int i = 0; i < eagle_feedback_N_imf_bins; i++) { + const double log10_m = feedback_props->log10_imf_min_mass_msun + i * dlog10; + const double m = exp10(log10_m); + + feedback_props->imf_mass_bin[i] = m; + feedback_props->imf_mass_bin_log10[i] = log10_m; + + feedback_props->imf[i] = A * pow(m, -alpha); + } + + eagle_imf_normalize(feedback_props); +} + +/** + * @brief Initialize a custom IMF. Behavior: + * - If both low_mass_slope and high_mass_slope are set (>0) and pivot set (>0): + * broken power-law with continuity at pivot. + * - Else if only high_mass_slope is set (>0): Chabrier low-mass lognormal + + * that high-mass slope at pivot=1 Msun. + * - Else: falls back to classic Chabrier. + */ +INLINE static void init_imf_custom(struct feedback_props *feedback_props, + const struct eagle_imf_options *opt) { + const double alpha_low = + (opt && opt->low_mass_slope > 0.0) ? opt->low_mass_slope : -1.0; + const double alpha_high = + (opt && opt->high_mass_slope > 0.0) ? opt->high_mass_slope : -1.0; + const double pivot_mass = + (opt && opt->pivot_mass_msun > 0.0) ? opt->pivot_mass_msun : -1.0; + + if (alpha_low > 0.0 && alpha_high > 0.0 && pivot_mass > 0.0) { + init_imf_kroupa(feedback_props, alpha_low, alpha_high, pivot_mass); + } else if (alpha_high > 0.0) { + /* Chabrier-like with custom high-mass slope */ + init_imf_chabrier( + feedback_props, alpha_high, + (opt && opt->pivot_mass_msun > 0.0) ? opt->pivot_mass_msun : 1.0, + (opt && opt->chabrier_m_c_msun > 0.0) ? opt->chabrier_m_c_msun : 0.079, + (opt && opt->chabrier_sigma_log10 > 0.0) ? opt->chabrier_sigma_log10 + : 0.69); + } else { + init_imf_chabrier(feedback_props, -1.0, -1.0, -1.0, -1.0); + } +} + +/** + * @brief Initialize IMF from options (e.g. read from YAML). Defaults to + * Chabrier when options are NULL or fields are unset. + */ +INLINE static void init_imf_from_options(struct feedback_props *feedback_props, + const struct eagle_imf_options *opt) { + enum eagle_imf_model model = (opt) ? opt->model : eagle_imf_model_chabrier; + switch (model) { + case eagle_imf_model_chabrier: + init_imf_chabrier(feedback_props, (opt ? opt->high_mass_slope : -1.0), + (opt ? opt->pivot_mass_msun : -1.0), + (opt ? opt->chabrier_m_c_msun : -1.0), + (opt ? opt->chabrier_sigma_log10 : -1.0)); + break; + case eagle_imf_model_kroupa: + init_imf_kroupa(feedback_props, (opt ? opt->low_mass_slope : -1.0), + (opt ? opt->high_mass_slope : -1.0), + (opt ? opt->pivot_mass_msun : -1.0)); + break; + case eagle_imf_model_salpeter: + init_imf_salpeter(feedback_props, (opt ? opt->high_mass_slope : -1.0)); + break; + case eagle_imf_model_custom: + default: + init_imf_custom(feedback_props, opt); + break; + } +} + +/** + * @brief Backward-compatible wrapper: initialize IMF with default (Chabrier) + * high-mass slope. + */ +INLINE static void init_imf(struct feedback_props *feedback_props) { + init_imf_chabrier(feedback_props, /*alpha_high=*/-1.0, /*pivot=*/-1.0, + /*m_c=*/-1.0, /*sigma_log10=*/-1.0); } /** diff --git a/src/feedback/EAGLE_kinetic/feedback.c b/src/feedback/EAGLE_kinetic/feedback.c index 193ec95820..54d4475bde 100644 --- a/src/feedback/EAGLE_kinetic/feedback.c +++ b/src/feedback/EAGLE_kinetic/feedback.c @@ -448,6 +448,42 @@ void feedback_props_init(struct feedback_props* fp, fp->log10_imf_max_mass_msun = log10(fp->imf_max_mass_msun); fp->log10_imf_min_mass_msun = log10(fp->imf_min_mass_msun); + /* Optional IMF configuration (default: Chabrier 2003) */ + { + fp->imf_model_code = 0; + fp->imf_high_mass_slope = -1.0; + fp->imf_low_mass_slope = -1.0; + fp->imf_pivot_mass_msun = -1.0; + fp->imf_chabrier_m_c_msun = -1.0; + fp->imf_chabrier_sigma_log10 = -1.0; + + char model_str[32] = {0}; + parser_get_opt_param_string(params, "EAGLEFeedback:IMF_Model", model_str, + "Chabrier"); + if (strcmp(model_str, "Chabrier") == 0) + fp->imf_model_code = 0; + else if (strcmp(model_str, "Kroupa") == 0) + fp->imf_model_code = 1; + else if (strcmp(model_str, "Salpeter") == 0) + fp->imf_model_code = 2; + else if (strcmp(model_str, "Custom") == 0) + fp->imf_model_code = 3; + else + error("Invalid EAGLEFeedback:IMF_Model '%s'", model_str); + + fp->imf_high_mass_slope = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_HighMassSlope", fp->imf_high_mass_slope); + fp->imf_low_mass_slope = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_LowMassSlope", fp->imf_low_mass_slope); + fp->imf_pivot_mass_msun = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_PivotMass", fp->imf_pivot_mass_msun); + fp->imf_chabrier_m_c_msun = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_ChabrierMc", fp->imf_chabrier_m_c_msun); + fp->imf_chabrier_sigma_log10 = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_ChabrierSigmaLog10", + fp->imf_chabrier_sigma_log10); + } + /* Properties of the SNII energy feedback model ------------------------- */ /* Are we sampling the SNII lifetimes for feedback or using a fixed delay? */ @@ -628,8 +664,20 @@ void feedback_props_init(struct feedback_props* fp, (X_H / m_p) * units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY); /* Initialise the IMF ------------------------------------------------- */ - - init_imf(fp); + { + struct eagle_imf_options opts; + opts.model = (fp->imf_model_code == 1) ? eagle_imf_model_kroupa + : (fp->imf_model_code == 2) ? eagle_imf_model_salpeter + : (fp->imf_model_code == 3) ? eagle_imf_model_custom + : eagle_imf_model_chabrier; + opts.high_mass_slope = fp->imf_high_mass_slope; + opts.low_mass_slope = fp->imf_low_mass_slope; + opts.pivot_mass_msun = fp->imf_pivot_mass_msun; + opts.chabrier_m_c_msun = fp->imf_chabrier_m_c_msun; + opts.chabrier_sigma_log10 = fp->imf_chabrier_sigma_log10; + + init_imf_from_options(fp, &opts); + } /* Calculate number of type II SN per unit solar mass based on our choice * of IMF and integration limits for type II SNe. diff --git a/src/feedback/EAGLE_kinetic/feedback_properties.h b/src/feedback/EAGLE_kinetic/feedback_properties.h index b5a9de747d..a6bd630ce0 100644 --- a/src/feedback/EAGLE_kinetic/feedback_properties.h +++ b/src/feedback/EAGLE_kinetic/feedback_properties.h @@ -214,6 +214,20 @@ struct feedback_props { */ double log10_imf_max_mass_msun; + /* ---- IMF model selection (for restart/restore) ---- */ + /*! Encoded IMF model: 0=Chabrier, 1=Kroupa, 2=Salpeter, 3=Custom */ + int imf_model_code; + /*! Optional high-mass slope (>0 means set) */ + double imf_high_mass_slope; + /*! Optional low-mass slope (>0 means set) */ + double imf_low_mass_slope; + /*! Optional pivot/break mass in Msun (>0 means set) */ + double imf_pivot_mass_msun; + /*! Optional Chabrier lognormal characteristic mass (Msun) */ + double imf_chabrier_m_c_msun; + /*! Optional Chabrier lognormal dispersion in log10 */ + double imf_chabrier_sigma_log10; + /* ------------ SNe feedback properties ------------ */ /*! Minimal stellar mass considered for SNII feedback (in solar masses) */ diff --git a/src/feedback/EAGLE_thermal/feedback.c b/src/feedback/EAGLE_thermal/feedback.c index fe799296d3..5af168bb93 100644 --- a/src/feedback/EAGLE_thermal/feedback.c +++ b/src/feedback/EAGLE_thermal/feedback.c @@ -451,7 +451,7 @@ void feedback_props_init(struct feedback_props* fp, fp->imf_max_mass_msun = parser_get_param_double(params, "EAGLEFeedback:IMF_max_mass_Msun"); fp->imf_min_mass_msun = - parser_get_param_double(params, "EAGLEFeedback:IMF_min_mass_Msun"); + parser_get_param_double(params, "EAGLEFeedback:IMF_min_mass_MSun"); /* Check that it makes sense. */ if (fp->imf_max_mass_msun < fp->imf_min_mass_msun) { @@ -461,6 +461,44 @@ void feedback_props_init(struct feedback_props* fp, fp->log10_imf_max_mass_msun = log10(fp->imf_max_mass_msun); fp->log10_imf_min_mass_msun = log10(fp->imf_min_mass_msun); + /* Optional IMF configuration (default: Chabrier 2003) */ + { + /* Initialize defaults */ + fp->imf_model_code = 0; + fp->imf_high_mass_slope = -1.0; + fp->imf_low_mass_slope = -1.0; + fp->imf_pivot_mass_msun = -1.0; + fp->imf_chabrier_m_c_msun = -1.0; + fp->imf_chabrier_sigma_log10 = -1.0; + + char model_str[32] = {0}; + parser_get_opt_param_string(params, "EAGLEFeedback:IMF_Model", model_str, + "Chabrier"); + + if (strcmp(model_str, "Chabrier") == 0) + fp->imf_model_code = 0; + else if (strcmp(model_str, "Kroupa") == 0) + fp->imf_model_code = 1; + else if (strcmp(model_str, "Salpeter") == 0) + fp->imf_model_code = 2; + else if (strcmp(model_str, "Custom") == 0) + fp->imf_model_code = 3; + else + error("Invalid EAGLEFeedback:IMF_Model '%s'", model_str); + + fp->imf_high_mass_slope = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_HighMassSlope", fp->imf_high_mass_slope); + fp->imf_low_mass_slope = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_LowMassSlope", fp->imf_low_mass_slope); + fp->imf_pivot_mass_msun = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_PivotMass", fp->imf_pivot_mass_msun); + fp->imf_chabrier_m_c_msun = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_ChabrierMc", fp->imf_chabrier_m_c_msun); + fp->imf_chabrier_sigma_log10 = parser_get_opt_param_double( + params, "EAGLEFeedback:IMF_ChabrierSigmaLog10", + fp->imf_chabrier_sigma_log10); + } + /* Properties of the SNII energy feedback model ------------------------- */ char model[64]; @@ -711,8 +749,20 @@ void feedback_props_init(struct feedback_props* fp, (X_H / m_p) * units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY); /* Initialise the IMF ------------------------------------------------- */ - - init_imf(fp); + { + struct eagle_imf_options opts; + opts.model = (fp->imf_model_code == 1) ? eagle_imf_model_kroupa + : (fp->imf_model_code == 2) ? eagle_imf_model_salpeter + : (fp->imf_model_code == 3) ? eagle_imf_model_custom + : eagle_imf_model_chabrier; + opts.high_mass_slope = fp->imf_high_mass_slope; + opts.low_mass_slope = fp->imf_low_mass_slope; + opts.pivot_mass_msun = fp->imf_pivot_mass_msun; + opts.chabrier_m_c_msun = fp->imf_chabrier_m_c_msun; + opts.chabrier_sigma_log10 = fp->imf_chabrier_sigma_log10; + + init_imf_from_options(fp, &opts); + } /* Calculate number of type II SN per unit solar mass based on our choice * of IMF and integration limits for type II SNe. diff --git a/src/feedback/EAGLE_thermal/feedback_properties.h b/src/feedback/EAGLE_thermal/feedback_properties.h index 5e1e38b11f..a89173c403 100644 --- a/src/feedback/EAGLE_thermal/feedback_properties.h +++ b/src/feedback/EAGLE_thermal/feedback_properties.h @@ -233,6 +233,20 @@ struct feedback_props { */ double log10_imf_max_mass_msun; + /* ---- IMF model selection ---- */ + /*! Encoded IMF model: 0=Chabrier, 1=Kroupa, 2=Salpeter, 3=Custom */ + int imf_model_code; + /*! Optional high-mass slope (>0 means set) */ + double imf_high_mass_slope; + /*! Optional low-mass slope (>0 means set) */ + double imf_low_mass_slope; + /*! Optional pivot/break mass in Msun (>0 means set) */ + double imf_pivot_mass_msun; + /*! Optional Chabrier lognormal characteristic mass (Msun) */ + double imf_chabrier_m_c_msun; + /*! Optional Chabrier lognormal dispersion in log10 */ + double imf_chabrier_sigma_log10; + /* ------------ SNe feedback properties ------------ */ /*! SNII feedback model: random, isotropic or minimum distance */ diff --git a/tools/imf_ks_normalisation.py b/tools/imf_ks_normalisation.py new file mode 100644 index 0000000000..f53db8c6d8 --- /dev/null +++ b/tools/imf_ks_normalisation.py @@ -0,0 +1,53 @@ +import numpy as np + +import fsps + + +fiducial_norm_constant = 1.65 # value used in Schaye+15 +fiducial_normalisation = 1.515e-4 * fiducial_norm_constant + +ages = np.linspace(0, 10, 100) + +sp = fsps.StellarPopulation( + imf_type=0, # Salpeter+55 + zcontinuous=1, + sfh=3, + logzsol=0.0, + add_agb_dust_model=False, + use_wr_spectra=False, + imf_upper_limit=100, + # imf_lower_limit=0.1, +) + +sp.set_tabular_sfh( + age=ages, + sfr=np.ones(ages.shape), +) + +salpeter_mag = sp.get_mags(tage=0.1, bands=['galex_fuv']) + +sp = fsps.StellarPopulation( + imf_type=1, # Chabrier+03 + zcontinuous=1, + sfh=3, + logzsol=0.0, + add_agb_dust_model=False, + use_wr_spectra=False, + imf_upper_limit=100, + # imf_lower_limit=0.1, +) + +sp.set_tabular_sfh( + age=ages, + sfr=np.ones(ages.shape), +) + +compare_mag = sp.get_mags(tage=0.1, bands=['galex_fuv']) + +new_norm_constant = salpeter_mag - compare_mag +new_normalisation = fiducial_normalisation * new_norm_constant + +print("Ratio (Custom / Salpeter):", 1 / new_norm_constant) +print("New normalisation:", new_normalisation) + +