From 865e509a3c380e06e4e128175ea22b126e458136 Mon Sep 17 00:00:00 2001 From: James Davies Date: Fri, 5 Apr 2024 14:40:29 +0200 Subject: [PATCH 1/5] make parameters for upper turnover, fix subsampler and test plots --- src/py21cmfast/inputs.py | 11 ++++ src/py21cmfast/src/21cmFAST.h | 1 + src/py21cmfast/src/Globals.h | 4 ++ src/py21cmfast/src/HaloBox.c | 82 ++++++++++++---------------- src/py21cmfast/src/interp_tables.c | 7 +-- src/py21cmfast/src/ps.c | 44 +++++++++++---- tests/test_c_interpolation_tables.py | 16 ++++-- 7 files changed, 98 insertions(+), 67 deletions(-) diff --git a/src/py21cmfast/inputs.py b/src/py21cmfast/inputs.py index d3df9ad13..b7a22c7fb 100644 --- a/src/py21cmfast/inputs.py +++ b/src/py21cmfast/inputs.py @@ -315,6 +315,12 @@ class GlobalParams(StructInstanceWrapper): AVG_BELOW_SAMPLER: int flag to add the average halo proeprties in each cell between the source mass limit and the mass sampled by the halo sampler + UPPER_STELLAR_TURNOVER_MASS: + The pivot mass associated with the optional upper mass power-law of the stellar-halo mass relation + (see FlagOptions.USE_UPPER_STELLAR_TURNOVER) + UPPER_STELLAR_TURNOVER_INDEX: + The power-law index associated with the optional upper mass power-law of the stellar-halo mass relation + (see FlagOptions.USE_UPPER_STELLAR_TURNOVER) """ def __init__(self, wrapped, ffi): @@ -757,6 +763,10 @@ class FlagOptions(StructWithDefaults): This is part of the perspective shift (see Davies & Furlanetto 2021) from counting photons/atoms in a sphere and flagging a central pixel to counting photons which we expect to reach the central pixel, and taking the ratio of atoms in the pixel. This flag simply turns off the filtering of N_rec grids, and takes the recombinations in the central cell. + USE_UPPER_STELLAR_TURNOVER: bool, optional + Whether to use an additional powerlaw in stellar mass fraction at high halo mass. The pivot mass scale and power-law index are + controlled by two global parameters, UPPER_STELLAR_TURNOVER_MASS and UPPER_STELLAR_TURNOVER_INDEX respectively. + This is currently only implemented in the halo model (USE_HALO_FIELD=True), and has no effect otherwise. """ _ffi = ffi @@ -778,6 +788,7 @@ class FlagOptions(StructWithDefaults): "FIXED_HALO_GRIDS": False, "CELL_RECOMB": False, "PHOTON_CONS_TYPE": 0, # Should these all be boolean? + "USE_UPPER_STELLAR_TURNOVER": True, } @property diff --git a/src/py21cmfast/src/21cmFAST.h b/src/py21cmfast/src/21cmFAST.h index 9b020ae7e..6f45f6703 100644 --- a/src/py21cmfast/src/21cmFAST.h +++ b/src/py21cmfast/src/21cmFAST.h @@ -92,6 +92,7 @@ struct FlagOptions{ bool FIXED_HALO_GRIDS; bool CELL_RECOMB; int PHOTON_CONS_TYPE; + bool USE_UPPER_STELLAR_TURNOVER; }; diff --git a/src/py21cmfast/src/Globals.h b/src/py21cmfast/src/Globals.h index 0508a6907..2bb12290f 100644 --- a/src/py21cmfast/src/Globals.h +++ b/src/py21cmfast/src/Globals.h @@ -91,6 +91,8 @@ struct GlobalParams{ double MIN_LOGPROB; int SAMPLE_METHOD; int AVG_BELOW_SAMPLER; + double UPPER_STELLAR_TURNOVER_MASS; + double UPPER_STELLAR_TURNOVER_INDEX; bool USE_ADIABATIC_FLUCTUATIONS; }; @@ -176,6 +178,8 @@ extern struct GlobalParams global_params = { .MIN_LOGPROB = -16, .SAMPLE_METHOD = 0, .AVG_BELOW_SAMPLER = 0, + .UPPER_STELLAR_TURNOVER_MASS = 2.8e11, + .UPPER_STELLAR_TURNOVER_INDEX = -0.61, .USE_ADIABATIC_FLUCTUATIONS = 1, }; diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index 22efd9da5..ea20d54ed 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/src/py21cmfast/src/HaloBox.c @@ -30,9 +30,10 @@ void set_halo_properties(float halo_mass, float M_turn_a, float M_turn_m, float fesc = fmin(norm_esc_var*pow(halo_mass/1e10,alpha_esc_var),1); //We don't want an upturn even with a negative ALPHA_STAR - if(astro_params_stoc->ALPHA_STAR > -0.61){ - fstar_mean = f10 * exp(-M_turn_a/halo_mass) * pow(2.6e11/1e10,astro_params_stoc->ALPHA_STAR); - fstar_mean /= pow(halo_mass/2.8e11,-astro_params_stoc->ALPHA_STAR) + pow(halo_mass/2.8e11,0.61); + if(flag_options_stoc->USE_UPPER_STELLAR_TURNOVER && (astro_params_stoc->ALPHA_STAR > global_params.UPPER_STELLAR_TURNOVER_INDEX)){ + fstar_mean = f10 * exp(-M_turn_a/halo_mass) * pow(global_params.UPPER_STELLAR_TURNOVER_MASS/1e10,astro_params_stoc->ALPHA_STAR); + fstar_mean /= pow(halo_mass/global_params.UPPER_STELLAR_TURNOVER_MASS,-astro_params_stoc->ALPHA_STAR) + + pow(halo_mass/global_params.UPPER_STELLAR_TURNOVER_MASS,-global_params.UPPER_STELLAR_TURNOVER_INDEX); } else{ fstar_mean = f10 * pow(halo_mass/1e10,fa) * exp(-M_turn_a/halo_mass); @@ -112,20 +113,29 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M double prefactor_sfr = RHOcrit * cosmo_params_stoc->OMb * norm_star / t_star / t_h; double prefactor_sfr_mini = RHOcrit * cosmo_params_stoc->OMb * norm_star_mini / t_star / t_h; - double Mlim_Fstar = Mass_limit_bisection(M_min, M_cell, alpha_star, norm_star); - double Mlim_Fesc = Mass_limit_bisection(M_min, M_cell, alpha_esc, norm_esc); + double Mlim_Fstar = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, alpha_star, norm_star); + double Mlim_Fesc = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, alpha_esc, norm_esc); - double Mlim_Fstar_mini = Mass_limit_bisection(M_min, M_cell, alpha_star_mini, norm_star_mini * pow(1e3,alpha_star_mini)); - double Mlim_Fesc_mini = Mass_limit_bisection(M_min, M_cell, alpha_esc, norm_esc_mini * pow(1e3,alpha_esc)); + double Mlim_Fstar_mini = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, alpha_star_mini, norm_star_mini * pow(1e3,alpha_star_mini)); + double Mlim_Fesc_mini = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, alpha_esc, norm_esc_mini * pow(1e3,alpha_esc)); double hm_avg=0, nion_avg=0, sfr_avg=0, sfr_avg_mini=0, wsfr_avg=0; double Mlim_a_avg=0, Mlim_m_avg=0, Mlim_r_avg=0; - //interptable variables - double min_density = -1; - double max_density = MAX_DELTAC_FRAC*Deltac; - - double delta_crit = get_delta_crit(user_params_stoc->HMF,sigma_cell,growth_z); + //find density grid limits + double min_density = 0.; + double max_density = 0.; + #pragma omp parallel num_threads(user_params_stoc->N_THREADS) + { + int i; + double dens; + #pragma omp for reduction(min:min_density) reduction(max:max_density) + for(i=0;idensity[i]; + if(dens > max_density) max_density = dens; + if(dens < min_density) min_density = dens; + } + } struct parameters_gsl_MF_integrals params = { .redshift = redshift, @@ -137,6 +147,7 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M double M_turn_a_nofb = flag_options_stoc->USE_MINI_HALOS ? atomic_cooling_threshold(redshift) : astro_params_stoc->M_TURN; LOG_DEBUG("Mean halo boxes || M = [%.2e %.2e] | Mcell = %.2e (s=%.2e) | z = %.2e | D = %.2e | cellvol = %.2e",M_min,M_max,M_cell,sigma_cell,redshift,growth_z,cell_volume); + LOG_DEBUG("Power law limits || Fstar = %.4e | Fesc = %.4e",Mlim_Fstar,Mlim_Fesc); //These tables are coarser than needed, an initial loop to find limits may help if(user_params_stoc->USE_INTERPOLATION_TABLES){ @@ -170,7 +181,7 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M #pragma omp parallel num_threads(user_params_stoc->N_THREADS) { int i; - double dens=0; + double dens; double mass=0, nion=0, sfr=0, h_count=0; double nion_mini=0, sfr_mini=0; double wsfr=0; @@ -198,41 +209,16 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M M_turn_m = M_turn_r > M_turn_m_nofb ? M_turn_r : M_turn_m_nofb; } - //ignore very low density NOTE:not using DELTA_MIN since it's perturbed (Eulerian) - if(dens <= -1){ - mass = 0.; - nion = 0.; - sfr = 0.; - h_count = 0; - } - //turn into one large halo if we exceed the critical - //Since these are perturbed (Eulerian) grids, I use the total cell mass (1+dens) - else if(dens>=MAX_DELTAC_FRAC*delta_crit){ - if(M_cell <= M_max){ - mass = M_cell * (1+dens) / cell_volume; - nion = prefactor_nion * M_cell * (1+dens) / cosmo_params_stoc->OMm / RHOcrit * pow(M_cell*(1+dens)/1e10,alpha_star) * pow(M_cell*(1+dens)/1e10,alpha_esc) / cell_volume; - sfr = prefactor_sfr * M_cell * (1+dens) / cosmo_params_stoc->OMm / RHOcrit * pow(M_cell*(1+dens)/1e10,alpha_star) / cell_volume; - h_count = 1; - } - else{ - //here the cell delta is above critical, but the integrals do not include the cell mass, so we set to zero - //NOTE: this function gives integral [M_min,M_limit] given a cell - mass = 0.; - nion = 0.; - sfr = 0.; - h_count = 0; - } - } - else{ - h_count = EvaluateNhalo(dens, growth_z, lnMmin, lnMmax, M_cell, sigma_cell, dens); - mass = EvaluateMcoll(dens, growth_z, lnMmin, lnMmax, M_cell, sigma_cell, dens); - nion = EvaluateNion_Conditional(dens,log10(M_turn_a),growth_z,M_min,M_max,M_cell,sigma_cell,Mlim_Fstar,Mlim_Fesc,false); - sfr = EvaluateSFRD_Conditional(dens,growth_z,M_min,M_max,M_cell,sigma_cell,log10(M_turn_a),Mlim_Fstar); - if(flag_options_stoc->USE_MINI_HALOS){ - sfr_mini = EvaluateSFRD_Conditional_MINI(dens,log10(M_turn_m),growth_z,M_min,M_max,M_cell,sigma_cell,log10(M_turn_a),Mlim_Fstar); - nion_mini = EvaluateNion_Conditional_MINI(dens,log10(M_turn_m),growth_z,M_min,M_max,M_cell,sigma_cell,log10(M_turn_a),Mlim_Fstar,Mlim_Fesc,false); - } + + h_count = EvaluateNhalo(dens, growth_z, lnMmin, lnMmax, M_cell, sigma_cell, dens); + mass = EvaluateMcoll(dens, growth_z, lnMmin, lnMmax, M_cell, sigma_cell, dens); + nion = EvaluateNion_Conditional(dens,log10(M_turn_a),growth_z,M_min,M_max,M_cell,sigma_cell,Mlim_Fstar,Mlim_Fesc,false); + sfr = EvaluateSFRD_Conditional(dens,growth_z,M_min,M_max,M_cell,sigma_cell,log10(M_turn_a),Mlim_Fstar); + if(flag_options_stoc->USE_MINI_HALOS){ + sfr_mini = EvaluateSFRD_Conditional_MINI(dens,log10(M_turn_m),growth_z,M_min,M_max,M_cell,sigma_cell,log10(M_turn_a),Mlim_Fstar); + nion_mini = EvaluateNion_Conditional_MINI(dens,log10(M_turn_m),growth_z,M_min,M_max,M_cell,sigma_cell,log10(M_turn_a),Mlim_Fstar,Mlim_Fesc,false); } + grids->halo_mass[i] = mass * prefactor_mass * (1+dens); grids->n_ion[i] = (nion*prefactor_nion + nion_mini*prefactor_nion_mini) * (1+dens); grids->halo_sfr[i] = (sfr*prefactor_sfr) * (1+dens); @@ -264,8 +250,8 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M dens, M_turn_m, M_turn_a, alpha_star_mini, 0., norm_star_mini, 1., Mlim_Fstar_mini, 0., user_params_stoc->INTEGRATION_METHOD_MINI)); - LOG_SUPER_DEBUG("ACG MTURN %.3e (%.3e) MCG MTURN %.3e (%.3e) REION %.3e",M_turn_a_nofb,M_turn_a,M_turn_m_nofb,M_turn_m,M_turn_r); } + LOG_SUPER_DEBUG("ACG MTURN %.3e (%.3e) MCG MTURN %.3e (%.3e) REION %.3e",M_turn_a_nofb,M_turn_a,M_turn_m_nofb,M_turn_m,M_turn_r); } hm_avg += grids->halo_mass[i]; diff --git a/src/py21cmfast/src/interp_tables.c b/src/py21cmfast/src/interp_tables.c index 25a2f22d0..b97e46299 100644 --- a/src/py21cmfast/src/interp_tables.c +++ b/src/py21cmfast/src/interp_tables.c @@ -392,7 +392,7 @@ void initialise_SFRD_Conditional_table(double min_density, double max_density, d } SFRD_conditional_table.x_min = min_density; SFRD_conditional_table.x_width = (max_density - min_density)/(NDELTA-1.); - Mlim_Fstar = Mass_limit_bisection(Mmin, Mmax, Alpha_star, Fstar10); + Mlim_Fstar = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, Alpha_star, Fstar10); if(minihalos){ if(!SFRD_conditional_table_MINI.allocated){ @@ -402,7 +402,7 @@ void initialise_SFRD_Conditional_table(double min_density, double max_density, d SFRD_conditional_table_MINI.x_width = (max_density - min_density)/(NDELTA-1.); SFRD_conditional_table_MINI.y_min = LOG10_MTURN_MIN; SFRD_conditional_table_MINI.y_width = (LOG10_MTURN_MAX - LOG10_MTURN_MIN)/(NMTURN-1.); - Mlim_Fstar_MINI = Mass_limit_bisection(Mmin, Mmax, Alpha_star_mini, Fstar7_MINI * pow(1e3, Alpha_star_mini)); + Mlim_Fstar_MINI = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, Alpha_star_mini, Fstar7_MINI * pow(1e3, Alpha_star_mini)); } double lnMmin = log(Mmin); @@ -416,7 +416,6 @@ void initialise_SFRD_Conditional_table(double min_density, double max_density, d for (i=0; iINTEGRATION_METHOD_ATOMIC)); @@ -772,7 +771,7 @@ double EvaluateNhalo(double condition, double growthf, double lnMmin, double lnM double EvaluateMcoll(double condition, double growthf, double lnMmin, double lnMmax, double M_cond, double sigma, double delta){ if(user_params_it->USE_INTERPOLATION_TABLES) return EvaluateRGTable1D(condition,&Mcoll_table); - return Nhalo_Conditional(growthf, lnMmax, lnMmax, M_cond, sigma, delta, user_params_it->INTEGRATION_METHOD_HALOS); + return Mcoll_Conditional(growthf, lnMmax, lnMmax, M_cond, sigma, delta, user_params_it->INTEGRATION_METHOD_HALOS); } //This one is always a table diff --git a/src/py21cmfast/src/ps.c b/src/py21cmfast/src/ps.c index d24aa762c..f3e9e7c5b 100644 --- a/src/py21cmfast/src/ps.c +++ b/src/py21cmfast/src/ps.c @@ -1550,9 +1550,15 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond .delta = delta, }; - //return 1 halo if delta is exceeded - if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)) - return 1./M_cond; + //return 1 halo AT THE CONDITION MASS if delta is exceeded + if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ + if(M_cond <= exp(lnM2)) + return 1./M_cond; + else + return 0.; + } + if(delta <= -1.) + return 0. return IntegratedNdM(lnM1,lnM2,params,-1, method); } @@ -1565,9 +1571,15 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond .delta = delta, }; - //return 100% of mass if delta is exceeded - if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)) - return 1.; + //return 100% of mass AT THE CONDITION MASS if delta is exceeded + if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ + if(M_cond <= exp(lnM2)) + return 1.; + else + return 0.; + } + if(delta <= -1.) + return 0. return IntegratedNdM(lnM1,lnM2,params,-2, method); } @@ -1594,8 +1606,14 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M //return 1 halo at the condition mass if delta is exceeded //NOTE: this will almost always be zero, due to the upper turover, // however this replaces an integral so it won't be slow - if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)) - return nion_fraction_mini(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) + if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ + if(M_cond <= exp(lnM2)) + return nion_fraction_mini(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) + else + return 0.; + } + if(delta <= -1.) + return 0. // LOG_ULTRA_DEBUG("params: D=%.2e Mtl=%.2e Mtu=%.2e as=%.2e ae=%.2e fs=%.2e fe=%.2e Ms=%.2e Me=%.2e hmf=%d sig=%.2e del=%.2e", // growthf,MassTurnover,MassTurnover_upper,Alpha_star,Alpha_esc,Fstar7,Fesc7,Mlim_Fstar,Mlim_Fesc,0,sigma2,delta2); @@ -1621,8 +1639,14 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond }; //return 1 halo at the condition mass if delta is exceeded - if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)) - return nion_fraction(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) + if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ + if(M_cond <= exp(lnM2)) + return nion_fraction(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) + else + return 0. + } + if(delta <= -1.) + return 0. // LOG_ULTRA_DEBUG("params: D=%.2e Mtl=%.2e as=%.2e ae=%.2e fs=%.2e fe=%.2e Ms=%.2e Me=%.2e sig=%.2e del=%.2e", // growthf,MassTurnover,Alpha_star,Alpha_esc,Fstar10,Fesc10,Mlim_Fstar,Mlim_Fesc,sigma2,delta2); diff --git a/tests/test_c_interpolation_tables.py b/tests/test_c_interpolation_tables.py index d86846861..0a0e4b08a 100644 --- a/tests/test_c_interpolation_tables.py +++ b/tests/test_c_interpolation_tables.py @@ -87,6 +87,7 @@ def test_sigma_table(name, plt): make_table_comparison_plot( mass_range, np.array([0]), + np.array([0]), sigma_table[:, None], dsigmasq_table[:, None], sigma_ref[:, None], @@ -389,6 +390,7 @@ def test_FgtrM_conditional_tables(name, R, plt): make_table_comparison_plot( edges_d[:-1], np.array([0]), + np.array([0]), fcoll_tables[:, None], dfcoll_tables[:, None], fcoll_integrals[:, None], @@ -524,6 +526,7 @@ def test_SFRD_z_tables(name, plt): sel_m = (xl * np.arange(6) / 6).astype(int) make_table_comparison_plot( z_array[:-1], + np.array([0]), edges_m[sel_m], SFRD_tables[:, None], SFRD_tables_mini[..., sel_m], @@ -665,6 +668,7 @@ def test_Nion_z_tables(name, plt): sel_m = (xl * np.arange(6) / 6).astype(int) make_table_comparison_plot( z_array[:-1], + np.array([0]), edges_m[sel_m], Nion_tables[:, None], Nion_tables_mini[..., sel_m], @@ -885,6 +889,7 @@ def test_Nion_conditional_tables(name, R, mini, intmethod, plt): make_table_comparison_plot( edges_d[:-1], edges_m[sel_m], + edges_m[sel_m], Nion_tb_plot, Nion_tables_mini[..., sel_m], Nion_il_plot, @@ -1047,6 +1052,7 @@ def test_SFRD_conditional_table(name, R, intmethod, plt): sel_m = (xl * np.arange(6) / 6).astype(int) make_table_comparison_plot( edges_d[:-1], + np.array([0]), edges_m[sel_m], SFRD_tables[:, None], SFRD_tables_mini[..., sel_m], @@ -1218,10 +1224,10 @@ def test_conditional_integral_methods(R, name, integrand, plt): ) -def make_table_comparison_plot(x1, x2, table_1d, table_2d, intgrl_1d, intgrl_2d, plt): +def make_table_comparison_plot(x1, tb1_z, tb2_z, table_1d, table_2d, intgrl_1d, intgrl_2d, plt): # rows = values,fracitonal diff, cols = 1d table, 2d table fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(16, 16)) - for i in range(x1.size): + for i in range(tb1_z.size): make_comparison_plot( x1, intgrl_1d[:, i], @@ -1229,12 +1235,12 @@ def make_table_comparison_plot(x1, x2, table_1d, table_2d, intgrl_1d, intgrl_2d, ax=axs[:, 0], xlab="delta", ylab="MF integral", - label_base=f"Mt = {x2[i]:.1e}" if x1.size > 1 else "", + label_base=f"Mt = {tb1_z[i]:.1e}" if tb1_z.size > 1 else "", logx=False, color=f"C{i:d}", ) - for i in range(x2.size): + for i in range(tb2_z.size): make_comparison_plot( x1, intgrl_2d[:, i], @@ -1242,7 +1248,7 @@ def make_table_comparison_plot(x1, x2, table_1d, table_2d, intgrl_1d, intgrl_2d, ax=axs[:, 1], xlab="delta", ylab="MF integral mini", - label_base=f"Mt = {x2[i]:.1e}" if x2.size > 1 else "", + label_base=f"Mt = {tb2_z[i]:.1e}" if tb2_z.size > 1 else "", logx=False, color=f"C{i:d}", ) From 01d748ba5b36fcb71f6575488900f7cd92a3293b Mon Sep 17 00:00:00 2001 From: James Davies Date: Mon, 8 Apr 2024 11:30:44 +0200 Subject: [PATCH 2/5] further bugfixes in fixed grids and crit delta limits --- src/py21cmfast/src/HaloBox.c | 8 ++------ src/py21cmfast/src/ps.c | 22 +++++++++++----------- tests/test_c_interpolation_tables.py | 4 +++- 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index ea20d54ed..36cbd27d6 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/src/py21cmfast/src/HaloBox.c @@ -149,11 +149,9 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M LOG_DEBUG("Mean halo boxes || M = [%.2e %.2e] | Mcell = %.2e (s=%.2e) | z = %.2e | D = %.2e | cellvol = %.2e",M_min,M_max,M_cell,sigma_cell,redshift,growth_z,cell_volume); LOG_DEBUG("Power law limits || Fstar = %.4e | Fesc = %.4e",Mlim_Fstar,Mlim_Fesc); - //These tables are coarser than needed, an initial loop to find limits may help + //These tables are coarser than needed, an initial loop for Mturn to find limits may help if(user_params_stoc->USE_INTERPOLATION_TABLES){ if(user_params_stoc->INTEGRATION_METHOD_ATOMIC == 1 || user_params_stoc->INTEGRATION_METHOD_MINI == 1){ - M_max = M_cell; //we need these to be the same for the "smoothness" assumption of the GL integration - lnMmax = lnMcell; initialise_GL(NGL_INT, lnMmin, lnMmax); } @@ -209,7 +207,6 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M M_turn_m = M_turn_r > M_turn_m_nofb ? M_turn_r : M_turn_m_nofb; } - h_count = EvaluateNhalo(dens, growth_z, lnMmin, lnMmax, M_cell, sigma_cell, dens); mass = EvaluateMcoll(dens, growth_z, lnMmin, lnMmax, M_cell, sigma_cell, dens); nion = EvaluateNion_Conditional(dens,log10(M_turn_a),growth_z,M_min,M_max,M_cell,sigma_cell,Mlim_Fstar,Mlim_Fesc,false); @@ -333,7 +330,7 @@ int get_box_averages(double redshift, double norm_esc, double alpha_esc, double initialise_GL(NGL_INT, lnMmin, lnMmax); //NOTE: we use the atomic method for all halo mass/count here - hm_expected = IntegratedNdM(lnMmin,lnMmax,params,2,user_params_stoc->INTEGRATION_METHOD_ATOMIC); + hm_expected = FgtrM_General(redshift,M_min) * prefactor_mass; nion_expected = Nion_General(redshift, lnMmin, lnMmax, M_turn_a, alpha_star, alpha_esc, norm_star, norm_esc, Mlim_Fstar, Mlim_Fesc) * prefactor_nion; sfr_expected = Nion_General(redshift, lnMmin, lnMmax, M_turn_a, alpha_star, 0., norm_star, 1., @@ -348,7 +345,6 @@ int get_box_averages(double redshift, double norm_esc, double alpha_esc, double 1., Mlim_Fstar_mini, 0.) * prefactor_sfr_mini; } - // hm_expected *= prefactor_mass; //for non-CMF, the factors are already there wsfr_expected = nion_expected / t_star / t_h; //same integral, different prefactors, different in the stochastic grids due to scatter averages[0] = hm_expected; diff --git a/src/py21cmfast/src/ps.c b/src/py21cmfast/src/ps.c index f3e9e7c5b..d5c7ea6ef 100644 --- a/src/py21cmfast/src/ps.c +++ b/src/py21cmfast/src/ps.c @@ -1552,13 +1552,13 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond //return 1 halo AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ - if(M_cond <= exp(lnM2)) + if(M_cond*0.99 <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) return 1./M_cond; else return 0.; } if(delta <= -1.) - return 0. + return 0.; return IntegratedNdM(lnM1,lnM2,params,-1, method); } @@ -1573,13 +1573,13 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond //return 100% of mass AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ - if(M_cond <= exp(lnM2)) + if(M_cond*0.99 <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) return 1.; else return 0.; } if(delta <= -1.) - return 0. + return 0.; return IntegratedNdM(lnM1,lnM2,params,-2, method); } @@ -1607,13 +1607,13 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M //NOTE: this will almost always be zero, due to the upper turover, // however this replaces an integral so it won't be slow if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ - if(M_cond <= exp(lnM2)) + if(M_cond*0.99 <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) return nion_fraction_mini(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) else return 0.; } - if(delta <= -1.) - return 0. + if(delta2 <= -1.) + return 0.; // LOG_ULTRA_DEBUG("params: D=%.2e Mtl=%.2e Mtu=%.2e as=%.2e ae=%.2e fs=%.2e fe=%.2e Ms=%.2e Me=%.2e hmf=%d sig=%.2e del=%.2e", // growthf,MassTurnover,MassTurnover_upper,Alpha_star,Alpha_esc,Fstar7,Fesc7,Mlim_Fstar,Mlim_Fesc,0,sigma2,delta2); @@ -1640,13 +1640,13 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond //return 1 halo at the condition mass if delta is exceeded if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ - if(M_cond <= exp(lnM2)) + if(M_cond*0.99 <= exp(lnM2)) return nion_fraction(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) else - return 0. + return 0.; } - if(delta <= -1.) - return 0. + if(delta2 <= -1.) + return 0.; // LOG_ULTRA_DEBUG("params: D=%.2e Mtl=%.2e as=%.2e ae=%.2e fs=%.2e fe=%.2e Ms=%.2e Me=%.2e sig=%.2e del=%.2e", // growthf,MassTurnover,Alpha_star,Alpha_esc,Fstar10,Fesc10,Mlim_Fstar,Mlim_Fesc,sigma2,delta2); diff --git a/tests/test_c_interpolation_tables.py b/tests/test_c_interpolation_tables.py index 0a0e4b08a..d66a221e6 100644 --- a/tests/test_c_interpolation_tables.py +++ b/tests/test_c_interpolation_tables.py @@ -169,8 +169,10 @@ def test_Massfunc_conditional_tables(name): up.INTEGRATION_METHOD_HALOS, ) + max_in_d = N_cmfi_cell[:, :1] + max_in_d[:,0] = 1. #fix delta=-1 N_cmfi_cell = ( - N_cmfi_cell / N_cmfi_cell[:, :1] + N_cmfi_cell / max_in_d ) # to get P(>M) since the y-axis is the lower integral limit mask_cell = N_cmfi_cell > np.exp( global_params.MIN_LOGPROB From 86ad74d1b73b90b7fa475a5157697a83e6701227 Mon Sep 17 00:00:00 2001 From: James Davies Date: Mon, 8 Apr 2024 12:00:01 +0200 Subject: [PATCH 3/5] set critical density mass-limit tolerance in integrals using FRACT_FLOAT_ERROR --- src/py21cmfast/src/ps.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/py21cmfast/src/ps.c b/src/py21cmfast/src/ps.c index d5c7ea6ef..92d3049ba 100644 --- a/src/py21cmfast/src/ps.c +++ b/src/py21cmfast/src/ps.c @@ -1552,7 +1552,7 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond //return 1 halo AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ - if(M_cond*0.99 <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) + if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) return 1./M_cond; else return 0.; @@ -1573,7 +1573,7 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond //return 100% of mass AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ - if(M_cond*0.99 <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) + if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) return 1.; else return 0.; @@ -1607,7 +1607,7 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M //NOTE: this will almost always be zero, due to the upper turover, // however this replaces an integral so it won't be slow if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ - if(M_cond*0.99 <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) + if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) return nion_fraction_mini(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) else return 0.; @@ -1640,7 +1640,7 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond //return 1 halo at the condition mass if delta is exceeded if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ - if(M_cond*0.99 <= exp(lnM2)) + if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) return nion_fraction(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) else return 0.; From 28b9a294cd9661b6c3dcd98e928b79fe96685e32 Mon Sep 17 00:00:00 2001 From: James Davies Date: Mon, 8 Apr 2024 12:44:51 +0200 Subject: [PATCH 4/5] fix massfunc table test maxima --- tests/test_c_interpolation_tables.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/test_c_interpolation_tables.py b/tests/test_c_interpolation_tables.py index d66a221e6..7d56e1dd0 100644 --- a/tests/test_c_interpolation_tables.py +++ b/tests/test_c_interpolation_tables.py @@ -170,7 +170,7 @@ def test_Massfunc_conditional_tables(name): ) max_in_d = N_cmfi_cell[:, :1] - max_in_d[:,0] = 1. #fix delta=-1 + max_in_d[edges_d[:-1] == -1, :] = 1.0 # fix delta=-1 where all entries are zero N_cmfi_cell = ( N_cmfi_cell / max_in_d ) # to get P(>M) since the y-axis is the lower integral limit @@ -1226,7 +1226,9 @@ def test_conditional_integral_methods(R, name, integrand, plt): ) -def make_table_comparison_plot(x1, tb1_z, tb2_z, table_1d, table_2d, intgrl_1d, intgrl_2d, plt): +def make_table_comparison_plot( + x1, tb1_z, tb2_z, table_1d, table_2d, intgrl_1d, intgrl_2d, plt +): # rows = values,fracitonal diff, cols = 1d table, 2d table fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(16, 16)) for i in range(tb1_z.size): From 44ad87b7b575f9fd264315d67caa72792658ed03 Mon Sep 17 00:00:00 2001 From: James Davies Date: Mon, 8 Apr 2024 13:48:21 +0200 Subject: [PATCH 5/5] move upper turnover parameters to AstroParams --- src/py21cmfast/inputs.py | 15 +++++++++------ src/py21cmfast/src/21cmFAST.h | 3 +++ src/py21cmfast/src/Globals.h | 4 ---- src/py21cmfast/src/HaloBox.c | 8 ++++---- 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/py21cmfast/inputs.py b/src/py21cmfast/inputs.py index b7a22c7fb..f86e75e9f 100644 --- a/src/py21cmfast/inputs.py +++ b/src/py21cmfast/inputs.py @@ -315,12 +315,6 @@ class GlobalParams(StructInstanceWrapper): AVG_BELOW_SAMPLER: int flag to add the average halo proeprties in each cell between the source mass limit and the mass sampled by the halo sampler - UPPER_STELLAR_TURNOVER_MASS: - The pivot mass associated with the optional upper mass power-law of the stellar-halo mass relation - (see FlagOptions.USE_UPPER_STELLAR_TURNOVER) - UPPER_STELLAR_TURNOVER_INDEX: - The power-law index associated with the optional upper mass power-law of the stellar-halo mass relation - (see FlagOptions.USE_UPPER_STELLAR_TURNOVER) """ def __init__(self, wrapped, ffi): @@ -991,6 +985,12 @@ class AstroParams(StructWithDefaults): Impact of the LW feedback on Mturn for minihaloes. Default is 22.8685 and 0.47 following Machacek+01, respectively. Latest simulations suggest 2.0 and 0.6. See Sec 2 of Muñoz+21 (2110.13919). A_VCB, BETA_VCB: float, optional Impact of the DM-baryon relative velocities on Mturn for minihaloes. Default is 1.0 and 1.8, and agrees between different sims. See Sec 2 of Muñoz+21 (2110.13919). + UPPER_STELLAR_TURNOVER_MASS: + The pivot mass associated with the optional upper mass power-law of the stellar-halo mass relation + (see FlagOptions.USE_UPPER_STELLAR_TURNOVER) + UPPER_STELLAR_TURNOVER_INDEX: + The power-law index associated with the optional upper mass power-law of the stellar-halo mass relation + (see FlagOptions.USE_UPPER_STELLAR_TURNOVER) """ _ffi = ffi @@ -1023,6 +1023,8 @@ class AstroParams(StructWithDefaults): "BETA_LW": 0.6, "A_VCB": 1.0, "BETA_VCB": 1.8, + "UPPER_STELLAR_TURNOVER_MASS": 11.447, # 2.8e11 + "UPPER_STELLAR_TURNOVER_INDEX": -0.61, } def __init__( @@ -1045,6 +1047,7 @@ def convert(self, key, val): "L_X", "L_X_MINI", "X_RAY_Tvir_MIN", + "UPPER_STELLAR_TURNOVER_MASS", ]: return 10**val else: diff --git a/src/py21cmfast/src/21cmFAST.h b/src/py21cmfast/src/21cmFAST.h index 6f45f6703..f771a778e 100644 --- a/src/py21cmfast/src/21cmFAST.h +++ b/src/py21cmfast/src/21cmFAST.h @@ -71,6 +71,9 @@ struct AstroParams{ float t_STAR; int N_RSD_STEPS; + + double UPPER_STELLAR_TURNOVER_MASS; + double UPPER_STELLAR_TURNOVER_INDEX; }; struct FlagOptions{ diff --git a/src/py21cmfast/src/Globals.h b/src/py21cmfast/src/Globals.h index 2bb12290f..0508a6907 100644 --- a/src/py21cmfast/src/Globals.h +++ b/src/py21cmfast/src/Globals.h @@ -91,8 +91,6 @@ struct GlobalParams{ double MIN_LOGPROB; int SAMPLE_METHOD; int AVG_BELOW_SAMPLER; - double UPPER_STELLAR_TURNOVER_MASS; - double UPPER_STELLAR_TURNOVER_INDEX; bool USE_ADIABATIC_FLUCTUATIONS; }; @@ -178,8 +176,6 @@ extern struct GlobalParams global_params = { .MIN_LOGPROB = -16, .SAMPLE_METHOD = 0, .AVG_BELOW_SAMPLER = 0, - .UPPER_STELLAR_TURNOVER_MASS = 2.8e11, - .UPPER_STELLAR_TURNOVER_INDEX = -0.61, .USE_ADIABATIC_FLUCTUATIONS = 1, }; diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index 36cbd27d6..fe333efeb 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/src/py21cmfast/src/HaloBox.c @@ -30,10 +30,10 @@ void set_halo_properties(float halo_mass, float M_turn_a, float M_turn_m, float fesc = fmin(norm_esc_var*pow(halo_mass/1e10,alpha_esc_var),1); //We don't want an upturn even with a negative ALPHA_STAR - if(flag_options_stoc->USE_UPPER_STELLAR_TURNOVER && (astro_params_stoc->ALPHA_STAR > global_params.UPPER_STELLAR_TURNOVER_INDEX)){ - fstar_mean = f10 * exp(-M_turn_a/halo_mass) * pow(global_params.UPPER_STELLAR_TURNOVER_MASS/1e10,astro_params_stoc->ALPHA_STAR); - fstar_mean /= pow(halo_mass/global_params.UPPER_STELLAR_TURNOVER_MASS,-astro_params_stoc->ALPHA_STAR) - + pow(halo_mass/global_params.UPPER_STELLAR_TURNOVER_MASS,-global_params.UPPER_STELLAR_TURNOVER_INDEX); + if(flag_options_stoc->USE_UPPER_STELLAR_TURNOVER && (astro_params_stoc->ALPHA_STAR > astro_params_stoc->UPPER_STELLAR_TURNOVER_INDEX)){ + fstar_mean = f10 * exp(-M_turn_a/halo_mass) * pow(astro_params_stoc->UPPER_STELLAR_TURNOVER_MASS/1e10,astro_params_stoc->ALPHA_STAR); + fstar_mean /= pow(halo_mass/astro_params_stoc->UPPER_STELLAR_TURNOVER_MASS,-astro_params_stoc->ALPHA_STAR) + + pow(halo_mass/astro_params_stoc->UPPER_STELLAR_TURNOVER_MASS,-astro_params_stoc->UPPER_STELLAR_TURNOVER_INDEX); } else{ fstar_mean = f10 * pow(halo_mass/1e10,fa) * exp(-M_turn_a/halo_mass);