From 924aceb69bc88f384d27f5f96f93c967fdc16640 Mon Sep 17 00:00:00 2001 From: christopherlovell Date: Thu, 21 Aug 2025 15:27:12 +0100 Subject: [PATCH 1/7] imf choice, and example submission script --- examples/EAGLE_low_z/EAGLE_6/imf_example.yml | 275 +++++++++++++++++++ src/feedback/EAGLE/imf.h | 251 ++++++++++++++--- 2 files changed, 491 insertions(+), 35 deletions(-) create mode 100644 examples/EAGLE_low_z/EAGLE_6/imf_example.yml diff --git a/examples/EAGLE_low_z/EAGLE_6/imf_example.yml b/examples/EAGLE_low_z/EAGLE_6/imf_example.yml new file mode 100644 index 0000000000..22564400af --- /dev/null +++ b/examples/EAGLE_low_z/EAGLE_6/imf_example.yml @@ -0,0 +1,275 @@ +# This file is identical to examples/EAGLE_low_z/EAGLE_25/eagle_25.yml +# with an added EAGLEFeedback.IMF block to demonstrate IMF configuration options +# in the EAGLE model. + +MetaData: + run_name: EAGLE-L0025N0376-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: + 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) + 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: 256 + 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] # Use DM as the primary FOF linking type + attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./EAGLE_ICs_25.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 + +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. + + # Added: IMF configuration block (optional). If omitted, defaults to Chabrier. + IMF: + Model: "Chabrier" # Allowed: Chabrier | Kroupa | Salpeter | Custom + # Optional overrides for Chabrier or to customize: + HighMassSlope: 2.3 # alpha for m > PivotMass + PivotMass: 1.0 # Msun; break between low/high-mass regimes + ChabrierMc: 0.079 # Msun; characteristic mass of lognormal + ChabrierSigmaLog10: 0.69 # dispersion in base-10 log mass + # Kroupa example: + # Model: "Kroupa" + # LowMassSlope: 1.3 # alpha below PivotMass + # HighMassSlope: 2.3 # alpha above PivotMass + # PivotMass: 0.5 # Msun + # Salpeter example: + # Model: "Salpeter" + # HighMassSlope: 2.35 # single power-law slope + # Custom broken power-law example: + # Model: "Custom" + # LowMassSlope: 1.0 + # HighMassSlope: 2.0 + # 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/imf.h b/src/feedback/EAGLE/imf.h index 564734bbfd..2a3e15986d 100644 --- a/src/feedback/EAGLE/imf.h +++ b/src/feedback/EAGLE/imf.h @@ -190,74 +190,255 @@ INLINE static double integrate_imf( return result * imf_log10_mass_bin_size * M_LN10; } +/* ================================================================ + * New: IMF models and options-based initialization (YAML-friendly) + * ================================================================ */ + /** - * @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 = +/** + * @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); +} + +/** + * @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); - /* 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; + 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); +} + +/** + * @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; + } +} - for (int i = 0; i < eagle_feedback_N_imf_bins; i++) - feedback_props->imf[i] /= norm; +/** + * @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); } /** From 1fac35da1cddc518b907996386219adc1c3967dd Mon Sep 17 00:00:00 2001 From: christopherlovell Date: Thu, 21 Aug 2025 17:48:28 +0100 Subject: [PATCH 2/7] add FSPS tool to calculate KS normalisation for alternative IMF --- tools/imf_ks_normalisation.py | 45 +++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 tools/imf_ks_normalisation.py diff --git a/tools/imf_ks_normalisation.py b/tools/imf_ks_normalisation.py new file mode 100644 index 0000000000..7471a6d727 --- /dev/null +++ b/tools/imf_ks_normalisation.py @@ -0,0 +1,45 @@ +import numpy as np + +import fsps + +sp = fsps.StellarPopulation( + imf_type=0, # Salpeter+55 + compute_vega_mags=False, + 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=np.array([0,1,2]), + sfr=np.array([1,1,1]), +) + +salpeter_mag = sp.get_mags(tage=0.1, bands=['galex_fuv']) + +sp = fsps.StellarPopulation( + imf_type=1, # Chabrier+03 + compute_vega_mags=False, + 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=np.array([0,1,2]), + sfr=np.array([1,1,1]), +) + +compare_mag = sp.get_mags(tage=0.1, bands=['galex_fuv']) + +print("Ratio (Custom / Salpeter):", 1 / (salpeter_mag - compare_mag)) + + From 2ba2a2980542878c38b86126714af191113b8d21 Mon Sep 17 00:00:00 2001 From: christopherlovell Date: Thu, 21 Aug 2025 17:54:50 +0100 Subject: [PATCH 3/7] update fsps parameters --- tools/imf_ks_normalisation.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tools/imf_ks_normalisation.py b/tools/imf_ks_normalisation.py index 7471a6d727..736dd86874 100644 --- a/tools/imf_ks_normalisation.py +++ b/tools/imf_ks_normalisation.py @@ -2,40 +2,40 @@ import fsps +ages = np.linspace(0, 10, 100) + sp = fsps.StellarPopulation( imf_type=0, # Salpeter+55 - compute_vega_mags=False, 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, + # imf_lower_limit=0.1, ) sp.set_tabular_sfh( - age=np.array([0,1,2]), - sfr=np.array([1,1,1]), + 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 - compute_vega_mags=False, 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, + # imf_lower_limit=0.1, ) sp.set_tabular_sfh( - age=np.array([0,1,2]), - sfr=np.array([1,1,1]), + age=ages, + sfr=np.ones(ages.shape), ) compare_mag = sp.get_mags(tage=0.1, bands=['galex_fuv']) From d48ad3a8aa0e15c85e3f80b66ba3403035b18f9a Mon Sep 17 00:00:00 2001 From: christopherlovell Date: Thu, 21 Aug 2025 17:58:16 +0100 Subject: [PATCH 4/7] imf variations updated in EAGLE feedback --- src/feedback/EAGLE/enrichment.h | 19 +++++- src/feedback/EAGLE/imf.h | 4 -- src/feedback/EAGLE_kinetic/feedback.c | 55 ++++++++++++++++- .../EAGLE_kinetic/feedback_properties.h | 14 +++++ src/feedback/EAGLE_thermal/feedback.c | 59 ++++++++++++++++++- .../EAGLE_thermal/feedback_properties.h | 14 +++++ 6 files changed, 155 insertions(+), 10 deletions(-) diff --git a/src/feedback/EAGLE/enrichment.h b/src/feedback/EAGLE/enrichment.h index c680784439..c24f73faf6 100644 --- a/src/feedback/EAGLE/enrichment.h +++ b/src/feedback/EAGLE/enrichment.h @@ -689,7 +689,24 @@ 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 2a3e15986d..564637dfc6 100644 --- a/src/feedback/EAGLE/imf.h +++ b/src/feedback/EAGLE/imf.h @@ -190,10 +190,6 @@ INLINE static double integrate_imf( return result * imf_log10_mass_bin_size * M_LN10; } -/* ================================================================ - * New: IMF models and options-based initialization (YAML-friendly) - * ================================================================ */ - /** * @brief Supported IMF model choices when reading from YAML or options. */ diff --git a/src/feedback/EAGLE_kinetic/feedback.c b/src/feedback/EAGLE_kinetic/feedback.c index 193ec95820..0875c74cb2 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,23 @@ 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..c83286f1df 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,23 @@ 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 */ From d3a1668d57f2dd53922c4e5452a20e75a7d63647 Mon Sep 17 00:00:00 2001 From: christopherlovell Date: Thu, 21 Aug 2025 17:58:40 +0100 Subject: [PATCH 5/7] rename example file --- .../{imf_example.yml => imf_example_6.yml} | 84 +++++++++---------- 1 file changed, 40 insertions(+), 44 deletions(-) rename examples/EAGLE_low_z/EAGLE_6/{imf_example.yml => imf_example_6.yml} (92%) diff --git a/examples/EAGLE_low_z/EAGLE_6/imf_example.yml b/examples/EAGLE_low_z/EAGLE_6/imf_example_6.yml similarity index 92% rename from examples/EAGLE_low_z/EAGLE_6/imf_example.yml rename to examples/EAGLE_low_z/EAGLE_6/imf_example_6.yml index 22564400af..75a737c535 100644 --- a/examples/EAGLE_low_z/EAGLE_6/imf_example.yml +++ b/examples/EAGLE_low_z/EAGLE_6/imf_example_6.yml @@ -1,9 +1,5 @@ -# This file is identical to examples/EAGLE_low_z/EAGLE_25/eagle_25.yml -# with an added EAGLEFeedback.IMF block to demonstrate IMF configuration options -# in the EAGLE model. - MetaData: - run_name: EAGLE-L0025N0376-Ref + run_name: EAGLE-L0006N0094-Ref # Define the system of units to use internally. InternalUnitSystem: @@ -15,7 +11,7 @@ InternalUnitSystem: # Structure finding options StructureFinding: - config_file_name: stf_input.cfg # Name of the STF config file. + 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). @@ -29,23 +25,26 @@ Cosmology: 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 @@ -61,13 +60,12 @@ Gravity: MAC: adaptive epsilon_fmm: 0.001 theta_cr: 0.7 # Opening angle (Multipole acceptance criterion) - mesh_side_length: 256 + 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). @@ -95,16 +93,17 @@ FOF: 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] # Use DM as the primary FOF linking type - attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types - + 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_25.hdf5 # The file to read + 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 @@ -116,7 +115,7 @@ EAGLEChemistry: # Solar abundances 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/ @@ -136,7 +135,7 @@ PS2020Cooling: 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 @@ -175,32 +174,29 @@ EAGLEFeedback: 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. - - # Added: IMF configuration block (optional). If omitted, defaults to Chabrier. - IMF: - Model: "Chabrier" # Allowed: Chabrier | Kroupa | Salpeter | Custom - # Optional overrides for Chabrier or to customize: - HighMassSlope: 2.3 # alpha for m > PivotMass - PivotMass: 1.0 # Msun; break between low/high-mass regimes - ChabrierMc: 0.079 # Msun; characteristic mass of lognormal - ChabrierSigmaLog10: 0.69 # dispersion in base-10 log mass - # Kroupa example: - # Model: "Kroupa" - # LowMassSlope: 1.3 # alpha below PivotMass - # HighMassSlope: 2.3 # alpha above PivotMass - # PivotMass: 0.5 # Msun - # Salpeter example: - # Model: "Salpeter" - # HighMassSlope: 2.35 # single power-law slope - # Custom broken power-law example: - # Model: "Custom" - # LowMassSlope: 1.0 - # HighMassSlope: 2.0 - # 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. + 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. From 41581170d9590acd503a778ffbd960832b81b75b Mon Sep 17 00:00:00 2001 From: christopherlovell Date: Thu, 21 Aug 2025 18:28:57 +0100 Subject: [PATCH 6/7] add full normalisation values --- tools/imf_ks_normalisation.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/tools/imf_ks_normalisation.py b/tools/imf_ks_normalisation.py index 736dd86874..f53db8c6d8 100644 --- a/tools/imf_ks_normalisation.py +++ b/tools/imf_ks_normalisation.py @@ -2,6 +2,10 @@ 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( @@ -40,6 +44,10 @@ compare_mag = sp.get_mags(tage=0.1, bands=['galex_fuv']) -print("Ratio (Custom / Salpeter):", 1 / (salpeter_mag - compare_mag)) +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) From 982dbcc49869d39b6d00512650942989ce12c150 Mon Sep 17 00:00:00 2001 From: christopherlovell Date: Thu, 21 Aug 2025 18:29:30 +0100 Subject: [PATCH 7/7] format updates --- src/feedback/EAGLE/enrichment.h | 11 ++-- src/feedback/EAGLE/imf.h | 95 ++++++++++++++------------- src/feedback/EAGLE_kinetic/feedback.c | 11 ++-- src/feedback/EAGLE_thermal/feedback.c | 11 ++-- 4 files changed, 63 insertions(+), 65 deletions(-) diff --git a/src/feedback/EAGLE/enrichment.h b/src/feedback/EAGLE/enrichment.h index c24f73faf6..c7e8bfc13b 100644 --- a/src/feedback/EAGLE/enrichment.h +++ b/src/feedback/EAGLE/enrichment.h @@ -692,13 +692,10 @@ void feedback_restore_tables(struct feedback_props* 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.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; diff --git a/src/feedback/EAGLE/imf.h b/src/feedback/EAGLE/imf.h index 564637dfc6..222f6a26b0 100644 --- a/src/feedback/EAGLE/imf.h +++ b/src/feedback/EAGLE/imf.h @@ -197,7 +197,7 @@ enum eagle_imf_model { eagle_imf_model_chabrier = 0, eagle_imf_model_kroupa = 1, eagle_imf_model_salpeter = 2, - eagle_imf_model_custom = 3 + eagle_imf_model_custom = 3 } __attribute__((packed)); /** @@ -206,28 +206,28 @@ enum eagle_imf_model { * slopes (phi ~ m^-alpha). */ struct eagle_imf_options { - enum eagle_imf_model model; /* pivot */ - double low_mass_slope; /* pivot */ + double low_mass_slope; /*log10_imf_max_mass_msun - - feedback_props->log10_imf_min_mass_msun) / - (double)(eagle_feedback_N_imf_bins - 1); +INLINE static void eagle_imf_allocate_arrays( + struct feedback_props *feedback_props, double *imf_log10_mass_bin_size) { + const double dlog10 = (feedback_props->log10_imf_max_mass_msun - + feedback_props->log10_imf_min_mass_msun) / + (double)(eagle_feedback_N_imf_bins - 1); if (swift_memalign("imf-tables", (void **)&feedback_props->imf, SWIFT_STRUCT_ALIGNMENT, @@ -248,14 +248,16 @@ INLINE static void eagle_imf_allocate_arrays(struct feedback_props *feedback_pro } /** - * @brief Internal helper: normalize IMF so that ∫ m φ(m) dm = 1 across [m_min, m_max]. + * @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 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; } /** @@ -265,13 +267,13 @@ INLINE static void eagle_imf_normalize(struct feedback_props *feedback_props) { */ INLINE static void init_imf_chabrier(struct feedback_props *feedback_props, double alpha_high_opt, - double pivot_mass_opt, - double m_c_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 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; + 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); @@ -315,17 +317,17 @@ INLINE static void init_imf_chabrier(struct feedback_props *feedback_props, * 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 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_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. */ + /* 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); @@ -376,23 +378,29 @@ INLINE static void init_imf_salpeter(struct feedback_props *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 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; + 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); + 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); } @@ -407,15 +415,13 @@ INLINE static void init_imf_from_options(struct feedback_props *feedback_props, 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)); + 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), + 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; @@ -430,7 +436,8 @@ INLINE static void init_imf_from_options(struct feedback_props *feedback_props, } /** - * @brief Backward-compatible wrapper: initialize IMF with default (Chabrier) high-mass slope. + * @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, diff --git a/src/feedback/EAGLE_kinetic/feedback.c b/src/feedback/EAGLE_kinetic/feedback.c index 0875c74cb2..54d4475bde 100644 --- a/src/feedback/EAGLE_kinetic/feedback.c +++ b/src/feedback/EAGLE_kinetic/feedback.c @@ -666,13 +666,10 @@ void feedback_props_init(struct feedback_props* fp, /* Initialise the IMF ------------------------------------------------- */ { 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.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; diff --git a/src/feedback/EAGLE_thermal/feedback.c b/src/feedback/EAGLE_thermal/feedback.c index c83286f1df..5af168bb93 100644 --- a/src/feedback/EAGLE_thermal/feedback.c +++ b/src/feedback/EAGLE_thermal/feedback.c @@ -751,13 +751,10 @@ void feedback_props_init(struct feedback_props* fp, /* Initialise the IMF ------------------------------------------------- */ { 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.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;