From eb4df53e062f422158e000c12e4fdc4ec42e64bc Mon Sep 17 00:00:00 2001 From: James Davies Date: Thu, 14 Nov 2024 15:35:22 +0100 Subject: [PATCH 1/4] fix trapezoidal bug, make error messages more clear --- src/py21cmfast/src/IonisationBox.c | 43 +++++++++++++++++++----------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/src/py21cmfast/src/IonisationBox.c b/src/py21cmfast/src/IonisationBox.c index 066bb6364..5bd68ac39 100644 --- a/src/py21cmfast/src/IonisationBox.c +++ b/src/py21cmfast/src/IonisationBox.c @@ -138,6 +138,7 @@ struct FilteredGrids{ void set_ionbox_constants(double redshift, double prev_redshift, CosmoParams *cosmo_params, UserParams *user_params, AstroParams *astro_params, FlagOptions *flag_options, struct IonBoxConstants *consts){ consts->redshift = redshift; + consts->prev_redshift = prev_redshift; //defaults for no photoncons consts->stored_redshift = redshift; consts->photoncons_adjustment_factor = 1.; @@ -485,15 +486,19 @@ void set_mean_fcoll(struct IonBoxConstants *c, IonizedBox *prev_box, IonizedBox *f_limit_acg = Fcoll_General(user_params_global->Z_HEAT_MAX, c->lnMmin, c->lnMmax_gl); //JD: the old parametrisation didn't have this limit before } - if(isfinite(curr_box->mean_f_coll)==0){ - LOG_ERROR("Mean collapse fraction is either infinite or NaN!"); + if(isfinite(curr_box->mean_f_coll)==0 || curr_box->mean_f_coll < 0){ + LOG_ERROR("Mean collapse fraction is invalid"); + LOG_ERROR("prev box %g prev intgrl %g curr intrgl %g --> %g", + prev_box->mean_f_coll, f_coll_prev, f_coll_curr, curr_box->mean_f_coll); Throw(InfinityorNaNError); } LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll: %e", curr_box->mean_f_coll); if (flag_options_global->USE_MINI_HALOS){ - if(isfinite(curr_box->mean_f_coll_MINI)==0){ - LOG_ERROR("Mean collapse fraction of MINI is either infinite or NaN!"); + if(isfinite(curr_box->mean_f_coll_MINI)==0 || curr_box->mean_f_coll_MINI < 0){ + LOG_ERROR("Mean collapse fraction is invalid"); + LOG_ERROR("prev box %g prev intgrl %g curr intrgl %g --> %g", + prev_box->mean_f_coll_MINI, f_coll_prev_mini, f_coll_curr_mini, curr_box->mean_f_coll_MINI); Throw(InfinityorNaNError); } LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll_MINI: %e", curr_box->mean_f_coll_MINI); @@ -672,6 +677,7 @@ void setup_integration_tables(struct FilteredGrids *fg_struct, struct IonBoxCons //previous redshift tables if needed if(need_prev && flag_options_global->USE_MINI_HALOS){ + //NOTE: we intentionally use the lower turnovers at this redshift, but should we be doing the same for the upper turnover? initialise_Nion_Conditional_spline(consts->prev_redshift,consts->mturn_a_nofb,prev_min_density,prev_max_density,consts->M_min,rspec.M_max_R,rspec.M_max_R, log10Mturn_min,log10Mturn_max,log10Mturn_min_MINI,log10Mturn_max_MINI, consts->alpha_star, consts->alpha_star_mini, @@ -801,12 +807,15 @@ void calculate_fcoll_grid(IonizedBox *box, IonizedBox *previous_ionize_box, stru //if (box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] < previous_ionize_box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]) // box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = previous_ionize_box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; f_coll_total += box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; - if(isfinite(f_coll_total)==0) { - LOG_ERROR("f_coll is either infinite or NaN!(%d,%d,%d)%g,%g,%g,%g,%g,%g,%g,%g,%g",\ - x,y,z,curr_dens,prev_dens,previous_ionize_box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)],\ - Splined_Fcoll, prev_Splined_Fcoll, curr_dens, prev_dens, log10_Mturnover, + if(isfinite(f_coll_total)==0){ + LOG_ERROR("f_coll is either infinite or NaN! %d %g (%d,%d,%d)?: dens %g, prev %g", + rspec->R_index,rspec->R,x,y,z,curr_dens,prev_dens); + LOG_ERROR("prev box %g intgrl %g curr %g --> %g l10mturn %g (%g)", + previous_ionize_box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)], + prev_Splined_Fcoll, Splined_Fcoll, + box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)], + log10_Mturnover, *((float *)fg_struct->log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z))); - Throw(InfinityorNaNError); } if (Splined_Fcoll_MINI > 1.) Splined_Fcoll_MINI = 1.; @@ -821,13 +830,15 @@ void calculate_fcoll_grid(IonizedBox *box, IonizedBox *previous_ionize_box, stru //if (box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] < previous_ionize_box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]) // box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = previous_ionize_box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; f_coll_MINI_total += box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; - if(isfinite(f_coll_MINI_total)==0) { - LOG_ERROR("f_coll_MINI is either infinite or NaN!(%d,%d,%d)%g,%g,%g,%g,%g,%g,%g",\ - x,y,z,curr_dens, prev_dens, previous_ionize_box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)],\ - Splined_Fcoll_MINI, prev_Splined_Fcoll_MINI, log10_Mturnover_MINI,\ - *((float *)fg_struct->log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z))); - LOG_DEBUG("%g,%g",previous_ionize_box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)],\ - previous_ionize_box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]); + if(isfinite(f_coll_MINI_total)==0){ + LOG_ERROR("f_coll_MINI is either infinite or NaN! %d %g (%d,%d,%d)?: dens %g, prev %g", + rspec->R_index,rspec->R,x,y,z,curr_dens,prev_dens); + LOG_ERROR("prev box %g intgrl %g curr int %g --> %g l10mturn %g (%g)", + previous_ionize_box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)], + prev_Splined_Fcoll_MINI, Splined_Fcoll_MINI, + box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)], + log10_Mturnover_MINI, + *((float *)fg_struct->log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z))); Throw(InfinityorNaNError); } } From e36e170af0988341f06487d153d76df56f036e9d Mon Sep 17 00:00:00 2001 From: James Davies Date: Thu, 21 Nov 2024 17:53:39 +0100 Subject: [PATCH 2/4] remove minihalo fesc dynamic default --- src/py21cmfast/wrapper/inputs.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/py21cmfast/wrapper/inputs.py b/src/py21cmfast/wrapper/inputs.py index b3215ceb8..79bff0ce1 100644 --- a/src/py21cmfast/wrapper/inputs.py +++ b/src/py21cmfast/wrapper/inputs.py @@ -712,6 +712,7 @@ class AstroParams(InputStruct): converter=float, ) F_ESC7_MINI = field( + default=-2.0, converter=float, transformer=logtransformer, ) @@ -775,14 +776,6 @@ def _F_STAR7_MINI_default(self): """ return self.F_STAR10 - 3 * self.ALPHA_STAR # -3*alpha since 1e7/1e10 = 1e-3 - # NOTE: Currently the default is not `None`, so this would normally do nothing. - # We need to examine the MCG/ACG connection to popII/popIII stars and - # discuss what this model should be. - @F_ESC7_MINI.default - def _F_ESC7_MINI_default(self): - """The stellar-to-halo mass ratio at 1e7 Solar Masses for Molecularly cooled galaxies.""" - return self.F_ESC10 - 3 * self.ALPHA_ESC # -3*alpha since 1e7/1e10 = 1e-3 - @ALPHA_STAR_MINI.default def _ALPHA_STAR_MINI_default(self): """ From c9e68d26afc133635ddccf1e1b3fdbe3f8118ffc Mon Sep 17 00:00:00 2001 From: James Davies Date: Fri, 15 Nov 2024 17:25:20 +0100 Subject: [PATCH 3/4] normalise mturn ffts properly --- src/py21cmfast/src/IonisationBox.c | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/py21cmfast/src/IonisationBox.c b/src/py21cmfast/src/IonisationBox.c index 5bd68ac39..942f72d96 100644 --- a/src/py21cmfast/src/IonisationBox.c +++ b/src/py21cmfast/src/IonisationBox.c @@ -651,7 +651,7 @@ void setup_integration_tables(struct FilteredGrids *fg_struct, struct IonBoxCons clip_and_get_extrema(fg_struct->log10_Mturnover_MINI_filtered,0.,LOG10_MTURN_MAX,&log10Mturn_min_MINI,&log10Mturn_max_MINI); } - LOG_ULTRA_DEBUG("Tb limits d (%.2e,%.2e), m (%.2e,%.2e) t (%.2e,%.2e) tm (%.2e,%.2e)", + LOG_SUPER_DEBUG("Tb limits d (%.2e,%.2e), m (%.2e,%.2e) t (%.2e,%.2e) tm (%.2e,%.2e)", min_density,max_density,consts->M_min,rspec.M_max_R,log10Mturn_min,log10Mturn_max, log10Mturn_min_MINI,log10Mturn_max_MINI); if(user_params_global->INTEGRATION_METHOD_ATOMIC == 1 || (flag_options_global->USE_MINI_HALOS && user_params_global->INTEGRATION_METHOD_MINI == 1)) @@ -831,7 +831,7 @@ void calculate_fcoll_grid(IonizedBox *box, IonizedBox *previous_ionize_box, stru // box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = previous_ionize_box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; f_coll_MINI_total += box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; if(isfinite(f_coll_MINI_total)==0){ - LOG_ERROR("f_coll_MINI is either infinite or NaN! %d %g (%d,%d,%d)?: dens %g, prev %g", + LOG_ERROR("f_coll_MINI is either infinite or NaN %d R=%g (%d,%d,%d)?: dens %g, prev %g", rspec->R_index,rspec->R,x,y,z,curr_dens,prev_dens); LOG_ERROR("prev box %g intgrl %g curr int %g --> %g l10mturn %g (%g)", previous_ionize_box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)], @@ -1312,6 +1312,10 @@ int ComputeIonizedBox(float redshift, float prev_redshift, UserParams *user_para user_params->N_THREADS, grid_struct->log10_Mturnover_MINI_unfiltered); dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, grid_struct->log10_Mturnover_unfiltered); + for (ct=0; ctlog10_Mturnover_unfiltered[ct] /= (float)HII_TOT_NUM_PIXELS; + grid_struct->log10_Mturnover_MINI_unfiltered[ct] /= (float)HII_TOT_NUM_PIXELS; + } } } if(flag_options->USE_TS_FLUCT){ From b647a89b404b2d935ad3b060f16a698067b89dbf Mon Sep 17 00:00:00 2001 From: James Davies Date: Wed, 4 Dec 2024 14:24:53 +0100 Subject: [PATCH 4/4] alter Nion_z table tolerance --- tests/test_c_interpolation_tables.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_c_interpolation_tables.py b/tests/test_c_interpolation_tables.py index 88efc05f4..443dea219 100644 --- a/tests/test_c_interpolation_tables.py +++ b/tests/test_c_interpolation_tables.py @@ -567,7 +567,7 @@ def test_Nion_z_tables(name, z_range, log10_mturn_range, plt): ylabels=["Nion", "Nion_mini"], ) - abs_tol = 1e-6 + abs_tol = 2e-6 print_failure_stats( Nion_tables, Nion_integrals,