diff --git a/src/py21cmfast/drivers/param_config.py b/src/py21cmfast/drivers/param_config.py index 5475de82..939e7de2 100644 --- a/src/py21cmfast/drivers/param_config.py +++ b/src/py21cmfast/drivers/param_config.py @@ -112,6 +112,18 @@ def _astro_params_validator(self, att, val): else: raise ValueError(msg) + @user_params.validator + def _user_params_validator(self, att, val): + if val is None: + return + # perform a very rudimentary check to see if we are underresolved and not using the linear approx + if val.BOX_LEN > val.DIM and not global_params.EVOLVE_DENSITY_LINEARLY: + warnings.warn( + "Resolution is likely too low for accurate evolved density fields\n It Is recommended" + + "that you either increase the resolution (DIM/BOX_LEN) or" + + "set the EVOLVE_DENSITY_LINEARLY flag to 1" + ) + def merge_keys(self): """The list of available structs in this instance.""" # Allow using **input_parameters diff --git a/src/py21cmfast/drivers/single_field.py b/src/py21cmfast/drivers/single_field.py index 16cf80c1..818bfeab 100644 --- a/src/py21cmfast/drivers/single_field.py +++ b/src/py21cmfast/drivers/single_field.py @@ -715,8 +715,11 @@ def compute_xray_source_field( R_outer = R_range[i].to("Mpc").value if zpp_avg[i] >= z_max: + box.filtered_sfr[i] = 0 + box.filtered_sfr_mini[i] = 0 + box.filtered_xray[i] = 0 + box.mean_log10_Mcrit_LW[i] = inputs.astro_params.M_TURN # minimum logger.debug(f"ignoring Radius {i} which is above Z_HEAT_MAX") - box.filtered_sfr[i, ...] = 0 continue hbox_interp = interp_halo_boxes( @@ -730,6 +733,7 @@ def compute_xray_source_field( box.filtered_sfr[i] = 0 box.filtered_sfr_mini[i] = 0 box.filtered_xray[i] = 0 + box.mean_log10_Mcrit_LW[i] = hbox_interp.log10_Mcrit_MCG_ave logger.debug(f"ignoring Radius {i} due to no stars") continue @@ -963,7 +967,7 @@ def compute_ionization_field( if not inputs.flag_options.USE_MINI_HALOS: previous_perturbed_field = PerturbedField( - inputs=inputs.clone(redshift=0.0), dummy=True + inputs=inputs.clone(redshift=0.0), initial=True ) elif previous_perturbed_field is None: # If we are beyond Z_HEAT_MAX, just make an empty box diff --git a/src/py21cmfast/src/BrightnessTemperatureBox.c b/src/py21cmfast/src/BrightnessTemperatureBox.c index 8c8fb27a..d9356c06 100644 --- a/src/py21cmfast/src/BrightnessTemperatureBox.c +++ b/src/py21cmfast/src/BrightnessTemperatureBox.c @@ -107,7 +107,7 @@ int ComputeBrightnessTemp(float redshift, UserParams *user_params, CosmoParams * sqrt( (0.15/(cosmo_params->OMm)/(cosmo_params->hlittle)/(cosmo_params->hlittle)) * (1.+redshift)/10.0 ); /////////////////////////////// END INITIALIZATION ///////////////////////////////////////////// - LOG_DEBUG("Performed Initialization."); + LOG_SUPER_DEBUG("Performed Initialization."); // ok, lets fill the delta_T box; which will be the same size as the bubble box #pragma omp parallel shared(const_factor,perturb_field,ionized_box,box,redshift,spin_temp,T_rad) \ @@ -141,7 +141,7 @@ int ComputeBrightnessTemp(float redshift, UserParams *user_params, CosmoParams * } } - LOG_DEBUG("Filled delta_T."); + LOG_SUPER_DEBUG("Filled delta_T."); if(isfinite(ave)==0) { LOG_ERROR("Average brightness temperature is infinite or NaN!"); @@ -211,7 +211,7 @@ int ComputeBrightnessTemp(float redshift, UserParams *user_params, CosmoParams * } } - LOG_DEBUG("Included velocities."); + LOG_SUPER_DEBUG("Included velocities."); if(isfinite(ave)==0) { @@ -233,7 +233,7 @@ int ComputeBrightnessTemp(float redshift, UserParams *user_params, CosmoParams * free(delta_T_RSD_LOS); fftwf_cleanup_threads(); fftwf_cleanup(); - LOG_DEBUG("Cleaned up."); + LOG_SUPER_DEBUG("Cleaned up."); } // End of try Catch(status){ diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index 0c920e88..fa4e8033 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/src/py21cmfast/src/HaloBox.c @@ -28,6 +28,7 @@ struct HaloBoxConstants{ double redshift; bool fix_mean; + bool scaling_median; double fstar_10; double alpha_star; @@ -84,9 +85,10 @@ struct HaloProperties{ void set_hbox_constants(double redshift, AstroParams *astro_params, FlagOptions *flag_options, struct HaloBoxConstants *consts){ consts->redshift = redshift; - //whether to fix *integrated* (not sampled) galaxy properties to the expected mean //Set on for the fixed grid case since we are missing halos above the cell mass consts->fix_mean = flag_options->FIXED_HALO_GRIDS; + //whether to fix *integrated* (not sampled) galaxy properties to the expected mean + consts->scaling_median = flag_options->HALO_SCALING_RELATIONS_MEDIAN; consts->fstar_10 = astro_params->F_STAR10; consts->alpha_star = astro_params->ALPHA_STAR; @@ -121,8 +123,8 @@ void set_hbox_constants(double redshift, AstroParams *astro_params, FlagOptions consts->mturn_m_nofb = 0.; if(flag_options->USE_MINI_HALOS){ - consts->mturn_m_nofb = lyman_werner_threshold(redshift, 0., consts->vcb_norel, astro_params); consts->vcb_norel = flag_options->FIX_VCB_AVG ? global_params.VAVG : 0; + consts->mturn_m_nofb = lyman_werner_threshold(redshift, 0., consts->vcb_norel, astro_params); } if(flag_options->FIXED_HALO_GRIDS || user_params_global->AVG_BELOW_SAMPLER){ @@ -207,8 +209,10 @@ double get_lx_on_sfr(double sfr, double metallicity, double lx_constant){ // return lx_on_sfr_Lehmer(metallicity); // return lx_on_sfr_Schechter(metallicity, lx_constant); // return lx_on_sfr_PL_Kaur(sfr,metallicity, lx_constant); - return lx_on_sfr_doublePL(metallicity, lx_constant); - // return lx_constant; + //HACK: new/old model switch with upperstellar flag + if(flag_options_global->USE_UPPER_STELLAR_TURNOVER) + return lx_on_sfr_doublePL(metallicity, lx_constant); + return lx_constant; } void get_halo_stellarmass(double halo_mass, double mturn_acg, double mturn_mcg, double star_rng, @@ -232,7 +236,8 @@ void get_halo_stellarmass(double halo_mass, double mturn_acg, double mturn_mcg, double sm_sample,sm_sample_mini; double baryon_ratio = cosmo_params_global->OMb / cosmo_params_global->OMm; - double stoc_adjustment_term = sigma_star * sigma_star / 2.; //adjustment to the mean for lognormal scatter + //adjustment to the mean for lognormal scatter + double stoc_adjustment_term = consts->scaling_median ? 0 : sigma_star * sigma_star / 2.; //We don't want an upturn even with a negative ALPHA_STAR if(flag_options_global->USE_UPPER_STELLAR_TURNOVER && (f_a > fu_a)){ @@ -278,7 +283,8 @@ void get_halo_sfr(double stellar_mass, double stellar_mass_mini, double sfr_rng, } sfr_mean = stellar_mass / (consts->t_star * consts->t_h); - double stoc_adjustment_term = sigma_sfr * sigma_sfr / 2.; //adjustment to the mean for lognormal scatter + //adjustment to the mean for lognormal scatter + double stoc_adjustment_term = consts->scaling_median ? 0 : sigma_sfr * sigma_sfr / 2.; sfr_sample = sfr_mean * exp(sfr_rng*sigma_sfr - stoc_adjustment_term); *sfr = sfr_sample; @@ -303,7 +309,9 @@ void get_halo_metallicity(double sfr, double stellar, double redshift, double *z void get_halo_xray(double sfr, double sfr_mini, double metallicity, double xray_rng, struct HaloBoxConstants *consts, double *xray_out){ double sigma_xray = consts->sigma_xray; - double stoc_adjustment_term = sigma_xray * sigma_xray / 2.; //adjustment to the mean for lognormal scatter + + //adjustment to the mean for lognormal scatter + double stoc_adjustment_term = consts->scaling_median ? 0 : sigma_xray * sigma_xray / 2.; double rng_factor = exp(xray_rng*consts->sigma_xray - stoc_adjustment_term); double lx_over_sfr = get_lx_on_sfr(sfr,metallicity,consts->l_x); @@ -877,7 +885,7 @@ void sum_halos_onto_grid(InitialConditions *ini_boxes, TsBox *previous_spin_temp } } total_n_halos = halos->n_halos - n_halos_cut; - LOG_SUPER_DEBUG("Cell 0 Halo: HM: %.2e SM: %.2e (%.2e) SF: %.2e (%.2e) X: %.2e NI: %.2e WS: %.2e ct : %d",grids->halo_mass[HII_R_INDEX(0,0,0)], + LOG_SUPER_DEBUG("Cell 0 Totals: HM: %.2e SM: %.2e (%.2e) SF: %.2e (%.2e) X: %.2e NI: %.2e WS: %.2e ct : %d",grids->halo_mass[HII_R_INDEX(0,0,0)], grids->halo_stars[HII_R_INDEX(0,0,0)],grids->halo_stars_mini[HII_R_INDEX(0,0,0)], grids->halo_sfr[HII_R_INDEX(0,0,0)],grids->halo_sfr_mini[HII_R_INDEX(0,0,0)], grids->halo_xray[HII_R_INDEX(0,0,0)],grids->n_ion[HII_R_INDEX(0,0,0)], diff --git a/src/py21cmfast/src/InitialConditions.c b/src/py21cmfast/src/InitialConditions.c index 4dfb859a..fab923e3 100644 --- a/src/py21cmfast/src/InitialConditions.c +++ b/src/py21cmfast/src/InitialConditions.c @@ -19,75 +19,10 @@ #include "indexing.h" #include "dft.h" #include "filtering.h" +#include "rng.h" #include "InitialConditions.h" -void seed_rng_threads(gsl_rng * rng_arr[], unsigned long long int seed){ - // setting tbe random seeds - gsl_rng * rseed = gsl_rng_alloc(gsl_rng_mt19937); // An RNG for generating seeds for multithreading - - gsl_rng_set(rseed, seed); - - unsigned int seeds[user_params_global->N_THREADS]; - - // For multithreading, seeds for the RNGs are generated from an initial RNG (based on the input random_seed) and then shuffled (Author: Fred Davies) - // int num_int = INT_MAX/16; - int num_int = INT_MAX/256; //JD: this was taking a few seconds per snapshot so i reduced the number TODO: init the RNG once - int i, thread_num; - unsigned int *many_ints = (unsigned int *)malloc((size_t)(num_int*sizeof(unsigned int))); // Some large number of possible integers - for (i=0; iN_THREADS, many_ints, num_int, sizeof(unsigned int)); // Populate the seeds array from the large list of integers - gsl_ran_shuffle(rseed, seeds, user_params_global->N_THREADS, sizeof(unsigned int)); // Shuffle the randomly selected integers - - int checker; - - checker = 0; - // seed the random number generators - for (thread_num = 0; thread_num < user_params_global->N_THREADS; thread_num++){ - switch (checker){ - case 0: - rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); - break; - case 1: - rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_gfsr4); - gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); - break; - case 2: - rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_cmrg); - gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); - break; - case 3: - rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_mrg); - gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); - break; - case 4: - rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_taus2); - gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); - break; - } // end switch - - checker += 1; - - if(checker==5) { - checker = 0; - } - } - - gsl_rng_free(rseed); - free(many_ints); -} - -void free_rng_threads(gsl_rng * rng_arr[]){ - int ii; - for(ii=0;iiN_THREADS;ii++){ - gsl_rng_free(rng_arr[ii]); - } -} - void adj_complex_conj(fftwf_complex *HIRES_box, UserParams *user_params, CosmoParams *cosmo_params){ /***** Adjust the complex conjugate relations for a real array *****/ diff --git a/src/py21cmfast/src/InitialConditions.h b/src/py21cmfast/src/InitialConditions.h index 181bd086..687bc64f 100644 --- a/src/py21cmfast/src/InitialConditions.h +++ b/src/py21cmfast/src/InitialConditions.h @@ -11,7 +11,4 @@ int ComputeInitialConditions( CosmoParams *cosmo_params, InitialConditions *boxes ); -void seed_rng_threads(gsl_rng * rng_arr[], unsigned long long int seed); -void free_rng_threads(gsl_rng * rng_arr[]); - #endif diff --git a/src/py21cmfast/src/IonisationBox.c b/src/py21cmfast/src/IonisationBox.c index 8a25c24b..5135928d 100644 --- a/src/py21cmfast/src/IonisationBox.c +++ b/src/py21cmfast/src/IonisationBox.c @@ -27,1477 +27,1366 @@ #include "interp_tables.h" #include "filtering.h" #include "bubble_helper_progs.h" +#include "InitialConditions.h" #include "IonisationBox.h" #define LOG10_MTURN_MAX ((double)(10)) //maximum mturn limit enforced on grids -int INIT_ERFC_INTERPOLATION = 1; int INIT_RECOMBINATIONS = 1; -double *ERFC_VALS, *ERFC_VALS_DIFF; - -float absolute_delta_z; - -float thistk; - -int ComputeIonizedBox(float redshift, float prev_redshift, UserParams *user_params, CosmoParams *cosmo_params, - AstroParams *astro_params, FlagOptions *flag_options, - PerturbedField *perturbed_field, - PerturbedField *previous_perturbed_field, - IonizedBox *previous_ionize_box, - TsBox *spin_temp, - HaloBox *halos, - InitialConditions *ini_boxes, - IonizedBox *box) { - - int status; - - Try{ // This Try brackets the whole function, so we don't indent. - LOG_DEBUG("input values:"); - LOG_DEBUG("redshift=%f, prev_redshift=%f", redshift, prev_redshift); - #if LOG_LEVEL >= DEBUG_LEVEL - writeUserParams(user_params); - writeCosmoParams(cosmo_params); - writeAstroParams(flag_options, astro_params); - writeFlagOptions(flag_options); - #endif - - // Makes the parameter structs visible to a variety of functions/macros - // Do each time to avoid Python garbage collection issues - Broadcast_struct_global_all(user_params,cosmo_params,astro_params,flag_options); - - omp_set_num_threads(user_params->N_THREADS); +//Parameters for the ionisation box calculations +struct IonBoxConstants{ + //redshift constants + double redshift; + double stored_redshift; + double prev_redshift; + double growth_factor; + double prev_growth_factor; + double photoncons_adjustment_factor; + double dz; + double fabs_dtdz; + + //compound/alternate flags + bool fix_mean; + bool filter_recombinations; + + //astro parameters + double fstar_10; + double alpha_star; + double fstar_7; + double alpha_star_mini; + double t_h; + double t_star_sec; + double fesc_10; + double alpha_esc; + double fesc_7; + + //astro calculated values + double vcb_norel; + double mturn_a_nofb; + double mturn_m_nofb; + double ion_eff_factor; + double ion_eff_factor_mini; + double ion_eff_factor_gl; + double ion_eff_factor_mini_gl; + double mfp_meandens; + double gamma_prefactor; + double gamma_prefactor_mini; + + double TK_nofluct; + double adia_TK_term; + + //power-law limits + double Mlim_Fstar; + double Mlim_Fesc; + double Mlim_Fstar_mini; + double Mlim_Fesc_mini; + + //HMF limits + double M_min; + double lnMmin; + double M_max_gl; + double lnMmax_gl; + double sigma_minmass; + + //dimensions + double pixel_length; + double pixel_mass; +}; + +//TODO: Something like these structs should also replace the globals in SpinTemperature.c + +//TODO: consider the case of having Radiusspec as array of structs (current) +// vs struct of arrays. The former allows you to pass a single struct into each function +// without an index so they don't need to know about other radii. +// The second is simpler to understand/work with in terms of allocation/scoping +struct RadiusSpec{ + //calculated and stored at the beginning + double R; + double M_max_R; + double ln_M_max_R; + double sigma_maxmass; + int R_index; + + //calculated within the loop + double f_coll_grid_mean; + double f_coll_grid_mean_MINI; +}; + +//this struct should hold all the pointers to grids which need to be filtered +struct FilteredGrids{ + //Always Used + fftwf_complex *deltax_unfiltered, *deltax_filtered; + + //Used when TS_FLUCT==True + fftwf_complex *xe_unfiltered, *xe_filtered; + + //Used when INHOMO_RECO==True and CELL_RECOMB=False + fftwf_complex *N_rec_unfiltered, *N_rec_filtered; + + //Used when USE_MINI_HALOS==True and USE_HALO_FIELD=False + fftwf_complex *prev_deltax_unfiltered, *prev_deltax_filtered; + fftwf_complex *log10_Mturnover_unfiltered, *log10_Mturnover_filtered; + fftwf_complex *log10_Mturnover_MINI_unfiltered, *log10_Mturnover_MINI_filtered; + + //Used when USE_HALO_FIELD=True + fftwf_complex *stars_unfiltered, *stars_filtered; + fftwf_complex *sfr_unfiltered, *sfr_filtered; +}; + +void set_ionbox_constants(double redshift, double prev_redshift, CosmoParams *cosmo_params, AstroParams *astro_params, + FlagOptions *flag_options, struct IonBoxConstants *consts){ + consts->redshift = redshift; + //defaults for no photoncons + consts->stored_redshift = redshift; + consts->photoncons_adjustment_factor = 1.; + + //dz is only used if inhomo_reco + if (prev_redshift < 1) + consts->dz = (1. + redshift) * (global_params.ZPRIME_STEP_FACTOR - 1.); + else + consts->dz = prev_redshift - redshift; + + //TODO: Figure out why we have the 1e15 here + consts->fabs_dtdz = fabs(dtdz(redshift))/1e15; //reduce to have good precision + + consts->growth_factor = dicke(redshift); + consts->prev_growth_factor = dicke(prev_redshift); + //whether to fix *integrated* (not sampled) galaxy properties to the expected mean + // constant for now, to be a flag later + consts->fix_mean = !flag_options->USE_HALO_FIELD; + consts->filter_recombinations = flag_options->INHOMO_RECO && !flag_options->CELL_RECOMB; + + consts->fstar_10 = astro_params->F_STAR10; + consts->alpha_star = astro_params->ALPHA_STAR; + + consts->fstar_7 = astro_params->F_STAR7_MINI; + consts->alpha_star_mini = astro_params->ALPHA_STAR_MINI; + + consts->t_h = t_hubble(redshift); + consts->t_star_sec = astro_params->t_STAR * consts->t_h; + + consts->alpha_esc = astro_params->ALPHA_ESC; + consts->fesc_10= astro_params->F_ESC10; + consts->fesc_7 = astro_params->F_ESC7_MINI; + + if(flag_options->PHOTON_CONS_TYPE == 2) + consts->alpha_esc = get_fesc_fit(redshift); + else if(flag_options->PHOTON_CONS_TYPE == 3) + consts->fesc_10 = get_fesc_fit(redshift); + + consts->mturn_a_nofb = flag_options->USE_MINI_HALOS ? atomic_cooling_threshold(redshift) : astro_params->M_TURN; + + consts->mturn_m_nofb = 0.; + if(flag_options->USE_MINI_HALOS){ + consts->vcb_norel = flag_options->FIX_VCB_AVG ? global_params.VAVG : 0; + consts->mturn_m_nofb = lyman_werner_threshold(redshift, 0., consts->vcb_norel, astro_params); + } - // Other parameters used in the code - int i,j,k,x,y,z, LAST_FILTER_STEP, first_step_R; - int counter, N_halos_in_cell; - unsigned long long ct; + if(consts->mturn_m_nofb < astro_params->M_TURN)consts->mturn_m_nofb = astro_params->M_TURN; + if(consts->mturn_a_nofb < astro_params->M_TURN)consts->mturn_a_nofb = astro_params->M_TURN; - float growth_factor, pixel_mass, cell_length_factor, M_MIN, prev_growth_factor; - float erfc_denom, res_xH, Splined_Fcoll, xHII_from_xrays, curr_dens, massofscaleR, ION_EFF_FACTOR; - float Splined_Fcoll_MINI, prev_dens, ION_EFF_FACTOR_MINI, prev_Splined_Fcoll, prev_Splined_Fcoll_MINI; - float ave_M_coll_cell, ave_N_min_cell; - double lnMmin,lnMmax,lnM_cond; + if(flag_options->FIXED_HALO_GRIDS || user_params_global->AVG_BELOW_SAMPLER){ + consts->Mlim_Fstar = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, consts->alpha_star, consts->fstar_10); + consts->Mlim_Fesc = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, consts->alpha_esc, consts->fesc_10); - float curr_vcb; + if(flag_options->USE_MINI_HALOS){ + consts->Mlim_Fstar_mini = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, consts->alpha_star_mini, + consts->fstar_7 * pow(1e3,consts->alpha_star_mini)); + consts->Mlim_Fesc_mini = Mass_limit_bisection(global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL, consts->alpha_esc, + consts->fesc_7 * pow(1e3,consts->alpha_esc)); + } + } - double global_xH, ST_over_PS, f_coll, R, stored_R, f_coll_min; - double ST_over_PS_MINI, f_coll_MINI, f_coll_min_MINI; + if(flag_options->USE_MASS_DEPENDENT_ZETA) { + consts->ion_eff_factor_gl = global_params.Pop2_ion * astro_params->F_STAR10 * consts->fesc_10; + consts->ion_eff_factor_mini_gl = global_params.Pop3_ion * astro_params->F_STAR7_MINI * astro_params->F_ESC7_MINI; + } + else { + consts->ion_eff_factor_gl = astro_params->HII_EFF_FACTOR; + consts->ion_eff_factor_mini_gl = 0.; + } - double t_ast, Gamma_R_prefactor, rec, dNrec, sigmaMmax; - double Gamma_R_prefactor_MINI; - float fabs_dtdz, ZSTEP, z_eff; + //The halo fields already have Fstar,Fesc,nion taken into account, so their global factor differs from the local one + if(flag_options->USE_HALO_FIELD) { + consts->ion_eff_factor = 1.; + consts->ion_eff_factor_mini = 1.; + } + else{ + consts->ion_eff_factor = consts->ion_eff_factor_gl; + consts->ion_eff_factor_mini = consts->ion_eff_factor_mini_gl; + } - int something_finite_or_infinite = 0; - int *overdense_int_boundexceeded_threaded = calloc(user_params->N_THREADS,sizeof(int)); + //MFP USED FOR THE EXPNENTIAL FILTER + //Yuxiang's evolving Rmax for MFP in ionised regions fit from Songaila+2010 + // if(flag_options->USE_EXP_FILTER){ + // if (redshift > 6) + // consts->mfp_meandens = 25.483241248322766 / cosmo_params->hlittle; + // else + // consts->mfp_meandens = 112 / cosmo_params->hlittle * pow( (1.+redshift) / 5. , -4.4); + // LOG_DEBUG("Set mfp = %.4e",consts->mfp_meandens); + // } + //Constant MFP + consts->mfp_meandens = 25.483241248322766 / cosmo_params->hlittle; - double ave_log10_Mturnover, ave_log10_Mturnover_MINI; + //set the minimum source mass + consts->M_min = minimum_source_mass(redshift,false,astro_params,flag_options); + consts->lnMmin = log(consts->M_min); + consts->lnMmax_gl = log(global_params.M_MAX_INTEGRAL); + consts->sigma_minmass = sigma_z0(consts->M_min); + + //global TK and adiabatic terms for temperature without the Ts Calculation + //final temperature = TK * (1+cT_ad*delta) + if(!flag_options->USE_TS_FLUCT){ + consts->TK_nofluct = T_RECFAST(redshift,0); + //finding the adiabatic index at the initial redshift from 2302.08506 to fix adiabatic fluctuations. + consts->adia_TK_term = cT_approx(redshift); + } + consts->pixel_length = user_params_global->BOX_LEN/(double)user_params_global->HII_DIM; + consts->pixel_mass = cosmo_params_global->OMm * RHOcrit * pow(consts->pixel_length,3); + + consts->gamma_prefactor = pow(1+redshift,2) * CMperMPC * SIGMA_HI * global_params.ALPHA_UVB \ + / (global_params.ALPHA_UVB+2.75) * N_b0 * consts->ion_eff_factor / 1.0e-12; + if(flag_options->USE_HALO_FIELD) + consts->gamma_prefactor /= RHOcrit * cosmo_params->OMb; //TODO: double-check these unit differences, HaloBox.halo_wsfr vs Nion_General units + else + consts->gamma_prefactor /= consts->t_star_sec; + consts->gamma_prefactor_mini = consts->gamma_prefactor * consts->ion_eff_factor_mini / consts->ion_eff_factor; + LOG_SUPER_DEBUG("Gamma Prefactor %.3e ion eff factor %.3e",consts->gamma_prefactor,consts->ion_eff_factor); + LOG_SUPER_DEBUG("Mini Gamma %.3e Mini ion %.3e",consts->gamma_prefactor_mini,consts->ion_eff_factor_mini); +} - float Mlim_Fstar, Mlim_Fesc; - float Mlim_Fstar_MINI, Mlim_Fesc_MINI; - float Mcrit_atom, log10_Mcrit_atom, log10_Mcrit_mol; - fftwf_complex *log10_Mturnover_unfiltered=NULL, *log10_Mturnover_filtered=NULL; - fftwf_complex *log10_Mturnover_MINI_unfiltered=NULL, *log10_Mturnover_MINI_filtered=NULL; - float log10_Mturnover, log10_Mturnover_MINI, Mcrit_LW, Mcrit_RE, Mturnover, Mturnover_MINI; +void allocate_fftw_grids(struct FilteredGrids **fg_struct){ + *fg_struct = malloc(sizeof(**fg_struct)); - float min_density, max_density; - float prev_min_density, prev_max_density; + (*fg_struct)->deltax_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->deltax_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - float stored_redshift, adjustment_factor; + if(flag_options_global->USE_MINI_HALOS && !flag_options_global->USE_HALO_FIELD){ + (*fg_struct)->prev_deltax_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->prev_deltax_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - float log10Mturn_min, log10Mturn_max; - float log10Mturn_min_MINI, log10Mturn_max_MINI; + (*fg_struct)->log10_Mturnover_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->log10_Mturnover_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->log10_Mturnover_MINI_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->log10_Mturnover_MINI_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + } - gsl_rng * r[user_params->N_THREADS]; + if(flag_options_global->USE_TS_FLUCT){ + (*fg_struct)->xe_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->xe_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + } -LOG_SUPER_DEBUG("initing heat"); - init_heat(); - float TK; - TK = T_RECFAST(redshift,0); - float cT_ad; //finding the adiabatic index at the initial redshift from 2302.08506 to fix adiabatic fluctuations. - cT_ad = cT_approx(redshift); -LOG_SUPER_DEBUG("inited heat"); + if(flag_options_global->INHOMO_RECO && !flag_options_global->CELL_RECOMB){ + (*fg_struct)->N_rec_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); // cumulative number of recombinations + (*fg_struct)->N_rec_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + } - init_ps(); + if(flag_options_global->USE_HALO_FIELD){ + (*fg_struct)->stars_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->stars_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); -LOG_SUPER_DEBUG("defined parameters"); + (*fg_struct)->sfr_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + (*fg_struct)->sfr_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + } + //TODO: check for null pointers and throw errors +} - double alpha_esc_var = astro_params->ALPHA_ESC; - double norm_esc_var = astro_params->F_ESC10; +void free_fftw_grids(struct FilteredGrids *fg_struct){ + fftwf_free(fg_struct->deltax_unfiltered); + fftwf_free(fg_struct->deltax_filtered); - if(flag_options->PHOTON_CONS_TYPE == 2){ - alpha_esc_var = get_fesc_fit(redshift); - LOG_DEBUG("PHOTONCONS: ALPHA ESC %.2e --> %.2e",astro_params->ALPHA_ESC,alpha_esc_var); - } - else if(flag_options->PHOTON_CONS_TYPE == 3){ - norm_esc_var = get_fesc_fit(redshift); - LOG_DEBUG("PHOTONCONS: F_ESC10 %.2e --> %.2e",astro_params->F_ESC10,norm_esc_var); - } + if(flag_options_global->USE_MINI_HALOS && !flag_options_global->USE_HALO_FIELD){ + fftwf_free(fg_struct->prev_deltax_unfiltered); + fftwf_free(fg_struct->prev_deltax_filtered); - //escape fractions taken into account in halo field - if(flag_options->USE_HALO_FIELD){ - ION_EFF_FACTOR = 1.; - ION_EFF_FACTOR_MINI = 1.; + fftwf_free(fg_struct->log10_Mturnover_unfiltered); + fftwf_free(fg_struct->log10_Mturnover_filtered); + fftwf_free(fg_struct->log10_Mturnover_MINI_unfiltered); + fftwf_free(fg_struct->log10_Mturnover_MINI_filtered); } - else if(flag_options->USE_MASS_DEPENDENT_ZETA) { - ION_EFF_FACTOR = global_params.Pop2_ion * astro_params->F_STAR10 * norm_esc_var; - ION_EFF_FACTOR_MINI = global_params.Pop3_ion * astro_params->F_STAR7_MINI * astro_params->F_ESC7_MINI; + if(flag_options_global->USE_TS_FLUCT) { + fftwf_free(fg_struct->xe_unfiltered); + fftwf_free(fg_struct->xe_filtered); } - else { - ION_EFF_FACTOR = astro_params->HII_EFF_FACTOR; - ION_EFF_FACTOR_MINI = 0.; + if (flag_options_global->INHOMO_RECO && !flag_options_global->CELL_RECOMB){ + fftwf_free(fg_struct->N_rec_unfiltered); + fftwf_free(fg_struct->N_rec_filtered); } - //Yuxiang's evolving Rmax for MFP in ionised regions - double exp_mfp; - // if(flag_options->USE_EXP_FILTER){ - // if (redshift > 6) - // exp_mfp = 25.483241248322766 / cosmo_params->hlittle; - // else - // exp_mfp = 112 / cosmo_params->hlittle * pow( (1.+redshift) / 5. , -4.4); - // } - //constant for now, the pow(4.4) makes kinks in reionisation - exp_mfp = 25.483241248322766 / cosmo_params->hlittle; - - // For recombinations - bool recomb_filter_flag = flag_options->INHOMO_RECO && !flag_options->CELL_RECOMB; - if(flag_options->INHOMO_RECO) { - if(INIT_RECOMBINATIONS) { - init_MHR(); - INIT_RECOMBINATIONS=0; - } + if(flag_options_global->USE_HALO_FIELD) { + fftwf_free(fg_struct->stars_unfiltered); + fftwf_free(fg_struct->stars_filtered); + fftwf_free(fg_struct->sfr_unfiltered); + fftwf_free(fg_struct->sfr_filtered); + } - if (prev_redshift < 1) //deal with first redshift - ZSTEP = (1. + redshift) * (global_params.ZPRIME_STEP_FACTOR - 1.); - else - ZSTEP = prev_redshift - redshift; + free(fg_struct); +} -#pragma omp parallel shared(box) private(ct) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (ct=0; ctGamma12_box[ct] = 0.0; - box->MFP_box[ct] = 0.0; +//fill fftwf boxes, do the r2c transform and normalise +//TODO: check for some invalid limit values to skip the clipping step +void prepare_box_for_filtering(float *input_box, fftwf_complex *output_c_box, double const_factor, double limit_min, double limit_max){ + int i,j,k; + unsigned long long int ct; + double curr_cell; + //NOTE: Meraxes just applies a pointer cast box = (fftwf_complex *) input. Figure out why this works, + // They pad the input by a factor of 2 to cover the complex part, but from the type I thought it would be stored [(r,c),(r,c)...] + // Not [(r,r,r,r....),(c,c,c....)] so the alignment should be wrong, right? + #pragma omp parallel for private(i,j,k,curr_cell) num_threads(user_params_global->N_THREADS) collapse(3) + for(i=0;iHII_DIM;i++){ + for(j=0;jHII_DIM;j++){ + for(k=0;kUSE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, output_c_box); + + //divide by pixel number in preparation for later inverse transform + #pragma omp parallel for private(ct) num_threads(user_params_global->N_THREADS) + for (ct=0; ctN_THREADS) +//Populating a previous box which has the required fields for the first snapshot: +// Since the values of all arrays should be passed in as zero, we only assign special +// first snapshot values here +void setup_first_z_prevbox(IonizedBox *previous_ionize_box, PerturbedField *previous_perturb, int n_radii){ + LOG_DEBUG("first redshift, do some initialization"); + int i,j,k; + unsigned long long int ct; + + //z_re_box is used in all cases + #pragma omp parallel shared(previous_ionize_box) private(i,j,k) num_threads(user_params_global->N_THREADS) { -#pragma omp for - for (ct=0; ctz_re_box[ct] = -1.0; + #pragma omp for + for (i=0; iHII_DIM; i++){ + for (j=0; jHII_DIM; j++){ + for (k=0; kz_re_box[HII_R_INDEX(i, j, k)] = -1.0; + } + } } } - LOG_SUPER_DEBUG("z_re_box init: "); - debugSummarizeBox(box->z_re_box, user_params->HII_DIM, user_params->NON_CUBIC_FACTOR, " "); - - fabs_dtdz = fabs(dtdz(redshift))/1e15; //reduce to have good precision - t_ast = astro_params->t_STAR * t_hubble(redshift); - - // Modify the current sampled redshift to a redshift which matches the expected filling factor given our astrophysical parameterisation. - // This is the photon non-conservation correction - if(flag_options->PHOTON_CONS_TYPE == 1) { - adjust_redshifts_for_photoncons(astro_params,flag_options,&redshift,&stored_redshift,&absolute_delta_z); - LOG_DEBUG("PhotonCons data:"); - LOG_DEBUG("original redshift=%f, updated redshift=%f delta-z = %f", stored_redshift, redshift, absolute_delta_z); - if(isfinite(redshift)==0 || isfinite(absolute_delta_z)==0) { - LOG_ERROR("Updated photon non-conservation redshift is either infinite or NaN!"); - LOG_ERROR("This can sometimes occur when reionisation stalls (i.e. extremely low"\ - "F_ESC or F_STAR or not enough sources)"); - Throw(PhotonConsError); + //previous Gamma12 is used for reionisation feedback when USE_MINI_HALOS + //previous delta and Fcoll are used for the trapezoidal integral when USE_MINI_HALOS + if(flag_options_global->USE_MINI_HALOS){ + previous_ionize_box->mean_f_coll = 0.0; + previous_ionize_box->mean_f_coll_MINI = 0.0; + #pragma omp parallel private(ct) num_threads(user_params_global->N_THREADS) + { + #pragma omp for + for(ct=0;ctdensity[ct] = -1.5; //TODO: figure out why 1.5, ensures no sources? + } } } +} - Splined_Fcoll = 0.; - Splined_Fcoll_MINI = 0.; +void calculate_mcrit_boxes(IonizedBox *prev_ionbox, TsBox *spin_temp, InitialConditions *ini_boxes, struct IonBoxConstants *consts, + fftwf_complex *log10_mcrit_acg, fftwf_complex *log10_mcrit_mcg, double *avg_mturn_acg, double *avg_mturn_mcg){ + double ave_log10_Mturnover = 0.; + double ave_log10_Mturnover_MINI = 0.; - double ArgBinWidth, InvArgBinWidth, erfc_arg_val, erfc_arg_min, erfc_arg_max; - int erfc_arg_val_index, ERFC_NUM_POINTS; + #pragma omp parallel num_threads(user_params_global->N_THREADS) + { + int x,y,z; + double Mcrit_RE, Mcrit_LW; + double curr_Mt, curr_Mt_MINI; + double curr_vcb = consts->vcb_norel; + #pragma omp for reduction(+:ave_log10_Mturnover,ave_log10_Mturnover_MINI) + for (x=0; xHII_DIM; x++){ + for (y=0; yHII_DIM; y++){ + for (z=0; zredshift, prev_ionbox->Gamma12_box[HII_R_INDEX(x, y, z)], prev_ionbox->z_re_box[HII_R_INDEX(x, y, z)]); - erfc_arg_val = 0.; - erfc_arg_val_index = 0; + if(user_params_global->USE_RELATIVE_VELOCITIES && !flag_options_global->FIX_VCB_AVG) + curr_vcb = ini_boxes->lowres_vcb[HII_R_INDEX(x,y,z)]; - // Setup an interpolation table for the error function, helpful for calcluating the collapsed fraction - // (only for the default model, i.e. mass-independent ionising efficiency) - erfc_arg_min = -15.0; - erfc_arg_max = 15.0; + Mcrit_LW = lyman_werner_threshold(consts->redshift, spin_temp->J_21_LW_box[HII_R_INDEX(x, y, z)], curr_vcb, astro_params_global); + if(Mcrit_LW != Mcrit_LW || Mcrit_LW == 0){ + LOG_ERROR("Mcrit error %d %d %d: M %.2e z %.2f J %.2e v %.2e",x,y,z,Mcrit_LW,consts->redshift, + spin_temp->J_21_LW_box[HII_R_INDEX(x, y, z)],curr_vcb); + Throw(ValueError); + } - ERFC_NUM_POINTS = 10000; + //JBM: this only accounts for effect 3 (largest on minihaloes). Effects 1 and 2 affect both minihaloes (MCGs) and regular ACGs, but they're smaller ~10%. See Sec 2 of Muñoz+21 (2110.13919) + curr_Mt = log10(fmax(Mcrit_RE,consts->mturn_a_nofb)); + curr_Mt_MINI = log10(fmax(Mcrit_RE,fmax(Mcrit_LW,consts->mturn_m_nofb))); - ArgBinWidth = (erfc_arg_max - erfc_arg_min)/((double)ERFC_NUM_POINTS - 1.); - InvArgBinWidth = 1./ArgBinWidth; + //To avoid allocating another box we directly assign turnover masses to the fftw grid + *((float *)log10_mcrit_acg + HII_R_FFT_INDEX(x,y,z)) = curr_Mt; + *((float *)log10_mcrit_mcg + HII_R_FFT_INDEX(x,y,z)) = curr_Mt_MINI; - if(!flag_options->USE_MASS_DEPENDENT_ZETA && INIT_ERFC_INTERPOLATION) { + ave_log10_Mturnover += curr_Mt; + ave_log10_Mturnover_MINI += curr_Mt_MINI; + } + } + } + } + *avg_mturn_acg = ave_log10_Mturnover/HII_TOT_NUM_PIXELS; + *avg_mturn_mcg = ave_log10_Mturnover_MINI/HII_TOT_NUM_PIXELS; +} - ERFC_VALS = calloc(ERFC_NUM_POINTS,sizeof(double)); - ERFC_VALS_DIFF = calloc(ERFC_NUM_POINTS,sizeof(double)); +// Determine the normalisation for the excursion set algorithm +// When USE_MINI_HALOS==True, we do a trapezoidal integration, where we take +// F_coll = f(z_current,Mturn_current) - f(z_previous,Mturn_current) + f(z_previous,Mturn_previous) +// all mturns are average log10 over the +// the `limit` outputs are set to the total value at the maximum redshift and current turnover, these form a +// lower limit on any grid cell +void set_mean_fcoll(struct IonBoxConstants *c, IonizedBox *prev_box, IonizedBox *curr_box, double mturn_acg, double mturn_mcg, double *f_limit_acg, double *f_limit_mcg){ + double f_coll_curr, f_coll_prev, f_coll_curr_mini, f_coll_prev_mini; + if(flag_options_global->USE_MASS_DEPENDENT_ZETA){ + f_coll_curr = Nion_General(c->redshift,c->lnMmin,c->lnMmax_gl,mturn_acg,c->alpha_star,c->alpha_esc, + c->fstar_10,c->fesc_10,c->Mlim_Fstar,c->Mlim_Fesc); + *f_limit_acg = Nion_General(global_params.Z_HEAT_MAX,c->lnMmin,c->lnMmax_gl,mturn_acg,c->alpha_star,c->alpha_esc, + c->fstar_10,c->fesc_10,c->Mlim_Fstar,c->Mlim_Fesc); + + if (flag_options_global->USE_MINI_HALOS){ + if (prev_box->mean_f_coll * c->ion_eff_factor_gl < 1e-4){ + //we don't have enough ionising radiation in the previous snapshot, just take the current value + curr_box->mean_f_coll = f_coll_curr; + } + else{ + f_coll_prev = Nion_General(c->prev_redshift,c->lnMmin,c->lnMmax_gl,mturn_acg,c->alpha_star,c->alpha_esc, + c->fstar_10,c->fesc_10,c->Mlim_Fstar,c->Mlim_Fesc); + curr_box->mean_f_coll = prev_box->mean_f_coll + f_coll_curr - f_coll_prev; + } + f_coll_curr_mini = Nion_General_MINI(c->redshift,c->lnMmin,c->lnMmax_gl,mturn_mcg,mturn_acg,c->alpha_star_mini, + c->alpha_esc,c->fstar_7,c->fesc_7,c->Mlim_Fstar_mini,c->Mlim_Fesc_mini); + if (prev_box->mean_f_coll_MINI * c->ion_eff_factor_gl < 1e-4){ + curr_box->mean_f_coll_MINI = f_coll_curr_mini; + } + else{ + f_coll_prev_mini = Nion_General_MINI(c->prev_redshift,c->lnMmin,c->lnMmax_gl,mturn_mcg,mturn_acg,c->alpha_star_mini, + c->alpha_esc,c->fstar_7,c->fesc_7,c->Mlim_Fstar_mini,c->Mlim_Fesc_mini); + curr_box->mean_f_coll_MINI = prev_box->mean_f_coll_MINI + f_coll_curr_mini - f_coll_prev_mini; + } + *f_limit_mcg = Nion_General_MINI(global_params.Z_HEAT_MAX,c->lnMmin,c->lnMmax_gl,mturn_mcg,mturn_acg,c->alpha_star_mini, + c->alpha_esc,c->fstar_7,c->fesc_7, + c->Mlim_Fstar_mini,c->Mlim_Fesc_mini); + } + else{ + curr_box->mean_f_coll = f_coll_curr; + curr_box->mean_f_coll_MINI = 0.; + } + } + else { + curr_box->mean_f_coll = Fcoll_General(c->redshift, c->lnMmin, c->lnMmax_gl); + *f_limit_acg = Fcoll_General(global_params.Z_HEAT_MAX, c->lnMmin, c->lnMmax_gl); //JD: the old parametrisation didn't have this limit before + } -#pragma omp parallel shared(ERFC_VALS,erfc_arg_min,ArgBinWidth) private(i,erfc_arg_val) num_threads(user_params->N_THREADS) - { -#pragma omp for - for(i=0;imean_f_coll)==0){ + LOG_ERROR("Mean collapse fraction is either infinite or NaN!"); + Throw(InfinityorNaNError); + } + LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll: %e", curr_box->mean_f_coll); - erfc_arg_val = erfc_arg_min + ArgBinWidth*(double)i; + 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!"); + Throw(InfinityorNaNError); + } + LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll_MINI: %e", curr_box->mean_f_coll_MINI); + } +} - ERFC_VALS[i] = splined_erfc(erfc_arg_val); +double set_fully_neutral_box(IonizedBox *box, TsBox *spin_temp, PerturbedField *perturbed_field, struct IonBoxConstants *consts){ + double global_xH; + unsigned long long int ct; + if(flag_options_global->USE_TS_FLUCT) { + #pragma omp parallel private(ct) num_threads(user_params_global->N_THREADS) + { + #pragma omp for reduction(+:global_xH) + for (ct=0; ctxH_box[ct] = 1.-spin_temp->x_e_box[ct]; // convert from x_e to xH + global_xH += box->xH_box[ct]; + box->temp_kinetic_all_gas[ct] = spin_temp->Tk_box[ct]; } } - -#pragma omp parallel shared(ERFC_VALS_DIFF,ERFC_VALS) private(i) num_threads(user_params->N_THREADS) + global_xH /= (double)HII_TOT_NUM_PIXELS; + } + else { + global_xH = 1. - xion_RECFAST(consts->redshift, 0); + #pragma omp parallel private(ct) num_threads(user_params_global->N_THREADS) { -#pragma omp for - for(i=0;i<(ERFC_NUM_POINTS-1);i++) { - ERFC_VALS_DIFF[i] = ERFC_VALS[i+1] - ERFC_VALS[i]; + #pragma omp for + for (ct=0; ctxH_box[ct] = global_xH; + box->temp_kinetic_all_gas[ct] = consts->TK_nofluct * (1.0 + consts->adia_TK_term * perturbed_field->density[ct]); // Is perturbed_field defined already here? we need it for cT. I'm also assuming we don't need to multiply by other z here. } } - - INIT_ERFC_INTERPOLATION = 0; } + return global_xH; +} -LOG_SUPER_DEBUG("erfc interpolation done"); - - ///////////////////////////////// BEGIN INITIALIZATION ////////////////////////////////// - - // perform a very rudimentary check to see if we are underresolved and not using the linear approx - if ((user_params->BOX_LEN > user_params->DIM) && !(global_params.EVOLVE_DENSITY_LINEARLY)){ - LOG_WARNING("Resolution is likely too low for accurate evolved density fields\n It Is recommended \ - that you either increase the resolution (DIM/Box_LEN) or set the EVOLVE_DENSITY_LINEARLY flag to 1\n"); +//TODO: SPEED TEST THE FOLLOWING ORDERS: +// (copy,copy,copy....) (filter,filter,filter,...) (transform,transform,...) +// (copy,filter,transform), (copy,filter,transform), (copy,filter,transform)... +// if the first is faster, consider reordering prepare_filter_grids() in the same way +// if the second is faster, make a function to do this for one grid +void copy_filter_transform(struct FilteredGrids *fg_struct, struct IonBoxConstants *consts, struct RadiusSpec rspec){ + memcpy(fg_struct->deltax_filtered, fg_struct->deltax_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + if(flag_options_global->USE_TS_FLUCT){ + memcpy(fg_struct->xe_filtered, fg_struct->xe_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); + } + if(consts->filter_recombinations){ + memcpy(fg_struct->N_rec_filtered, fg_struct->N_rec_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); + } + if(flag_options_global->USE_HALO_FIELD){ + memcpy(fg_struct->stars_filtered, fg_struct->stars_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); + memcpy(fg_struct->sfr_filtered, fg_struct->sfr_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); + } + else{ + if(flag_options_global->USE_MINI_HALOS){ + memcpy(fg_struct->prev_deltax_filtered, fg_struct->prev_deltax_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + memcpy(fg_struct->log10_Mturnover_MINI_filtered, fg_struct->log10_Mturnover_MINI_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + memcpy(fg_struct->log10_Mturnover_filtered, fg_struct->log10_Mturnover_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + } } - // initialize power spectrum - growth_factor = dicke(redshift); - prev_growth_factor = dicke(prev_redshift); - - fftwf_complex *deltax_unfiltered, *deltax_unfiltered_original, *deltax_filtered; - fftwf_complex *xe_unfiltered, *xe_filtered, *N_rec_unfiltered, *N_rec_filtered; - fftwf_complex *prev_deltax_unfiltered, *prev_deltax_filtered; - - //new halo property grids - fftwf_complex *stars_unfiltered,*stars_filtered; - fftwf_complex *sfr_unfiltered,*sfr_filtered; - - //NOTE FOR REFACTOR: These don't need to be allocated/filtered if (USE_HALO FIELD && CELL_RECOMB) - // Also, I don't think deltax_unfiltered_original is useful at all since we have the PerturbedField - deltax_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - deltax_unfiltered_original = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - deltax_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - - if (flag_options->USE_MINI_HALOS){ - prev_deltax_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - prev_deltax_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + if(rspec.R_index > 0){ + double R = rspec.R; + filter_box(fg_struct->deltax_filtered, 1, global_params.HII_FILTER, R, 0.); + if (flag_options_global->USE_TS_FLUCT) { + filter_box(fg_struct->xe_filtered, 1, global_params.HII_FILTER, R, 0.); + } + if (consts->filter_recombinations) { + filter_box(fg_struct->N_rec_filtered, 1, global_params.HII_FILTER, R, 0.); + } + if (flag_options_global->USE_HALO_FIELD) { + int filter_hf = flag_options_global->USE_EXP_FILTER ? 3 : global_params.HII_FILTER; + filter_box(fg_struct->stars_filtered, 1, filter_hf, R, consts->mfp_meandens); + filter_box(fg_struct->sfr_filtered, 1, filter_hf, R, consts->mfp_meandens); + } + else{ + if(flag_options_global->USE_MINI_HALOS){ + filter_box(fg_struct->prev_deltax_filtered, 1, global_params.HII_FILTER, R, 0.); + filter_box(fg_struct->log10_Mturnover_MINI_filtered, 1, global_params.HII_FILTER, R, 0.); + filter_box(fg_struct->log10_Mturnover_filtered, 1, global_params.HII_FILTER, R, 0.); + } + } } - if(flag_options->USE_TS_FLUCT) { - xe_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - xe_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + // Perform FFTs + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->deltax_filtered); + if (flag_options_global->USE_HALO_FIELD) { + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->stars_filtered); + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->sfr_filtered); } - if (recomb_filter_flag){ - N_rec_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); // cumulative number of recombinations - N_rec_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + else{ + if(flag_options_global->USE_MINI_HALOS){ + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->prev_deltax_filtered); + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->log10_Mturnover_MINI_filtered); + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->log10_Mturnover_filtered); + } } - - // Calculate the density field for this redshift if the initial conditions/cosmology are changing - if(flag_options->PHOTON_CONS_TYPE == 1) { - adjustment_factor = dicke(redshift)/dicke(stored_redshift); + if (flag_options_global->USE_TS_FLUCT) { + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->xe_filtered); } - else { - adjustment_factor = 1.; + if (consts->filter_recombinations) { + dft_c2r_cube(user_params_global->USE_FFTW_WISDOM, user_params_global->HII_DIM, HII_D_PARA, user_params_global->N_THREADS, fg_struct->N_rec_filtered); } +} -#pragma omp parallel shared(deltax_unfiltered,perturbed_field,adjustment_factor) private(i,j,k) num_threads(user_params->N_THREADS) +//After filtering the grids, we need to clip them to physical values and take the extrema for some interpolation tables +void clip_and_get_extrema(fftwf_complex * grid, double lower_limit, double upper_limit, double *grid_min, double *grid_max){ + double min_buf,max_buf; + //setting the extrema to the first cell to guarantee it's in the range + min_buf = *((float *) grid + HII_R_FFT_INDEX(0, 0, 0)); + max_buf = *((float *) grid + HII_R_FFT_INDEX(0, 0, 0)); + #pragma omp parallel num_threads(user_params_global->N_THREADS) { -#pragma omp for - for (i=0; iHII_DIM; i++){ - for (j=0; jHII_DIM; j++){ - for (k=0; kdensity[HII_R_INDEX(i,j,k)])*adjustment_factor; + int x,y,z; + float curr; + #pragma omp for reduction(max:max_buf) reduction(min:min_buf) + for (x = 0; x < user_params_global->HII_DIM; x++) { + for (y = 0; y < user_params_global->HII_DIM; y++) { + for (z = 0; z < HII_D_PARA; z++) { + // delta cannot be less than -1 + curr = *((float *) grid + HII_R_FFT_INDEX(x, y, z)); + *((float *) grid + HII_R_FFT_INDEX(x, y, z)) = fmaxf(fmin(curr,upper_limit),lower_limit); + + if (curr < min_buf) + min_buf = curr; + if (curr > max_buf) + max_buf = curr; } } } } + *grid_min = min_buf; + *grid_max = max_buf; +} - LOG_SUPER_DEBUG("density field calculated"); - // keep the unfiltered density field in an array, to save it for later - memcpy(deltax_unfiltered_original, deltax_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - - i=0; +//TODO: maybe put the grid clipping outside this function +void setup_integration_tables(struct FilteredGrids *fg_struct, struct IonBoxConstants *consts, struct RadiusSpec rspec, bool need_prev){ + double min_density, max_density, prev_min_density, prev_max_density; + double log10Mturn_min, log10Mturn_max, log10Mturn_min_MINI, log10Mturn_max_MINI; + + //TODO: instead of putting a random upper limit, put a proper flag for switching of one/both sides of the clipping + clip_and_get_extrema(fg_struct->deltax_filtered,-1,1e6,&min_density,&max_density); + min_density -= 0.001; + max_density += 0.001; + + if (flag_options_global->USE_MASS_DEPENDENT_ZETA){ + if (flag_options_global->USE_MINI_HALOS){ + // do the same for prev + clip_and_get_extrema(fg_struct->prev_deltax_filtered,-1,1e6,&prev_min_density,&prev_max_density); + clip_and_get_extrema(fg_struct->log10_Mturnover_filtered,0.,LOG10_MTURN_MAX,&log10Mturn_min,&log10Mturn_max); + clip_and_get_extrema(fg_struct->log10_Mturnover_MINI_filtered,0.,LOG10_MTURN_MAX,&log10Mturn_min_MINI,&log10Mturn_max_MINI); + } - // Newer setup to be performed in parallel - // TODO: examine this RNG model, seed properly with seed_rng_threads() from ICs - int thread_num; - for(thread_num = 0; thread_num < user_params->N_THREADS; thread_num++){ - // Original defaults with gsl_rng_mt19937 and SEED = 0, thus start with this and iterate for all other threads by their thread number - r[thread_num] = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r[thread_num], thread_num); + LOG_ULTRA_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)) + initialise_GL(consts->lnMmin, rspec.ln_M_max_R); + if(user_params_global->USE_INTERPOLATION_TABLES){ + //Buffers to avoid both zero bin widths and max cell segfault in 2D interptables + prev_min_density -= 0.001; + prev_max_density += 0.001; + log10Mturn_min = log10Mturn_min * 0.99; + log10Mturn_max = log10Mturn_max * 1.01; + log10Mturn_min_MINI = log10Mturn_min_MINI * 0.99; + log10Mturn_max_MINI = log10Mturn_max_MINI * 1.01; + + //current redshift tables (automatically handles minihalo case) + initialise_Nion_Conditional_spline(consts->redshift,consts->mturn_a_nofb,min_density,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, + consts->alpha_esc,consts->fstar_10, + consts->fesc_10,consts->Mlim_Fstar,consts->Mlim_Fesc,consts->fstar_7, + consts->fesc_7,consts->Mlim_Fstar_mini,consts->Mlim_Fesc_mini, + user_params_global->INTEGRATION_METHOD_ATOMIC, user_params_global->INTEGRATION_METHOD_MINI, + flag_options_global->USE_MINI_HALOS,false); + + //previous redshift tables if needed + if(need_prev && flag_options_global->USE_MINI_HALOS){ + 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, + consts->alpha_esc,consts->fstar_10, + consts->fesc_10,consts->Mlim_Fstar,consts->Mlim_Fesc,consts->fstar_7, + consts->fstar_7,consts->Mlim_Fstar_mini,consts->Mlim_Fesc_mini, + user_params_global->INTEGRATION_METHOD_ATOMIC, user_params_global->INTEGRATION_METHOD_MINI, + flag_options_global->USE_MINI_HALOS,true); + } + } } - - pixel_mass = RtoM(L_FACTOR*user_params->BOX_LEN/(float)(user_params->HII_DIM)); - cell_length_factor = L_FACTOR; - - if(flag_options->USE_HALO_FIELD && (global_params.FIND_BUBBLE_ALGORITHM == 2) && ((user_params->BOX_LEN/(float)(user_params->HII_DIM) < 1))) { - cell_length_factor = 1.; + else { + //This was previously one table for all R, which can be done with the EPS mass function (and some others) + //TODO: I don't expect this to be a bottleneck, but we can look into re-making the 2/3D ERFC tables if needed + if(user_params_global->USE_INTERPOLATION_TABLES) + initialise_FgtrM_delta_table(min_density, max_density, consts->redshift, consts->growth_factor, consts->sigma_minmass, rspec.sigma_maxmass); } +} - if (prev_redshift < 1){ -LOG_DEBUG("first redshift, do some initialization"); - previous_ionize_box->z_re_box = (float *) calloc(HII_TOT_NUM_PIXELS, sizeof(float)); -#pragma omp parallel shared(previous_ionize_box) private(i,j,k) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (i=0; iHII_DIM; i++){ - for (j=0; jHII_DIM; j++){ - for (k=0; kz_re_box[HII_R_INDEX(i, j, k)] = -1.0; +//TODO: We should speed test different configurations, separating grids, parallel sections etc. +// See the note above copy_filter_transform() for the general idea +// If we separate by grid we can reuse the clipping function above +void calculate_fcoll_grid(IonizedBox *box, IonizedBox *previous_ionize_box, struct FilteredGrids *fg_struct, struct IonBoxConstants *consts, + struct RadiusSpec *rspec){ + double f_coll_total,f_coll_MINI_total; + //TODO: make proper error tracking through the parallel region + bool error_flag; + + int fc_r_idx; + fc_r_idx = flag_options_global->USE_MINI_HALOS ? rspec->R_index : 0; + #pragma omp parallel num_threads(user_params_global->N_THREADS) + { + int x,y,z; + double curr_dens; + double Splined_Fcoll,Splined_Fcoll_MINI; + double log10_Mturnover,log10_Mturnover_MINI; + double prev_dens,prev_Splined_Fcoll,prev_Splined_Fcoll_MINI; + log10_Mturnover = log10(consts->mturn_a_nofb); //is only overwritten with minihalos + #pragma omp for reduction(+:f_coll_total,f_coll_MINI_total) + for (x = 0; x < user_params_global->HII_DIM; x++) { + for (y = 0; y < user_params_global->HII_DIM; y++) { + for (z = 0; z < HII_D_PARA; z++) { + //clip the filtered grids to physical values + // delta cannot be less than -1 + *((float *) fg_struct->deltax_filtered + HII_R_FFT_INDEX(x, y, z)) = fmaxf( + *((float *) fg_struct->deltax_filtered + HII_R_FFT_INDEX(x, y, z)), -1. + FRACT_FLOAT_ERR); + + // cannot be less than zero + if (consts->filter_recombinations) { + *((float *) fg_struct->N_rec_filtered + HII_R_FFT_INDEX(x, y, z)) = \ + fmaxf(*((float *) fg_struct->N_rec_filtered + HII_R_FFT_INDEX(x, y, z)), 0.0); } - } - } - } - if (flag_options->INHOMO_RECO) - previous_ionize_box->dNrec_box = (float *) calloc(HII_TOT_NUM_PIXELS, sizeof(float)); - } - //set the minimum source mass - M_MIN = minimum_source_mass(redshift,false,astro_params,flag_options); - lnMmin = log(M_MIN); - lnMmax = log(global_params.M_MAX_INTEGRAL); - LOG_SUPER_DEBUG("minimum source mass has been set: %f", M_MIN); - //Find the mass limits and average turnovers - if (flag_options->USE_MASS_DEPENDENT_ZETA) { - if (flag_options->USE_MINI_HALOS){ - ave_log10_Mturnover = 0.; - ave_log10_Mturnover_MINI = 0.; - - // this is the first z, and the previous_ionize_box are empty - if (prev_redshift < 1){ - previous_ionize_box->Gamma12_box = (float *) calloc(HII_TOT_NUM_PIXELS, sizeof(float)); - // really painful to get the length... - counter = 1; - R=fmax(global_params.R_BUBBLE_MIN, (cell_length_factor*user_params->BOX_LEN/(float)user_params->HII_DIM)); - while ((R - fmin(astro_params->R_BUBBLE_MAX, L_FACTOR*user_params->BOX_LEN)) <= FRACT_FLOAT_ERR ){ - if(R >= fmin(astro_params->R_BUBBLE_MAX, L_FACTOR*user_params->BOX_LEN)) { - stored_R = R/(global_params.DELTA_R_HII_FACTOR); + // x_e has to be between zero and unity + if (flag_options_global->USE_TS_FLUCT) { + *((float *) fg_struct->xe_filtered + HII_R_FFT_INDEX(x, y, z)) = \ + fmaxf(*((float *) fg_struct->xe_filtered + HII_R_FFT_INDEX(x, y, z)), 0.); + *((float *) fg_struct->xe_filtered + HII_R_FFT_INDEX(x, y, z)) = \ + fminf(*((float *) fg_struct->xe_filtered + HII_R_FFT_INDEX(x, y, z)), 0.999); } - R*= global_params.DELTA_R_HII_FACTOR; - counter += 1; - } - previous_ionize_box->Fcoll = (float *) calloc(HII_TOT_NUM_PIXELS*counter, sizeof(float)); - previous_ionize_box->Fcoll_MINI = (float *) calloc(HII_TOT_NUM_PIXELS*counter, sizeof(float)); - previous_ionize_box->mean_f_coll = 0.0; - previous_ionize_box->mean_f_coll_MINI = 0.0; - -#pragma omp parallel shared(prev_deltax_unfiltered) private(i,j,k) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (i=0; iHII_DIM; i++){ - for (j=0; jHII_DIM; j++){ - for (k=0; kUSE_HALO_FIELD) { + *((float *)fg_struct->stars_filtered + HII_R_FFT_INDEX(x,y,z)) = fmaxf( + *((float *)fg_struct->stars_filtered + HII_R_FFT_INDEX(x,y,z)) , 0.0); + *((float *)fg_struct->sfr_filtered + HII_R_FFT_INDEX(x,y,z)) = fmaxf( + *((float *)fg_struct->sfr_filtered + HII_R_FFT_INDEX(x,y,z)) , 0.0); + + //Ionising photon output + Splined_Fcoll = *((float *)fg_struct->stars_filtered + HII_R_FFT_INDEX(x,y,z)); + //Minihalos are taken care of already + Splined_Fcoll_MINI = 0.; + //The smoothing done with minihalos corrects for sudden changes in M_crit + //Nion_smoothed(z,Mcrit) = Nion(z,Mcrit) + (Nion(z_prev,Mcrit_prev) - Nion(z_prev,Mcrit)) + prev_Splined_Fcoll = 0.; + prev_Splined_Fcoll_MINI = 0.; + } + else { + curr_dens = *((float *) fg_struct->deltax_filtered + HII_R_FFT_INDEX(x, y, z)); + if (flag_options_global->USE_MASS_DEPENDENT_ZETA){ + if (flag_options_global->USE_MINI_HALOS){ + log10_Mturnover = *((float *)fg_struct->log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z)); + log10_Mturnover_MINI = *((float *)fg_struct->log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)); + + Splined_Fcoll_MINI = EvaluateNion_Conditional_MINI(curr_dens,log10_Mturnover_MINI,consts->growth_factor,consts->M_min, + rspec->M_max_R,rspec->M_max_R,rspec->sigma_maxmass,consts->mturn_a_nofb, + consts->Mlim_Fstar_mini,consts->Mlim_Fesc,false); + + + if (previous_ionize_box->mean_f_coll_MINI * consts->ion_eff_factor_mini_gl + + previous_ionize_box->mean_f_coll * consts->ion_eff_factor_gl > 1e-4){ + prev_dens = *((float *)fg_struct->prev_deltax_filtered + HII_R_FFT_INDEX(x,y,z)); + prev_Splined_Fcoll = EvaluateNion_Conditional(prev_dens,log10_Mturnover,consts->prev_growth_factor, + consts->M_min,rspec->M_max_R,rspec->M_max_R, + rspec->sigma_maxmass,consts->Mlim_Fstar,consts->Mlim_Fesc,true); + prev_Splined_Fcoll_MINI = EvaluateNion_Conditional_MINI(prev_dens,log10_Mturnover_MINI,consts->prev_growth_factor,consts->M_min, + rspec->M_max_R,rspec->M_max_R,rspec->sigma_maxmass,consts->mturn_a_nofb, + consts->Mlim_Fstar_mini,consts->Mlim_Fesc_mini,true); + } + else{ + prev_Splined_Fcoll = 0.; + prev_Splined_Fcoll_MINI = 0.; + } } + Splined_Fcoll = EvaluateNion_Conditional(curr_dens,log10_Mturnover,consts->growth_factor, + consts->M_min,rspec->M_max_R,rspec->M_max_R, + rspec->sigma_maxmass,consts->Mlim_Fstar,consts->Mlim_Fesc,false); + } + else{ + Splined_Fcoll = EvaluateFcoll_delta(curr_dens,consts->growth_factor,consts->sigma_minmass,rspec->sigma_maxmass); } } - } - } - else{ -#pragma omp parallel shared(prev_deltax_unfiltered,previous_perturbed_field) private(i,j,k) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (i=0; iHII_DIM; i++){ - for (j=0; jHII_DIM; j++){ - for (k=0; kdensity[HII_R_INDEX(i,j,k)]; - } + // save the value of the collasped fraction into the Fcoll array + //TODO: each of these grids are clipped before filtering, after filtering, + // after the Fcoll integration, and after trapezoidal integration here + // we should figure which of these clips are necessary/useful + if (flag_options_global->USE_MINI_HALOS && !flag_options_global->USE_HALO_FIELD){ + if (Splined_Fcoll > 1.) Splined_Fcoll = 1.; + if (Splined_Fcoll < 0.) Splined_Fcoll = 1e-40; + if (prev_Splined_Fcoll > 1.) prev_Splined_Fcoll = 1.; + if (prev_Splined_Fcoll < 0.) prev_Splined_Fcoll = 1e-40; + box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = \ + previous_ionize_box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] + Splined_Fcoll - prev_Splined_Fcoll; + + if (box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] > 1.) box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1.; + //if (box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] <0.) box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1e-40; + //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, + *((float *)fg_struct->log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z))); + Throw(InfinityorNaNError); + } + + if (Splined_Fcoll_MINI > 1.) Splined_Fcoll_MINI = 1.; + if (Splined_Fcoll_MINI < 0.) Splined_Fcoll_MINI = 1e-40; + if (prev_Splined_Fcoll_MINI > 1.) prev_Splined_Fcoll_MINI = 1.; + if (prev_Splined_Fcoll_MINI < 0.) prev_Splined_Fcoll_MINI = 1e-40; + box->Fcoll_MINI[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)] + Splined_Fcoll_MINI - prev_Splined_Fcoll_MINI; + + if (box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] >1.) box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1.; + //if (box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] <0.) box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1e-40; + //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)]); + Throw(InfinityorNaNError); } } + else{ + box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = Splined_Fcoll; + f_coll_total += Splined_Fcoll; + } } } + } + } // end loop through Fcoll box + rspec->f_coll_grid_mean = f_coll_total / HII_TOT_NUM_PIXELS; + rspec->f_coll_grid_mean_MINI = f_coll_MINI_total / HII_TOT_NUM_PIXELS; +} -LOG_SUPER_DEBUG("previous density field calculated"); +int setup_radii(struct RadiusSpec **rspec_array, struct IonBoxConstants *consts){ + double maximum_radius = fmin(astro_params_global->R_BUBBLE_MAX, + L_FACTOR*user_params_global->BOX_LEN); - // fields added for minihalos - Mcrit_atom = atomic_cooling_threshold(redshift); - log10_Mcrit_atom = log10(Mcrit_atom); - log10_Mcrit_mol = log10(lyman_werner_threshold(redshift, 0.,0., astro_params)); - log10_Mturnover_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - log10_Mturnover_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - log10_Mturnover_MINI_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - log10_Mturnover_MINI_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + double cell_length_factor = L_FACTOR; + //TODO: figure out why this is used in such a specific case + if(flag_options_global->USE_HALO_FIELD && (global_params.FIND_BUBBLE_ALGORITHM == 2) && (consts->pixel_length < 1)) + cell_length_factor = 1.; - if (!log10_Mturnover_unfiltered || !log10_Mturnover_filtered || !log10_Mturnover_MINI_unfiltered || !log10_Mturnover_MINI_filtered){// || !Mcrit_RE_grid || !Mcrit_LW_grid) - LOG_ERROR("Error allocating memory for Mturnover or Mturnover_MINI boxes"); - Throw(MemoryAllocError); - } -LOG_SUPER_DEBUG("Calculating and outputting Mcrit boxes for atomic and molecular halos..."); - -#pragma omp parallel shared(redshift,previous_ionize_box,spin_temp,Mcrit_atom,log10_Mturnover_unfiltered,log10_Mturnover_MINI_unfiltered)\ - private(x,y,z,Mcrit_RE,Mcrit_LW,Mturnover,Mturnover_MINI,log10_Mturnover,log10_Mturnover_MINI,curr_vcb) num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(+:ave_log10_Mturnover,ave_log10_Mturnover_MINI) - for (x=0; xHII_DIM; x++){ - for (y=0; yHII_DIM; y++){ - for (z=0; zGamma12_box[HII_R_INDEX(x, y, z)], previous_ionize_box->z_re_box[HII_R_INDEX(x, y, z)]); - if (flag_options->FIX_VCB_AVG){ //with this flag we ignore reading vcb box - curr_vcb = global_params.VAVG; + double minimum_radius = fmax(global_params.R_BUBBLE_MIN,cell_length_factor*consts->pixel_length); + + //minimum number such that min_R*delta^N > max_R + int n_radii = (int)(log(maximum_radius/minimum_radius)/log(global_params.DELTA_R_HII_FACTOR) + 1); + *rspec_array = malloc(sizeof(**rspec_array)*n_radii); + + //We want the following behaviour from our radius Values: + // The smallest radius is the cell size or global min + // The largest radius is the box size of global max + // Each step is set by multiplying by the same factor + //This is not possible for most sets of these three parameters + // so we let the first step (largest -> second largest) be different, + // finding the other radii by stepping *up* from the minimum + int i; + for(i=0;i maximum_radius - FRACT_FLOAT_ERR){ + (*rspec_array)[i].R = maximum_radius; + n_radii = i+1; //also ends the loop after this iteration + } + (*rspec_array)[i].M_max_R = RtoM((*rspec_array)[i].R); + (*rspec_array)[i].ln_M_max_R = log((*rspec_array)[i].M_max_R); + (*rspec_array)[i].sigma_maxmass = sigma_z0((*rspec_array)[i].M_max_R); + } + LOG_DEBUG("set max radius: %f", (*rspec_array)[n_radii-1].R); + return n_radii; +} + +void find_ionised_regions(IonizedBox *box, IonizedBox *previous_ionize_box, PerturbedField *perturbed_field, TsBox *spin_temp, struct RadiusSpec rspec, + struct IonBoxConstants *consts, struct FilteredGrids *fg_struct, double f_limit_acg, double f_limit_mcg){ + double mean_fix_term_acg = 1.; + double mean_fix_term_mcg = 1.; + int fc_r_idx; + fc_r_idx = flag_options_global->USE_MINI_HALOS ? rspec.R_index : 0; + + LOG_SUPER_DEBUG("global mean fcoll (mini) %.3e (%.3e) box mean fcoll %.3e (%.3e) ratio %.3e (%.3e)", + box->mean_f_coll,box->mean_f_coll_MINI, + rspec.f_coll_grid_mean,rspec.f_coll_grid_mean_MINI, + box->mean_f_coll/rspec.f_coll_grid_mean,box->mean_f_coll/rspec.f_coll_grid_mean); + if(consts->fix_mean){ + mean_fix_term_acg = box->mean_f_coll/rspec.f_coll_grid_mean; + if(flag_options_global->USE_MINI_HALOS){ + mean_fix_term_mcg = box->mean_f_coll_MINI/rspec.f_coll_grid_mean_MINI; + } + } + else{ + //if we don't fix the mean, make the mean_f_coll in the output reflect the box + //since currently it is the global expected mean from the Unconditional MF + box->mean_f_coll = rspec.f_coll_grid_mean; + box->mean_f_coll_MINI = rspec.f_coll_grid_mean_MINI; + } + + #pragma omp parallel num_threads(user_params_global->N_THREADS) + { + int x,y,z; + double curr_dens; + double curr_fcoll; + double curr_fcoll_mini; + double rec,xHII_from_xrays; + double res_xH; + #pragma omp for + for(x = 0; x < user_params_global->HII_DIM; x++) { + for(y = 0; y < user_params_global->HII_DIM; y++) { + for(z = 0; z < HII_D_PARA; z++) { + //Use unfiltered density for CELL_RECOMB case, since the "Fcoll" represents photons + // reaching the central cell rather than photons in the entire sphere + //TODO: test using the filtered density here with CELL_RECOMB, since it's what Davies+22 does, + // It's not clear how the MFP filter should count recombinations/local density. + //If we don't filter on the last step, might as well use the original field to prevent aliasing? + if(rspec.R_index==0) + curr_dens = perturbed_field->density[HII_R_INDEX(x,y,z)]*consts->photoncons_adjustment_factor; + else + curr_dens = *((float *)fg_struct->deltax_filtered + HII_R_FFT_INDEX(x,y,z)); + + curr_fcoll = box->Fcoll[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; + curr_fcoll = mean_fix_term_acg * curr_fcoll; + + //Since the halo boxes give ionising photon output, this term accounts for the local density of absorbers + // We have separated the source/absorber filtering in the halo model so this is necessary + if(flag_options_global->USE_HALO_FIELD) + curr_fcoll *= 1 / (RHOcrit*cosmo_params_global->OMb*(1+curr_dens)); + + //MINIHALOS are already included in the halo model + curr_fcoll_mini = 0.; + if (flag_options_global->USE_MINI_HALOS && !flag_options_global->USE_HALO_FIELD){ + curr_fcoll_mini = box->Fcoll_MINI[fc_r_idx * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; + curr_fcoll_mini = mean_fix_term_mcg * curr_fcoll_mini; + } + + if (flag_options_global->USE_MASS_DEPENDENT_ZETA) { + if (curr_fcoll < f_limit_acg) curr_fcoll = f_limit_acg; + if (flag_options_global->USE_MINI_HALOS){ + if (curr_fcoll_mini < f_limit_mcg) curr_fcoll_mini = f_limit_mcg; + } + } + + if (flag_options_global->INHOMO_RECO) { + if(flag_options_global->CELL_RECOMB) + rec = previous_ionize_box->dNrec_box[HII_R_INDEX(x,y,z)]; + else + rec = (*((float *) fg_struct->N_rec_filtered + HII_R_FFT_INDEX(x, y, z))); // number of recombinations per mean baryon + + rec /= (1. + curr_dens); // number of recombinations per baryon inside cell/filter + } + else { + rec = 0.; + } + + // adjust the denominator of the collapse fraction for the residual electron fraction in the neutral medium + if (flag_options_global->USE_TS_FLUCT){ + xHII_from_xrays = *((float *)fg_struct->xe_filtered + HII_R_FFT_INDEX(x,y,z)); + } else { + xHII_from_xrays = 0.; + } + +#if LOG_LEVEL >= SUPER_DEBUG_LEVEL + if(x+y+z == 0){ + LOG_SUPER_DEBUG("Cell 0: R=%.1f | d %.4e | fcoll %.4e (%.4e) | rec %.4e | X %.4e", + rspec.R,curr_dens,curr_fcoll,curr_fcoll_mini,rec,xHII_from_xrays); + } +#endif + + // check if fully ionized! + if ( (curr_fcoll * consts->ion_eff_factor + curr_fcoll_mini * consts->ion_eff_factor_mini > (1. - xHII_from_xrays)*(1.0+rec)) ){ //IONIZED!! + // if this is the first crossing of the ionization barrier for this cell (largest R), record the gamma + // this assumes photon-starved growth of HII regions... breaks down post EoR + if (flag_options_global->INHOMO_RECO && (box->xH_box[HII_R_INDEX(x,y,z)] > FRACT_FLOAT_ERR)){ + if(flag_options_global->USE_HALO_FIELD){ + //since ION_EFF_FACTOR==1 here, gamma_prefactor is the same for ACG and MCG + box->Gamma12_box[HII_R_INDEX(x,y,z)] = rspec.R * consts->gamma_prefactor / (1+curr_dens) \ + * (*((float *)fg_struct->sfr_filtered + HII_R_FFT_INDEX(x,y,z))); } else{ - if(user_params->USE_RELATIVE_VELOCITIES){ - curr_vcb = ini_boxes->lowres_vcb[HII_R_INDEX(x,y,z)]; - } - else{ //set vcb to a constant, either zero or vavg. - curr_vcb = 0.0; - } + box->Gamma12_box[HII_R_INDEX(x,y,z)] = rspec.R * (consts->gamma_prefactor * curr_fcoll + consts->gamma_prefactor_mini * curr_fcoll_mini); } + box->MFP_box[HII_R_INDEX(x,y,z)] = rspec.R; + } - Mcrit_LW = lyman_werner_threshold(redshift, spin_temp->J_21_LW_box[HII_R_INDEX(x, y, z)], curr_vcb, astro_params); - if(Mcrit_LW != Mcrit_LW || Mcrit_LW == 0){ - LOG_ERROR("Mcrit error %d %d %d: M %.2e z %.2f J %.2e v %.2e",x,y,z,Mcrit_LW,redshift,spin_temp->J_21_LW_box[HII_R_INDEX(x, y, z)],curr_vcb); - Throw(ValueError); - } + // keep track of the first time this cell is ionized (earliest time) + if (previous_ionize_box->z_re_box[HII_R_INDEX(x,y,z)] < 0){ + box->z_re_box[HII_R_INDEX(x,y,z)] = consts->redshift; //TODO: stored_redshift??? + } else{ + box->z_re_box[HII_R_INDEX(x,y,z)] = previous_ionize_box->z_re_box[HII_R_INDEX(x,y,z)]; + } - //JBM: this only accounts for effect 3 (largest on minihaloes). Effects 1 and 2 affect both minihaloes (MCGs) and regular ACGs, but they're smaller ~10%. See Sec 2 of Muñoz+21 (2110.13919) + // FLAG CELL(S) AS IONIZED + if (global_params.FIND_BUBBLE_ALGORITHM == 2) // center method + box->xH_box[HII_R_INDEX(x,y,z)] = 0; + else if (global_params.FIND_BUBBLE_ALGORITHM == 1) // sphere method + update_in_sphere(box->xH_box, user_params_global->HII_DIM, HII_D_PARA, rspec.R/(user_params_global->BOX_LEN), \ + x/(user_params_global->HII_DIM+0.0), y/(user_params_global->HII_DIM+0.0), z/(HII_D_PARA+0.0)); + } // end ionized + // If not fully ionized, then assign partial ionizations + else if (rspec.R_index==0 && (box->xH_box[HII_R_INDEX(x, y, z)] > TINY)) { + //NOTE: Previously there was an RNG model here which multiplied Fcoll by a sampled + // Poisson/ term if 1/5 < M_coll / M_min < 5. This only ever affected the + // old parametrisation due to the M_min term. + + res_xH = 1. - curr_fcoll * consts->ion_eff_factor - curr_fcoll_mini * consts->ion_eff_factor_mini; + // put the partial ionization here because we need to exclude xHII_from_xrays... + if (flag_options_global->USE_TS_FLUCT){ + box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = ComputePartiallyIoinizedTemperature(spin_temp->Tk_box[HII_R_INDEX(x,y,z)], res_xH); + } + else{ + box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = ComputePartiallyIoinizedTemperature( + consts->TK_nofluct*(1 + consts->adia_TK_term*perturbed_field->density[HII_R_INDEX(x,y,z)]), res_xH); + } + res_xH -= xHII_from_xrays; + // and make sure fraction doesn't blow up for underdense pixels + if (res_xH < 0) + res_xH = 0; + else if (res_xH > 1) + res_xH = 1; - //*((float *)Mcrit_RE_grid + HII_R_FFT_INDEX(x,y,z)) = Mcrit_RE; - //*((float *)Mcrit_LW_grid + HII_R_FFT_INDEX(x,y,z)) = Mcrit_LW; - Mturnover = Mcrit_RE > Mcrit_atom ? Mcrit_RE : Mcrit_atom; - Mturnover = fmax(Mturnover,astro_params->M_TURN); - Mturnover_MINI = Mcrit_RE > Mcrit_LW ? Mcrit_RE : Mcrit_LW; - Mturnover_MINI = fmax(Mturnover_MINI,astro_params->M_TURN); - log10_Mturnover = log10(Mturnover); - log10_Mturnover_MINI = log10(Mturnover_MINI); + box->xH_box[HII_R_INDEX(x, y, z)] = res_xH; - *((float *)log10_Mturnover_unfiltered + HII_R_FFT_INDEX(x,y,z)) = log10_Mturnover; - *((float *)log10_Mturnover_MINI_unfiltered + HII_R_FFT_INDEX(x,y,z)) = log10_Mturnover_MINI; + } // end partial ionizations at last filtering step + } // k + } // j + } // i + } +} - ave_log10_Mturnover += log10_Mturnover; - ave_log10_Mturnover_MINI += log10_Mturnover_MINI; +void set_ionized_temperatures(IonizedBox *box, PerturbedField *perturbed_field, TsBox *spin_temp, struct IonBoxConstants *consts){ + int x,y,z; + #pragma omp parallel private(x,y,z) num_threads(user_params_global->N_THREADS) + { + float thistk; + #pragma omp for + for (x=0; xHII_DIM; x++){ + for (y=0; yHII_DIM; y++){ + for (z=0; zz_re_box[HII_R_INDEX(x,y,z)]>0) && (box->xH_box[HII_R_INDEX(x,y,z)] < TINY)){ + //TODO: do we want to use the photoncons redshift here or the original one? + box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = ComputeFullyIoinizedTemperature(box->z_re_box[HII_R_INDEX(x,y,z)], \ + consts->stored_redshift, perturbed_field->density[HII_R_INDEX(x,y,z)]); + // Below sometimes (very rare though) can happen when the density drops too fast and to below T_HI + if (flag_options_global->USE_TS_FLUCT){ + if (box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] < spin_temp->Tk_box[HII_R_INDEX(x,y,z)]) + box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = spin_temp->Tk_box[HII_R_INDEX(x,y,z)]; + } + else{ + thistk = consts->TK_nofluct*(1 + consts->adia_TK_term * \ + perturbed_field->density[HII_R_INDEX(x,y,z)]); + if (box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] < thistk) + box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = thistk; } } } } - - box->log10_Mturnover_ave = ave_log10_Mturnover/(double) HII_TOT_NUM_PIXELS; - box->log10_Mturnover_MINI_ave = ave_log10_Mturnover_MINI/(double) HII_TOT_NUM_PIXELS; - Mturnover = pow(10., box->log10_Mturnover_ave); - Mturnover_MINI = pow(10., box->log10_Mturnover_MINI_ave); - Mlim_Fstar_MINI = Mass_limit_bisection(M_MIN, 1e16, astro_params->ALPHA_STAR_MINI, astro_params->F_STAR7_MINI * pow(1e3,astro_params->ALPHA_STAR_MINI)); - Mlim_Fesc_MINI = Mass_limit_bisection(M_MIN, 1e16, alpha_esc_var, astro_params->F_ESC7_MINI * pow(1e3, alpha_esc_var)); - LOG_DEBUG("average log10 turnover masses are %.2f and %.2f for ACGs and MCGs", box->log10_Mturnover_ave, box->log10_Mturnover_MINI_ave); } - else{ - Mturnover = astro_params->M_TURN; - Mcrit_atom = Mturnover; //for table init - box->log10_Mturnover_ave = log10(Mturnover); - box->log10_Mturnover_MINI_ave = log10(Mturnover); - } - Mlim_Fstar = Mass_limit_bisection(M_MIN, 1e16, astro_params->ALPHA_STAR, astro_params->F_STAR10); - Mlim_Fesc = Mass_limit_bisection(M_MIN, 1e16, alpha_esc_var, norm_esc_var); } - if(user_params->USE_INTERPOLATION_TABLES) { - if(user_params->INTEGRATION_METHOD_ATOMIC == 2 || user_params->INTEGRATION_METHOD_MINI == 2){ - initialiseSigmaMInterpTable(fmin(MMIN_FAST,M_MIN),1e20); - } - else{ - initialiseSigmaMInterpTable(M_MIN,1e20); + for (x=0; xHII_DIM; x++){ + for (y=0; yHII_DIM; y++){ + for (z=0; ztemp_kinetic_all_gas[HII_R_INDEX(x,y,z)])==0){ + LOG_ERROR("Tk after fully ioinzation is either infinite or a Nan. Something has gone wrong "\ + "in the temperature calculation: z_re=%.4f, redshift=%.4f, curr_dens=%.4e", box->z_re_box[HII_R_INDEX(x,y,z)], consts->stored_redshift, + perturbed_field->density[HII_R_INDEX(x,y,z)]); + Throw(InfinityorNaNError); + } + } } - - } - - -LOG_SUPER_DEBUG("sigma table has been initialised"); - - // check for WDM - - if (global_params.P_CUTOFF && ( M_MIN < M_J_WDM())){ - LOG_WARNING("The default Jeans mass of %e Msun is smaller than the scale supressed by the effective pressure of WDM.", M_MIN); - M_MIN = M_J_WDM(); - LOG_WARNING("Setting a new effective Jeans mass from WDM pressure supression of %e Msun", M_MIN); } +} - // ARE WE USING A DISCRETE HALO FIELD (identified in the ICs with FindHaloes.c and evolved with PerturbHaloField.c) - if(flag_options->USE_HALO_FIELD) { - double nion_avg=0; - double wsfr_avg=0; - stars_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - stars_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); +void set_recombination_rates(IonizedBox *box, IonizedBox *previous_ionize_box, PerturbedField *perturbed_field, struct IonBoxConstants *consts){ + bool finite_error; + #pragma omp parallel num_threads(user_params_global->N_THREADS) + { + int x,y,z; + double curr_dens,dNrec; + double z_eff; + #pragma omp for + for(x = 0; x < user_params_global->HII_DIM; x++) { + for(y = 0; y < user_params_global->HII_DIM; y++) { + for(z = 0; z < HII_D_PARA; z++) { + // use the original density and redshift for the snapshot (not the adjusted redshift) + // Only want to use the adjusted redshift for the ionisation field + //NOTE: but the structure field wasn't adjusted, this seems wrong + // curr_dens = 1.0 + (perturbed_field->density[HII_R_INDEX(x, y, z)]) / adjustment_factor; + curr_dens = 1.0 + (perturbed_field->density[HII_R_INDEX(x, y, z)]); + z_eff = pow(curr_dens, 1.0 / 3.0); + z_eff *= (1 + consts->stored_redshift); + + dNrec = splined_recombination_rate(z_eff - 1., box->Gamma12_box[HII_R_INDEX(x, y, z)]) * + consts->fabs_dtdz * consts->dz * (1. - box->xH_box[HII_R_INDEX(x, y, z)]); + + if (isfinite(dNrec) == 0) { + finite_error = true; + } - sfr_unfiltered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - sfr_filtered = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); + box->dNrec_box[HII_R_INDEX(x, y, z)] = + previous_ionize_box->dNrec_box[HII_R_INDEX(x, y, z)] + dNrec; -#pragma omp parallel shared(stars_unfiltered,sfr_unfiltered,halos) private(i,j,k) num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(+:nion_avg,wsfr_avg) - for (i=0; iHII_DIM; i++){ - for (j=0; jHII_DIM; j++){ - for (k=0; kn_ion[HII_R_INDEX(i,j,k)]; - *((float *)sfr_unfiltered + HII_R_FFT_INDEX(i,j,k)) = halos->whalo_sfr[HII_R_INDEX(i,j,k)]; - nion_avg += halos->n_ion[HII_R_INDEX(i,j,k)]; - wsfr_avg += halos->whalo_sfr[HII_R_INDEX(i,j,k)]; +#if LOG_LEVEL >= SUPER_DEBUG_LEVEL + if(x+y+z == 0){ + LOG_SUPER_DEBUG("Cell 0: d %.4e | G12 %.4e | xH %.4e ==> dNrec %.4e Nrec (%.4e --> %.4e)", + curr_dens,box->Gamma12_box[HII_R_INDEX(x,y,z)], + box->xH_box[HII_R_INDEX(x,y,z)],dNrec, + previous_ionize_box->dNrec_box[HII_R_INDEX(x, y, z)], + box->dNrec_box[HII_R_INDEX(x, y, z)]); } +#endif } } } - LOG_DEBUG("HaloBox, avg n_ion %.3e avg w_sfr %.3e",nion_avg/HII_TOT_NUM_PIXELS,wsfr_avg/HII_TOT_NUM_PIXELS); - } // end of the USE_HALO_FIELD option + } - // lets check if we are going to bother with computing the inhmogeneous field at all... - global_xH = 0.0; + if(finite_error) { + LOG_ERROR("Recombinations have returned either an infinite or NaN value."); + Throw(InfinityorNaNError); + } +} - if(user_params->INTEGRATION_METHOD_ATOMIC == 1 || (flag_options->USE_MINI_HALOS && user_params->INTEGRATION_METHOD_MINI == 1)) - initialise_GL(lnMmin, lnMmax); - - // Determine the normalisation for the excursion set algorithm - if (flag_options->USE_MASS_DEPENDENT_ZETA) { - if (flag_options->USE_MINI_HALOS){ - if (previous_ionize_box->mean_f_coll * ION_EFF_FACTOR < 1e-4){ - box->mean_f_coll = Nion_General(redshift,lnMmin,lnMmax,Mturnover,astro_params->ALPHA_STAR,alpha_esc_var, - astro_params->F_STAR10,norm_esc_var,Mlim_Fstar,Mlim_Fesc); - } - else{ - box->mean_f_coll = previous_ionize_box->mean_f_coll + \ - Nion_General(redshift,lnMmin,lnMmax,Mturnover,astro_params->ALPHA_STAR,alpha_esc_var, - astro_params->F_STAR10,norm_esc_var,Mlim_Fstar,Mlim_Fesc) - \ - Nion_General(prev_redshift,lnMmin,lnMmax,Mturnover,astro_params->ALPHA_STAR,alpha_esc_var, - astro_params->F_STAR10,norm_esc_var,Mlim_Fstar,Mlim_Fesc); - } - if (previous_ionize_box->mean_f_coll_MINI * ION_EFF_FACTOR_MINI < 1e-4){ - box->mean_f_coll_MINI = Nion_General_MINI(redshift,lnMmin,lnMmax,Mturnover_MINI,Mcrit_atom, - astro_params->ALPHA_STAR_MINI,alpha_esc_var,astro_params->F_STAR7_MINI, - astro_params->F_ESC7_MINI,Mlim_Fstar_MINI,Mlim_Fesc_MINI); - } - else{ - box->mean_f_coll_MINI = previous_ionize_box->mean_f_coll_MINI + \ - Nion_General_MINI(redshift,lnMmin,lnMmax,Mturnover_MINI,Mcrit_atom,astro_params->ALPHA_STAR_MINI, - alpha_esc_var,astro_params->F_STAR7_MINI,astro_params->F_ESC7_MINI - ,Mlim_Fstar_MINI,Mlim_Fesc_MINI) - \ - Nion_General_MINI(prev_redshift,lnMmin,lnMmax,Mturnover_MINI,Mcrit_atom,astro_params->ALPHA_STAR_MINI, - alpha_esc_var,astro_params->F_STAR7_MINI,astro_params->F_ESC7_MINI, - Mlim_Fstar_MINI,Mlim_Fesc_MINI); - } - f_coll_min = Nion_General(global_params.Z_HEAT_MAX,lnMmin,lnMmax,Mturnover,astro_params->ALPHA_STAR, - alpha_esc_var,astro_params->F_STAR10,norm_esc_var,Mlim_Fstar,Mlim_Fesc); - f_coll_min_MINI = Nion_General_MINI(global_params.Z_HEAT_MAX,lnMmin,lnMmax,Mturnover_MINI,Mcrit_atom, - astro_params->ALPHA_STAR_MINI,alpha_esc_var,astro_params->F_STAR7_MINI, - astro_params->F_ESC7_MINI,Mlim_Fstar_MINI,Mlim_Fesc_MINI); - } - else{ - box->mean_f_coll = Nion_General(redshift,lnMmin,lnMmax,Mturnover,astro_params->ALPHA_STAR,alpha_esc_var, - astro_params->F_STAR10,norm_esc_var,Mlim_Fstar,Mlim_Fesc); - box->mean_f_coll_MINI = 0.; +int ComputeIonizedBox(float redshift, float prev_redshift, UserParams *user_params, CosmoParams *cosmo_params, + AstroParams *astro_params, FlagOptions *flag_options, + PerturbedField *perturbed_field, + PerturbedField *previous_perturbed_field, + IonizedBox *previous_ionize_box, + TsBox *spin_temp, + HaloBox *halos, + InitialConditions *ini_boxes, + IonizedBox *box) { - f_coll_min = Nion_General(global_params.Z_HEAT_MAX,lnMmin,lnMmax,Mturnover,astro_params->ALPHA_STAR,alpha_esc_var, - astro_params->F_STAR10,norm_esc_var,Mlim_Fstar,Mlim_Fesc); - } - } - else { - LOG_DEBUG("Setting mean collapse fraction"); - box->mean_f_coll = Fcoll_General(redshift, lnMmin, lnMmax); - } + int status; - if(isfinite(box->mean_f_coll)==0) { - LOG_ERROR("Mean collapse fraction is either infinite or NaN!"); -// Throw(ParameterError); - Throw(InfinityorNaNError); - } -LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll: %e", box->mean_f_coll); + Try{ // This Try brackets the whole function, so we don't indent. + LOG_DEBUG("input values:"); + LOG_DEBUG("redshift=%f, prev_redshift=%f", redshift, prev_redshift); + #if LOG_LEVEL >= SUPER_DEBUG_LEVEL + writeUserParams(user_params); + writeCosmoParams(cosmo_params); + writeAstroParams(flag_options, astro_params); + writeFlagOptions(flag_options); + #endif - if (flag_options->USE_MINI_HALOS){ - if(isfinite(box->mean_f_coll_MINI)==0) { - LOG_ERROR("Mean collapse fraction of MINI is either infinite or NaN!"); -// Throw(ParameterError); - Throw(InfinityorNaNError); - } -LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll_MINI: %e", box->mean_f_coll_MINI); - } + // Makes the parameter structs visible to a variety of functions/macros + // Do each time to avoid Python garbage collection issues + Broadcast_struct_global_all(user_params,cosmo_params,astro_params,flag_options); + omp_set_num_threads(user_params->N_THREADS); - if (box->mean_f_coll * ION_EFF_FACTOR + box->mean_f_coll_MINI * ION_EFF_FACTOR_MINI< global_params.HII_ROUND_ERR){ // way too small to ionize anything... - // printf( "The mean collapse fraction is %e, which is much smaller than the effective critical collapse fraction of %e\n I will just declare everything to be neutral\n", mean_f_coll, f_coll_crit); + unsigned long long ct; - // find the neutral fraction - if(flag_options->USE_TS_FLUCT) { -#pragma omp parallel shared(box,spin_temp) private(ct) num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(+:global_xH) - for (ct=0; ctxH_box[ct] = 1.-spin_temp->x_e_box[ct]; // convert from x_e to xH - global_xH += box->xH_box[ct]; - box->temp_kinetic_all_gas[ct] = spin_temp->Tk_box[ct]; - } - } - global_xH /= (double)HII_TOT_NUM_PIXELS; - } - else { - global_xH = 1. - xion_RECFAST(redshift, 0); - -#pragma omp parallel shared(box,global_xH,TK,perturbed_field,cT_ad) private(ct) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (ct=0; ctxH_box[ct] = global_xH; - box->temp_kinetic_all_gas[ct] = TK * (1.0 + cT_ad * perturbed_field->density[ct]); // Is perturbed_field defined already here? we need it for cT. I'm also assuming we don't need to multiply by other z here. - } - } + double global_xH; + + //TODO: check if this is used in this file with TS fluctuations + init_heat(); + init_ps(); + + struct IonBoxConstants ionbox_constants; + set_ionbox_constants(redshift,prev_redshift,cosmo_params,astro_params,flag_options,&ionbox_constants); + + //boxes which aren't guaranteed to have every element assigned to need to be initialised + if(flag_options->INHOMO_RECO) { + if(INIT_RECOMBINATIONS) { + init_MHR(); + INIT_RECOMBINATIONS=0; } } - else { - // Take the ionisation fraction from the X-ray ionisations from Ts.c (only if the calculate spin temperature flag is set) - if (flag_options->USE_TS_FLUCT) { -#pragma omp parallel shared(xe_unfiltered, spin_temp) private(i, j, k) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (i = 0; i < user_params->HII_DIM; i++) { - for (j = 0; j < user_params->HII_DIM; j++) { - for (k = 0; k < HII_D_PARA; k++) { - *((float *) xe_unfiltered + HII_R_FFT_INDEX(i, j, k)) = spin_temp->x_e_box[HII_R_INDEX(i, j, k)]; - } - } - } - } - } - LOG_SUPER_DEBUG("calculated ionization fraction"); - - if (recomb_filter_flag) { -#pragma omp parallel shared(N_rec_unfiltered, previous_ionize_box) private(i, j, k) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (i = 0; i < user_params->HII_DIM; i++) { - for (j = 0; j < user_params->HII_DIM; j++) { - for (k = 0; k < HII_D_PARA; k++) { - *((float *) N_rec_unfiltered + - HII_R_FFT_INDEX(i, j, k)) = previous_ionize_box->dNrec_box[HII_R_INDEX(i, j, k)]; - } - } - } - } + #pragma omp parallel shared(box) private(ct) num_threads(user_params->N_THREADS) + { + #pragma omp for + for (ct=0; ctz_re_box[ct] = -1.0; } - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, deltax_unfiltered); - LOG_SUPER_DEBUG("FFTs performed"); + } - if(flag_options->USE_MINI_HALOS){ - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, prev_deltax_unfiltered); - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, log10_Mturnover_MINI_unfiltered); - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, log10_Mturnover_unfiltered); - LOG_SUPER_DEBUG("MINI HALO ffts performed"); + LOG_SUPER_DEBUG("z_re_box init: "); + debugSummarizeBox(box->z_re_box, user_params->HII_DIM, user_params->NON_CUBIC_FACTOR, " "); + + // Modify the current sampled redshift to a redshift which matches the expected filling factor given our astrophysical parameterisation. + // This is the photon non-conservation correction + float absolute_delta_z = 0.; + float redshift_pc, stored_redshift_pc; + if(flag_options->PHOTON_CONS_TYPE == 1) { + redshift_pc = redshift; + adjust_redshifts_for_photoncons(astro_params,flag_options,&redshift_pc,&stored_redshift_pc,&absolute_delta_z); + ionbox_constants.redshift = redshift_pc; + ionbox_constants.stored_redshift = stored_redshift_pc; + ionbox_constants.photoncons_adjustment_factor = dicke(redshift_pc)/dicke(stored_redshift_pc); + //TODO: if we want this to work with minihalos we need to change prev_redshift too + LOG_DEBUG("PhotonCons data:"); + LOG_DEBUG("original redshift=%f, updated redshift=%f delta-z = %f", stored_redshift_pc, redshift_pc, absolute_delta_z); + if(isfinite(redshift_pc)==0 || isfinite(absolute_delta_z)==0) { + LOG_ERROR("Updated photon non-conservation redshift is either infinite or NaN!"); + LOG_ERROR("This can sometimes occur when reionisation stalls (i.e. extremely low"\ + "F_ESC or F_STAR or not enough sources)"); + Throw(PhotonConsError); } + } + ///////////////////////////////// BEGIN INITIALIZATION ////////////////////////////////// + + struct RadiusSpec *radii_spec; + int n_radii; + n_radii = setup_radii(&radii_spec,&ionbox_constants); + + //CONSTRUCT GRIDS OUTSIDE R LOOP HERE + //if we don't have a previous ionised box, make a fake one here + if (prev_redshift < 1) + setup_first_z_prevbox(previous_ionize_box,previous_perturbed_field,n_radii); + struct FilteredGrids *grid_struct; + allocate_fftw_grids(&grid_struct); + + //Find the mass limits and average turnovers + double Mturnover_global_avg, Mturnover_global_avg_MINI; + if (flag_options->USE_MASS_DEPENDENT_ZETA){ if (flag_options->USE_HALO_FIELD){ - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, stars_unfiltered); - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, sfr_unfiltered); - LOG_SUPER_DEBUG("HALO_FIELD ffts performed"); + //Here these are only used for the global calculations + box->log10_Mturnover_ave = halos->log10_Mcrit_ACG_ave; + box->log10_Mturnover_MINI_ave = halos->log10_Mcrit_MCG_ave; + Mturnover_global_avg = pow(10., halos->log10_Mcrit_ACG_ave); + Mturnover_global_avg_MINI = pow(10., halos->log10_Mcrit_MCG_ave); } - - if(flag_options->USE_TS_FLUCT) { - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, xe_unfiltered); - LOG_SUPER_DEBUG("Ts ffts performed"); + else if (flag_options->USE_MINI_HALOS){ + LOG_SUPER_DEBUG("Calculating and outputting Mcrit boxes for atomic and molecular halos..."); + calculate_mcrit_boxes(previous_ionize_box, + spin_temp, + ini_boxes, + &ionbox_constants, + grid_struct->log10_Mturnover_unfiltered, + grid_struct->log10_Mturnover_MINI_unfiltered, + &(box->log10_Mturnover_ave), + &(box->log10_Mturnover_MINI_ave)); + + Mturnover_global_avg = pow(10., box->log10_Mturnover_ave); + Mturnover_global_avg_MINI = pow(10., box->log10_Mturnover_MINI_ave); + LOG_DEBUG("average log10 turnover masses are %.2f and %.2f for ACGs and MCGs", box->log10_Mturnover_ave, box->log10_Mturnover_MINI_ave); } - - if (recomb_filter_flag) { - dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, N_rec_unfiltered); + else{ + Mturnover_global_avg = astro_params->M_TURN; + box->log10_Mturnover_ave = log10(Mturnover_global_avg); + box->log10_Mturnover_MINI_ave = log10(Mturnover_global_avg); } + } + // lets check if we are going to bother with computing the inhmogeneous field at all... + global_xH = 0.0; - // remember to add the factor of VOLUME/TOT_NUM_PIXELS when converting from - // real space to k-space - // Note: we will leave off factor of VOLUME, in anticipation of the inverse FFT below -#pragma omp parallel shared(deltax_unfiltered,xe_unfiltered,N_rec_unfiltered,prev_deltax_unfiltered,\ - log10_Mturnover_unfiltered,log10_Mturnover_MINI_unfiltered,stars_unfiltered) \ - private(ct) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (ct=0; ctUSE_TS_FLUCT) { xe_unfiltered[ct] /= (double)HII_TOT_NUM_PIXELS; } - if (recomb_filter_flag){ N_rec_unfiltered[ct] /= (double)HII_TOT_NUM_PIXELS; } - if(flag_options->USE_HALO_FIELD) { - stars_unfiltered[ct] /= (double)HII_TOT_NUM_PIXELS; - sfr_unfiltered[ct] /= (double)HII_TOT_NUM_PIXELS; - } - if(flag_options->USE_MINI_HALOS){ - prev_deltax_unfiltered[ct] /= (HII_TOT_NUM_PIXELS+0.0); - log10_Mturnover_unfiltered[ct] /= (HII_TOT_NUM_PIXELS+0.0); - log10_Mturnover_MINI_unfiltered[ct] /= (HII_TOT_NUM_PIXELS+0.0); - } + //HMF integral initialisation + if(user_params->USE_INTERPOLATION_TABLES) { + if(user_params->INTEGRATION_METHOD_ATOMIC == 2 || user_params->INTEGRATION_METHOD_MINI == 2) + initialiseSigmaMInterpTable(fmin(MMIN_FAST,ionbox_constants.M_min),1e20); + else + initialiseSigmaMInterpTable(ionbox_constants.M_min,1e20); + } + + if(user_params->INTEGRATION_METHOD_ATOMIC == 1 || (flag_options->USE_MINI_HALOS && user_params->INTEGRATION_METHOD_MINI == 1)) + initialise_GL(ionbox_constants.lnMmin, ionbox_constants.lnMmax_gl); + + double f_limit_acg; + double f_limit_mcg; + set_mean_fcoll(&ionbox_constants,previous_ionize_box,box,Mturnover_global_avg,Mturnover_global_avg_MINI,&f_limit_acg,&f_limit_mcg); + double exp_global_hii = box->mean_f_coll * ionbox_constants.ion_eff_factor_gl \ + + box->mean_f_coll_MINI * ionbox_constants.ion_eff_factor_mini_gl; + + //TODO: change this from an if-else to an early-exit / cleanup call + if(exp_global_hii < global_params.HII_ROUND_ERR){ // way too small to ionize anything... + LOG_DEBUG("Mean collapsed fraciton %.3e too small to ionize, stopping early",exp_global_hii); + global_xH = set_fully_neutral_box(box,spin_temp,perturbed_field,&ionbox_constants); + } + else{ + //DO THE R2C TRANSFORMS + //TODO: add debug average printing to these boxes + //TODO: put a flag for to turn off clipping instead of putting the wide limits + prepare_box_for_filtering(perturbed_field->density,grid_struct->deltax_unfiltered,ionbox_constants.photoncons_adjustment_factor,-1.,1e6); + if(flag_options->USE_HALO_FIELD) { + prepare_box_for_filtering(halos->n_ion,grid_struct->stars_unfiltered,1.,0.,1e20); + prepare_box_for_filtering(halos->whalo_sfr,grid_struct->sfr_unfiltered,1.,0.,1e20); + } + else{ + if(flag_options->USE_MINI_HALOS){ + prepare_box_for_filtering(previous_perturbed_field->density,grid_struct->prev_deltax_unfiltered,1.,-1,1e6); + //since the turnover mass boxes were assigned separately (they needed more complex functions)... + dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_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); } } + if(flag_options->USE_TS_FLUCT){ + prepare_box_for_filtering(spin_temp->x_e_box,grid_struct->xe_unfiltered,1.,0,1.); + } - LOG_SUPER_DEBUG("deltax unfiltered calculated"); + if(ionbox_constants.filter_recombinations){ + prepare_box_for_filtering(previous_ionize_box->dNrec_box,grid_struct->N_rec_unfiltered,1.,0,1e20); + } + LOG_SUPER_DEBUG("FFTs performed"); // ************************************************************************************* // // ***************** LOOP THROUGH THE FILTER RADII (in Mpc) *************************** // // ************************************************************************************* // // set the max radius we will use, making sure we are always sampling the same values of radius // (this avoids aliasing differences w redshift) - R=fmax(global_params.R_BUBBLE_MIN, (cell_length_factor*user_params->BOX_LEN/(float)user_params->HII_DIM)); - while ((R - fmin(astro_params->R_BUBBLE_MAX, L_FACTOR * user_params->BOX_LEN)) <= FRACT_FLOAT_ERR) { - R *= global_params.DELTA_R_HII_FACTOR; - if (R >= fmin(astro_params->R_BUBBLE_MAX, L_FACTOR * user_params->BOX_LEN)) { - stored_R = R / (global_params.DELTA_R_HII_FACTOR); + int R_ct; + struct RadiusSpec curr_radius; + for(R_ct=n_radii;R_ct--;){ + curr_radius = radii_spec[R_ct]; + + //TODO: As far as I can tell, This was the previous behaviour with the while loop + // So if the cell size is smaller than the minimum mass (rare) we still filter the last step + // and don't assign any partial ionisaitons + if(ionbox_constants.M_min > RtoM(curr_radius.R)){ + LOG_DEBUG("Radius %.2e Mass %.2e smaller than minimum %.2e, stopping...", + radii_spec[R_ct].R,radii_spec[R_ct].M_max_R,ionbox_constants.M_min); + break; } - } - - LOG_DEBUG("set max radius: %f", R); + LOG_ULTRA_DEBUG("while loop for until RtoM(R)=%f reaches M_MIN=%f", RtoM(curr_radius.R), ionbox_constants.M_min); - R=fmin(astro_params->R_BUBBLE_MAX, L_FACTOR*user_params->BOX_LEN); + //do all the filtering and inverse transform + copy_filter_transform(grid_struct,&ionbox_constants,curr_radius); - LAST_FILTER_STEP = 0; - - first_step_R = 1; - - counter = 0; - - while (!LAST_FILTER_STEP && (M_MIN < RtoM(R)) ){ - LOG_ULTRA_DEBUG("while loop for until RtoM(R)=%f reaches M_MIN=%f", RtoM(R), M_MIN); - - // Check if we are the last filter step - if ( ((R/(global_params.DELTA_R_HII_FACTOR) - cell_length_factor*(user_params->BOX_LEN)/(float)(user_params->HII_DIM)) <= FRACT_FLOAT_ERR) || \ - ((R/(global_params.DELTA_R_HII_FACTOR) - global_params.R_BUBBLE_MIN) <= FRACT_FLOAT_ERR) ) { - LAST_FILTER_STEP = 1; - R = fmax(cell_length_factor*user_params->BOX_LEN/(double)(user_params->HII_DIM), global_params.R_BUBBLE_MIN); - } - - // Copy all relevant quantities from memory into new arrays to be smoothed and FFT'd. - if (flag_options->USE_TS_FLUCT) { - memcpy(xe_filtered, xe_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); - } - if (recomb_filter_flag) { - memcpy(N_rec_filtered, N_rec_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); - } - if (flag_options->USE_HALO_FIELD) { - memcpy(stars_filtered, stars_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); - memcpy(sfr_filtered, sfr_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); - } - - memcpy(deltax_filtered, deltax_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - - if(flag_options->USE_MINI_HALOS){ - memcpy(prev_deltax_filtered, prev_deltax_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - memcpy(log10_Mturnover_MINI_filtered, log10_Mturnover_MINI_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - memcpy(log10_Mturnover_filtered, log10_Mturnover_unfiltered, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS); - } - - if (!LAST_FILTER_STEP || - ((R - cell_length_factor * (user_params->BOX_LEN / (double) (user_params->HII_DIM))) > - FRACT_FLOAT_ERR)) { - if (flag_options->USE_TS_FLUCT) { - filter_box(xe_filtered, 1, global_params.HII_FILTER, R, 0.); - } - if (recomb_filter_flag) { - filter_box(N_rec_filtered, 1, global_params.HII_FILTER, R, 0.); - } - if (flag_options->USE_HALO_FIELD) { - int filter_hf = flag_options->USE_EXP_FILTER ? 3 : global_params.HII_FILTER; - filter_box(stars_filtered, 1, filter_hf, R, exp_mfp); - filter_box(sfr_filtered, 1, filter_hf, R, exp_mfp); - } - filter_box(deltax_filtered, 1, global_params.HII_FILTER, R, 0.); - if(flag_options->USE_MINI_HALOS){ - filter_box(prev_deltax_filtered, 1, global_params.HII_FILTER, R, 0.); - filter_box(log10_Mturnover_MINI_filtered, 1, global_params.HII_FILTER, R, 0.); - filter_box(log10_Mturnover_filtered, 1, global_params.HII_FILTER, R, 0.); - } - } - - // Perform FFTs - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, deltax_filtered); - - if(flag_options->USE_MINI_HALOS){ - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, prev_deltax_filtered); - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, log10_Mturnover_MINI_filtered); - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, log10_Mturnover_filtered); - } - - if (flag_options->USE_HALO_FIELD) { - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, stars_filtered); - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, sfr_filtered); - } - - if (flag_options->USE_TS_FLUCT) { - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, xe_filtered); - } - - if (recomb_filter_flag) { - dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, N_rec_filtered); - } - - // Check if this is the last filtering scale. If so, we don't need deltax_unfiltered anymore. - // We will re-read it to get the real-space field, which we will use to set the residual neutral fraction - ST_over_PS = 0; - ST_over_PS_MINI = 0; - f_coll = 0; - f_coll_MINI = 0; - massofscaleR = RtoM(R); - lnM_cond = log(massofscaleR); - - sigmaMmax = sigma_z0(massofscaleR); + bool need_prev_ion = flag_options->USE_MINI_HALOS && \ + (previous_ionize_box->mean_f_coll_MINI * ionbox_constants.ion_eff_factor_mini_gl + \ + previous_ionize_box->mean_f_coll * ionbox_constants.ion_eff_factor_gl > 1e-4); if (!flag_options->USE_HALO_FIELD) { - if(user_params->INTEGRATION_METHOD_ATOMIC == 1 || (flag_options->USE_MINI_HALOS && user_params->INTEGRATION_METHOD_MINI == 1)) - initialise_GL(lnMmin, lnM_cond); - if (flag_options->USE_MASS_DEPENDENT_ZETA) { - - min_density = max_density = 0.0; - -#pragma omp parallel shared(deltax_filtered) private(x, y, z) num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(max:max_density) reduction(min:min_density) - for (x = 0; x < user_params->HII_DIM; x++) { - for (y = 0; y < user_params->HII_DIM; y++) { - for (z = 0; z < HII_D_PARA; z++) { - // delta cannot be less than -1 - *((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)) = fmaxf( - *((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)), -1. + FRACT_FLOAT_ERR); - - if (*((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)) < min_density) { - min_density = *((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)); - } - if (*((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)) > max_density) { - max_density = *((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)); - } - } - } - } - } - - if (flag_options->USE_MINI_HALOS){ - // do the same for prev - prev_min_density = prev_max_density = 0.0; -#pragma omp parallel shared(prev_deltax_filtered) private(x, y, z) num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(max:prev_max_density) reduction(min:prev_min_density) - for (x=0; xHII_DIM; x++){ - for (y=0; yHII_DIM; y++){ - for (z=0; z prev_max_density ) { - prev_max_density = *((float *)prev_deltax_filtered + HII_R_FFT_INDEX(x,y,z)); - } - } - } - } - } - - // do the same for logM - log10Mturn_min = 999; - log10Mturn_max = 0.0; - log10Mturn_min_MINI = 999; - log10Mturn_max_MINI = 0.0; - -#pragma omp parallel shared(log10_Mturnover_filtered,log10_Mturnover_MINI_filtered,log10_Mcrit_atom,log10_Mcrit_mol) private(x, y, z) num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(max:log10Mturn_max,log10Mturn_max_MINI) reduction(min:log10Mturn_min,log10Mturn_min_MINI) - for (x=0; xHII_DIM; x++){ - for (y=0; yHII_DIM; y++){ - for (z=0; z LOG10_MTURN_MAX) - *((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z)) = LOG10_MTURN_MAX; - // Mturnover cannot be less than Mcrit_mol - if (*((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)) < log10_Mcrit_mol) - *((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)) = log10_Mcrit_mol; - if (*((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)) > LOG10_MTURN_MAX) - *((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)) = LOG10_MTURN_MAX; - - if (*((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z)) < log10Mturn_min) - log10Mturn_min = *((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z)); - if (*((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z)) > log10Mturn_max) - log10Mturn_max = *((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z)); - if (*((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)) < log10Mturn_min_MINI) - log10Mturn_min_MINI = *((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)); - if (*((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)) > log10Mturn_max_MINI) - log10Mturn_max_MINI = *((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)); - } - } - } - } - - if(user_params->USE_INTERPOLATION_TABLES) { - log10Mturn_min = log10Mturn_min *0.99; - log10Mturn_max = log10Mturn_max *1.01; - log10Mturn_min_MINI = log10Mturn_min_MINI *0.99; - log10Mturn_max_MINI = log10Mturn_max_MINI *1.01; - } - } - - LOG_ULTRA_DEBUG("Tb limits d (%.2e,%.2e), m (%.2e,%.2e) t (%.2e,%.2e) tm (%.2e,%.2e)", - min_density,max_density,M_MIN,massofscaleR,log10Mturn_min,log10Mturn_max, - log10Mturn_min_MINI,log10Mturn_max_MINI); - if(user_params->USE_INTERPOLATION_TABLES){ - //To avoid both zero bin widths and densest cell segfault - min_density -= 0.001; - max_density += 0.001; - prev_min_density -= 0.001; - prev_max_density += 0.001; - initialise_Nion_Conditional_spline(redshift,Mcrit_atom,min_density,max_density,M_MIN,massofscaleR,massofscaleR, - log10Mturn_min,log10Mturn_max,log10Mturn_min_MINI,log10Mturn_max_MINI, - astro_params->ALPHA_STAR, astro_params->ALPHA_STAR_MINI, - alpha_esc_var,astro_params->F_STAR10, - norm_esc_var,Mlim_Fstar,Mlim_Fesc,astro_params->F_STAR7_MINI, - astro_params->F_ESC7_MINI, Mlim_Fstar_MINI, Mlim_Fesc_MINI, - user_params->INTEGRATION_METHOD_ATOMIC, user_params->INTEGRATION_METHOD_MINI, - flag_options->USE_MINI_HALOS,false); - - if(flag_options->USE_MINI_HALOS && (previous_ionize_box->mean_f_coll_MINI * ION_EFF_FACTOR_MINI + previous_ionize_box->mean_f_coll * ION_EFF_FACTOR > 1e-4)){ - initialise_Nion_Conditional_spline(prev_redshift,Mcrit_atom,prev_min_density,prev_max_density,M_MIN,massofscaleR,massofscaleR, - log10Mturn_min,log10Mturn_max,log10Mturn_min_MINI,log10Mturn_max_MINI, - astro_params->ALPHA_STAR, astro_params->ALPHA_STAR_MINI, - alpha_esc_var,astro_params->F_STAR10, - norm_esc_var,Mlim_Fstar,Mlim_Fesc,astro_params->F_STAR7_MINI, - astro_params->F_ESC7_MINI,Mlim_Fstar_MINI, Mlim_Fesc_MINI, - user_params->INTEGRATION_METHOD_ATOMIC, user_params->INTEGRATION_METHOD_MINI, - flag_options->USE_MINI_HALOS,true); - } - } - } - else { - erfc_denom = 2. * (pow(sigma_z0(M_MIN), 2) - pow(sigma_z0(massofscaleR), 2)); - if (erfc_denom < 0) { // our filtering scale has become too small - break; - } - erfc_denom = sqrt(erfc_denom); - erfc_denom = 1. / (growth_factor * erfc_denom); - } - } - LOG_SUPER_DEBUG("Initialised tables"); - // Determine the global averaged f_coll for the overall normalisation - - // Reset value of int check to see if we are over-stepping our interpolation table - for (i = 0; i < user_params->N_THREADS; i++) { - overdense_int_boundexceeded_threaded[i] = 0; - } - - // renormalize the collapse fraction so that the mean matches ST, - // since we are using the evolved (non-linear) density field -#pragma omp parallel shared(deltax_filtered,N_rec_filtered,xe_filtered,overdense_int_boundexceeded_threaded,erfc_denom,erfc_arg_min,\ - erfc_arg_max,InvArgBinWidth,ArgBinWidth,ERFC_VALS_DIFF,ERFC_VALS,log10_Mturnover_filtered,log10Mturn_min, \ - log10_Mturnover_MINI_filtered,prev_deltax_filtered,previous_ionize_box,ION_EFF_FACTOR,\ - box,counter,stars_filtered,massofscaleR,sigmaMmax,\ - M_MIN,growth_factor,Mlim_Fstar,Mlim_Fesc,Mcrit_atom,Mlim_Fstar_MINI,Mlim_Fesc_MINI,prev_growth_factor) \ - private(x,y,z,curr_dens,Splined_Fcoll,Splined_Fcoll_MINI,erfc_arg_val,erfc_arg_val_index,log10_Mturnover,\ - log10_Mturnover_MINI,prev_dens,prev_Splined_Fcoll,prev_Splined_Fcoll_MINI) \ - num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(+:f_coll,f_coll_MINI) - for (x = 0; x < user_params->HII_DIM; x++) { - for (y = 0; y < user_params->HII_DIM; y++) { - for (z = 0; z < HII_D_PARA; z++) { - - // delta cannot be less than -1 - *((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)) = fmaxf( - *((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)), -1. + FRACT_FLOAT_ERR); - - // cannot be less than zero - if (recomb_filter_flag) { - *((float *) N_rec_filtered + HII_R_FFT_INDEX(x, y, z)) = fmaxf(*((float *) N_rec_filtered + HII_R_FFT_INDEX(x, y, z)), 0.0); - } - - // x_e has to be between zero and unity - if (flag_options->USE_TS_FLUCT) { - *((float *) xe_filtered + HII_R_FFT_INDEX(x, y, z)) = fmaxf(*((float *) xe_filtered + HII_R_FFT_INDEX(x, y, z)), 0.); - *((float *) xe_filtered + HII_R_FFT_INDEX(x, y, z)) = fminf(*((float *) xe_filtered + HII_R_FFT_INDEX(x, y, z)), 0.999); - } - - // stellar mass & sfr cannot be less than zero - if(flag_options->USE_HALO_FIELD) { - *((float *)stars_filtered + HII_R_FFT_INDEX(x,y,z)) = fmaxf( - *((float *)stars_filtered + HII_R_FFT_INDEX(x,y,z)) , 0.0); - *((float *)sfr_filtered + HII_R_FFT_INDEX(x,y,z)) = fmaxf( - *((float *)sfr_filtered + HII_R_FFT_INDEX(x,y,z)) , 0.0); - - //Ionising photon output - Splined_Fcoll = *((float *)stars_filtered + HII_R_FFT_INDEX(x,y,z)); - //Minihalos are taken care of already - Splined_Fcoll_MINI = 0.; - //The smoothing done with minihalos corrects for sudden changes in M_crit - //Nion_smoothed(z,Mcrit) = Nion(z,Mcrit) + (Nion(z_prev,Mcrit_prev) - Nion(z_prev,Mcrit)) - prev_Splined_Fcoll = 0.; - prev_Splined_Fcoll_MINI = 0.; - } - else { - curr_dens = *((float *) deltax_filtered + HII_R_FFT_INDEX(x, y, z)); - if (flag_options->USE_MASS_DEPENDENT_ZETA) { - if (flag_options->USE_MINI_HALOS){ - log10_Mturnover = *((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z)); - log10_Mturnover_MINI = *((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z)); - - Splined_Fcoll_MINI = EvaluateNion_Conditional_MINI(curr_dens,log10_Mturnover_MINI,growth_factor,M_MIN, - massofscaleR,massofscaleR,sigmaMmax,Mcrit_atom,Mlim_Fstar,Mlim_Fesc,false); - - prev_dens = *((float *)prev_deltax_filtered + HII_R_FFT_INDEX(x,y,z)); - - if (previous_ionize_box->mean_f_coll_MINI * ION_EFF_FACTOR_MINI + previous_ionize_box->mean_f_coll * ION_EFF_FACTOR > 1e-4){ - prev_Splined_Fcoll = EvaluateNion_Conditional(prev_dens,log10_Mturnover,prev_growth_factor,M_MIN,massofscaleR, - massofscaleR,sigmaMmax,Mlim_Fstar,Mlim_Fesc,true); - prev_Splined_Fcoll_MINI = EvaluateNion_Conditional_MINI(prev_dens,log10_Mturnover_MINI,prev_growth_factor,M_MIN, - massofscaleR,massofscaleR,sigmaMmax,Mcrit_atom,Mlim_Fstar,Mlim_Fesc,true); - } - else{ - prev_Splined_Fcoll = 0.; - prev_Splined_Fcoll_MINI = 0.; - } - } - else{ - log10_Mturnover = log10(astro_params->M_TURN); - } - Splined_Fcoll = EvaluateNion_Conditional(curr_dens,log10_Mturnover,growth_factor,M_MIN,massofscaleR,massofscaleR,sigmaMmax,Mlim_Fstar,Mlim_Fesc,false); - } - else{ - erfc_arg_val = (Deltac - curr_dens) * erfc_denom; - if (erfc_arg_val < erfc_arg_min || erfc_arg_val > erfc_arg_max) { - Splined_Fcoll = splined_erfc(erfc_arg_val); - } else { - erfc_arg_val_index = (int) floor((erfc_arg_val - erfc_arg_min) * InvArgBinWidth); - - Splined_Fcoll = ERFC_VALS[erfc_arg_val_index] + \ - (erfc_arg_val - (erfc_arg_min + ArgBinWidth * (double) erfc_arg_val_index)) * ERFC_VALS_DIFF[erfc_arg_val_index] *InvArgBinWidth; - } - } - } - // save the value of the collasped fraction into the Fcoll array - if (flag_options->USE_MINI_HALOS && !flag_options->USE_HALO_FIELD){ - if (Splined_Fcoll > 1.) Splined_Fcoll = 1.; - if (Splined_Fcoll < 0.) Splined_Fcoll = 1e-40; - if (prev_Splined_Fcoll > 1.) prev_Splined_Fcoll = 1.; - if (prev_Splined_Fcoll < 0.) prev_Splined_Fcoll = 1e-40; - 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)] + Splined_Fcoll - prev_Splined_Fcoll; - - if (box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] >1.) box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1.; - //if (box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] <0.) box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1e-40; - //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 += box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; - if(isfinite(f_coll)==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[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)],\ - Splined_Fcoll, prev_Splined_Fcoll, curr_dens, prev_dens, \ - log10_Mturnover, *((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z))); -// Throw(ParameterError); - Throw(InfinityorNaNError); - } - - if (Splined_Fcoll_MINI > 1.) Splined_Fcoll_MINI = 1.; - if (Splined_Fcoll_MINI < 0.) Splined_Fcoll_MINI = 1e-40; - if (prev_Splined_Fcoll_MINI > 1.) prev_Splined_Fcoll_MINI = 1.; - if (prev_Splined_Fcoll_MINI < 0.) prev_Splined_Fcoll_MINI = 1e-40; - 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)] + Splined_Fcoll_MINI - prev_Splined_Fcoll_MINI; - - if (box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] >1.) box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1.; - //if (box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] <0.) box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = 1e-40; - //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 += box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; - if(isfinite(f_coll_MINI)==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[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)],\ - Splined_Fcoll_MINI, prev_Splined_Fcoll_MINI, log10_Mturnover_MINI,\ - *((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z))); - LOG_DEBUG("%g,%g",previous_ionize_box->Fcoll[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)]); - Throw(InfinityorNaNError); - } - } - else{ - box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)] = Splined_Fcoll; - f_coll += Splined_Fcoll; - } - } - } - } - } // end loop through Fcoll box - - for (i = 0; i < user_params->N_THREADS; i++) { - if (overdense_int_boundexceeded_threaded[i] == 1) { - LOG_ERROR("I have overstepped my allocated memory for one of the interpolation tables for the nion_splines"); -// Throw(ParameterError); - Throw(TableEvaluationError); - } - } - if(isfinite(f_coll)==0) { - LOG_ERROR("f_coll is either infinite or NaN!"); -// Throw(ParameterError); - Throw(InfinityorNaNError); + setup_integration_tables(grid_struct,&ionbox_constants,curr_radius,need_prev_ion); } - f_coll /= (double) HII_TOT_NUM_PIXELS; - - if(isfinite(f_coll_MINI)==0) { - LOG_ERROR("f_coll_MINI is either infinite or NaN!"); -// Throw(ParameterError); - Throw(InfinityorNaNError); - } - - f_coll_MINI /= (double) HII_TOT_NUM_PIXELS; + calculate_fcoll_grid(box,previous_ionize_box,grid_struct,&ionbox_constants,&curr_radius); // To avoid ST_over_PS becoming nan when f_coll = 0, I set f_coll = FRACT_FLOAT_ERR. - if (flag_options->USE_MASS_DEPENDENT_ZETA) { - if (f_coll <= f_coll_min) f_coll = f_coll_min; + // TODO: This was the previous behaviour, but is this right? + // setting the *total* to the minimum for the adjustment factor, + // then clipping the grid in the loop below? + if (flag_options_global->USE_MASS_DEPENDENT_ZETA) { + if (curr_radius.f_coll_grid_mean <= f_limit_acg) curr_radius.f_coll_grid_mean = f_limit_acg; if (flag_options->USE_MINI_HALOS){ - if (f_coll_MINI <= f_coll_min_MINI) f_coll_MINI = f_coll_min_MINI; + if (curr_radius.f_coll_grid_mean_MINI <= f_limit_mcg) curr_radius.f_coll_grid_mean_MINI = f_limit_mcg; } } else { - if (f_coll <= FRACT_FLOAT_ERR) f_coll = FRACT_FLOAT_ERR; - } - - if(flag_options->USE_HALO_FIELD){ - ST_over_PS = 1.; - ST_over_PS_MINI = 1.; - } - else{ - ST_over_PS = box->mean_f_coll/f_coll; - ST_over_PS_MINI = box->mean_f_coll_MINI/f_coll_MINI; - LOG_SUPER_DEBUG("global mean fcoll %.4e box mean fcoll %.4e ratio %.4e",box->mean_f_coll,f_coll,ST_over_PS); - LOG_SUPER_DEBUG("MINI: global mean fcoll %.4e box mean fcoll %.4e ratio %.4e",box->mean_f_coll_MINI,f_coll_MINI,ST_over_PS_MINI); - } - - Gamma_R_prefactor = (R*CMperMPC) * SIGMA_HI * global_params.ALPHA_UVB / (global_params.ALPHA_UVB+2.75) * N_b0 * ION_EFF_FACTOR / 1.0e-12; - Gamma_R_prefactor_MINI = (R*CMperMPC) * SIGMA_HI * global_params.ALPHA_UVB / (global_params.ALPHA_UVB+2.75) * N_b0 * ION_EFF_FACTOR_MINI / 1.0e-12; - if(flag_options->PHOTON_CONS_TYPE == 1) { - // Used for recombinations, which means we want to use the original redshift not the adjusted redshift - Gamma_R_prefactor *= pow(1+stored_redshift, 2); - Gamma_R_prefactor_MINI *= pow(1+stored_redshift, 2); - } - else { - Gamma_R_prefactor *= pow(1+redshift, 2); - Gamma_R_prefactor_MINI *= pow(1+redshift, 2); + if (curr_radius.f_coll_grid_mean <= FRACT_FLOAT_ERR) curr_radius.f_coll_grid_mean = FRACT_FLOAT_ERR; } - //With the halo field, we use the filtered, f_esc and N_ion weighted star formation rate, which should be equivalent to - // `Fcoll` * OMb * RHOcrit * (1+delta) in the no halo field case. we also need a factor of 1/(1+delta) later on - // to match that in the recombination (`Fcoll` is effectively fesc*star per baryon, whereas the filtered grids are fesc*SFRD) - //So the halo option needs an extra density term and the nonhalo option needs the SFR term - if(flag_options->USE_HALO_FIELD){ - Gamma_R_prefactor /= RHOcrit * cosmo_params->OMb; - Gamma_R_prefactor_MINI /= RHOcrit * cosmo_params->OMb; - } - else{ - //Minihalos already included - Gamma_R_prefactor /= t_ast; - Gamma_R_prefactor_MINI /= t_ast; - } - - if (global_params.FIND_BUBBLE_ALGORITHM != 2 && global_params.FIND_BUBBLE_ALGORITHM != 1) { // center method - LOG_ERROR("Incorrect choice of find bubble algorithm: %i", - global_params.FIND_BUBBLE_ALGORITHM); - Throw(ValueError); - } - -#pragma omp parallel shared(deltax_filtered,N_rec_filtered,xe_filtered,box,ST_over_PS,pixel_mass,M_MIN,r,f_coll_min,Gamma_R_prefactor,\ - ION_EFF_FACTOR,ION_EFF_FACTOR_MINI,LAST_FILTER_STEP,counter,ST_over_PS_MINI,f_coll_min_MINI,Gamma_R_prefactor_MINI,TK,cT_ad,perturbed_field) \ - private(x,y,z,curr_dens,Splined_Fcoll,f_coll,ave_M_coll_cell,ave_N_min_cell,N_halos_in_cell,rec,xHII_from_xrays,\ - Splined_Fcoll_MINI,f_coll_MINI, res_xH) \ - num_threads(user_params->N_THREADS) - { -#pragma omp for - for (x = 0; x < user_params->HII_DIM; x++) { - for (y = 0; y < user_params->HII_DIM; y++) { - for (z = 0; z < HII_D_PARA; z++) { - - //TODO: There may be a reason to use unfiltered density for CELL_RECOMB case, - // since the "Fcoll" represents photons reaching the central cell rather than - // photons in the entire sphere (See issue #434) - // if(flag_options->CELL_RECOMB) - // curr_dens = perturbed_field->density[HII_R_INDEX(x,y,z)]; - // else - curr_dens = *((float *)deltax_filtered + HII_R_FFT_INDEX(x,y,z)); - - Splined_Fcoll = box->Fcoll[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; - - //Since the halo boxes give ionising photon output, this term accounts for the local density of absorbers - // We have separated the source/absorber filtering in the halo model so this is necessary - if(flag_options->USE_HALO_FIELD){ - Splined_Fcoll *= 1 / (RHOcrit*cosmo_params->OMb*(1+curr_dens)); - } - - f_coll = ST_over_PS * Splined_Fcoll; - - //MINIHALOS are already included in the halo model - if (flag_options->USE_MINI_HALOS && !flag_options->USE_HALO_FIELD){ - Splined_Fcoll_MINI = box->Fcoll_MINI[counter * HII_TOT_NUM_PIXELS + HII_R_INDEX(x,y,z)]; - f_coll_MINI = ST_over_PS_MINI * Splined_Fcoll_MINI; - } - else{ - f_coll_MINI = 0.; - } - - if (LAST_FILTER_STEP){ - ave_M_coll_cell = (f_coll + f_coll_MINI) * pixel_mass * (1. + curr_dens); - ave_N_min_cell = ave_M_coll_cell / M_MIN; // ave # of M_MIN halos in cell - if(user_params->NO_RNG) { - N_halos_in_cell = 1.; - } - else { - N_halos_in_cell = (int) gsl_ran_poisson(r[omp_get_thread_num()], - global_params.N_POISSON); - } - } - - if (flag_options->USE_MASS_DEPENDENT_ZETA) { - if (f_coll <= f_coll_min) f_coll = f_coll_min; - if (flag_options->USE_MINI_HALOS){ - if (f_coll_MINI <= f_coll_min_MINI) f_coll_MINI = f_coll_min_MINI; - } - } - - if (flag_options->INHOMO_RECO) { - if(flag_options->CELL_RECOMB) - rec = previous_ionize_box->dNrec_box[HII_R_INDEX(x,y,z)]; - else - rec = (*((float *) N_rec_filtered + HII_R_FFT_INDEX(x, y, z))); // number of recombinations per mean baryon - - rec /= (1. + curr_dens); // number of recombinations per baryon inside cell/filter - } - else { - rec = 0.; - } - - // adjust the denominator of the collapse fraction for the residual electron fraction in the neutral medium - if (flag_options->USE_TS_FLUCT){ - xHII_from_xrays = *((float *)xe_filtered + HII_R_FFT_INDEX(x,y,z)); - } else { - xHII_from_xrays = 0.; - } - - if(x+y+z == 0 && !flag_options->USE_HALO_FIELD){ - //reusing variables (i know its not log10) - if(flag_options->USE_MINI_HALOS){ - log10_Mturnover = pow(10,*((float *)log10_Mturnover_filtered + HII_R_FFT_INDEX(x,y,z))); - log10_Mturnover_MINI = pow(10,*((float *)log10_Mturnover_MINI_filtered + HII_R_FFT_INDEX(x,y,z))); - } - else{ - log10_Mturnover = astro_params->M_TURN; - } - LOG_SUPER_DEBUG("Cell 0: R=%.1f | d %.4e | fcoll (s %.4e f %.4e i %.4e) | rec %.4e | X %.4e", - R,curr_dens,Splined_Fcoll,f_coll,\ - Nion_ConditionalM(growth_factor,lnMmin,lnM_cond,massofscaleR,sigmaMmax,curr_dens, - log10_Mturnover, - astro_params->ALPHA_STAR,astro_params->ALPHA_ESC,astro_params->F_STAR10, - astro_params->F_ESC10,Mlim_Fstar,Mlim_Fesc,user_params->INTEGRATION_METHOD_ATOMIC),rec,xHII_from_xrays); - if(flag_options->USE_MINI_HALOS){ - LOG_SUPER_DEBUG("Mini (s %.4e f %.4e i %.4e)",Splined_Fcoll_MINI,f_coll_MINI,\ - Nion_ConditionalM_MINI(growth_factor,lnMmin,lnM_cond,massofscaleR,sigmaMmax,curr_dens, - log10_Mturnover_MINI,Mcrit_atom, - astro_params->ALPHA_STAR_MINI,astro_params->ALPHA_ESC,astro_params->F_STAR7_MINI, - astro_params->F_ESC7_MINI,Mlim_Fstar,Mlim_Fesc,user_params->INTEGRATION_METHOD_MINI)); - } - } - - // check if fully ionized! - if ( (f_coll * ION_EFF_FACTOR + f_coll_MINI * ION_EFF_FACTOR_MINI > (1. - xHII_from_xrays)*(1.0+rec)) ){ //IONIZED!! - // if this is the first crossing of the ionization barrier for this cell (largest R), record the gamma - // this assumes photon-starved growth of HII regions... breaks down post EoR - if (flag_options->INHOMO_RECO && (box->xH_box[HII_R_INDEX(x,y,z)] > FRACT_FLOAT_ERR) ){ - if(flag_options->USE_HALO_FIELD){ - box->Gamma12_box[HII_R_INDEX(x,y,z)] = Gamma_R_prefactor / (1+curr_dens) * (*((float *)sfr_filtered + HII_R_FFT_INDEX(x,y,z))); - } - else{ - box->Gamma12_box[HII_R_INDEX(x,y,z)] = Gamma_R_prefactor * f_coll + Gamma_R_prefactor_MINI * f_coll_MINI; - } - box->MFP_box[HII_R_INDEX(x,y,z)] = R; - } - - // keep track of the first time this cell is ionized (earliest time) - if (previous_ionize_box->z_re_box[HII_R_INDEX(x,y,z)] < 0){ - box->z_re_box[HII_R_INDEX(x,y,z)] = redshift; - } else{ - box->z_re_box[HII_R_INDEX(x,y,z)] = previous_ionize_box->z_re_box[HII_R_INDEX(x,y,z)]; - } + find_ionised_regions(box,previous_ionize_box,perturbed_field,spin_temp,curr_radius,&ionbox_constants, + grid_struct,f_limit_acg,f_limit_mcg); - // FLAG CELL(S) AS IONIZED - if (global_params.FIND_BUBBLE_ALGORITHM == 2) // center method - box->xH_box[HII_R_INDEX(x,y,z)] = 0; - if (global_params.FIND_BUBBLE_ALGORITHM == 1) // sphere method - update_in_sphere(box->xH_box, user_params->HII_DIM, HII_D_PARA, R/(user_params->BOX_LEN), \ - x/(user_params->HII_DIM+0.0), y/(user_params->HII_DIM+0.0), z/(HII_D_PARA+0.0)); - } // end ionized - // If not fully ionized, then assign partial ionizations - else if (LAST_FILTER_STEP && (box->xH_box[HII_R_INDEX(x, y, z)] > TINY)) { - - if (f_coll>1) f_coll=1; - if (f_coll_MINI>1) f_coll_MINI=1; - - if (!flag_options->USE_HALO_FIELD){ - if(ave_N_min_cell < global_params.N_POISSON) { - f_coll = N_halos_in_cell * ( ave_M_coll_cell / (float)global_params.N_POISSON ) / (pixel_mass*(1. + curr_dens)); - if (flag_options->USE_MINI_HALOS){ - f_coll_MINI = f_coll * (f_coll_MINI * ION_EFF_FACTOR_MINI) / (f_coll * ION_EFF_FACTOR + f_coll_MINI * ION_EFF_FACTOR_MINI); - f_coll = f_coll - f_coll_MINI; - } - else{ - f_coll_MINI = 0.; - } - } - - if (ave_M_coll_cell < (M_MIN / 5.)) { - f_coll = 0.; - f_coll_MINI = 0.; - } - } - - if (f_coll>1) f_coll=1; - if (f_coll_MINI>1) f_coll_MINI=1; - res_xH = 1. - f_coll * ION_EFF_FACTOR - f_coll_MINI * ION_EFF_FACTOR_MINI; - // put the partial ionization here because we need to exclude xHII_from_xrays... - if (flag_options->USE_TS_FLUCT){ - box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = ComputePartiallyIoinizedTemperature(spin_temp->Tk_box[HII_R_INDEX(x,y,z)], res_xH); - } - else{ - box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = ComputePartiallyIoinizedTemperature(TK*(1 + cT_ad*perturbed_field->density[HII_R_INDEX(x,y,z)]), res_xH); - } - res_xH -= xHII_from_xrays; - - // and make sure fraction doesn't blow up for underdense pixels - if (res_xH < 0) - res_xH = 0; - else if (res_xH > 1) - res_xH = 1; - - box->xH_box[HII_R_INDEX(x, y, z)] = res_xH; - - } // end partial ionizations at last filtering step - } // k - } // j - } // i - } - - LOG_SUPER_DEBUG("z_re_box after R=%f: ", R); +#if LOG_LEVEL >= ULTRA_DEBUG_LEVEL + LOG_ULTRA_DEBUG("z_re_box after R=%f: ", curr_radius.R); debugSummarizeBox(box->z_re_box, user_params->HII_DIM, user_params->NON_CUBIC_FACTOR, " "); - - - if (first_step_R) { - R = stored_R; - first_step_R = 0; - } else { - R /= (global_params.DELTA_R_HII_FACTOR); - } - if (flag_options->USE_MINI_HALOS) - counter += 1; +#endif } + set_ionized_temperatures(box,perturbed_field,spin_temp,&ionbox_constants); + // find the neutral fraction +#if LOG_LEVEL >= DEBUG_LEVEL + global_xH = 0; -#pragma omp parallel shared(box,spin_temp,redshift,deltax_unfiltered_original,TK) private(x,y,z,thistk) num_threads(user_params->N_THREADS) + #pragma omp parallel shared(box) private(ct) num_threads(user_params->N_THREADS) { -#pragma omp for - for (x=0; xHII_DIM; x++){ - for (y=0; yHII_DIM; y++){ - for (z=0; zz_re_box[HII_R_INDEX(x,y,z)]>0) && (box->xH_box[HII_R_INDEX(x,y,z)] < TINY)){ - box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = ComputeFullyIoinizedTemperature(box->z_re_box[HII_R_INDEX(x,y,z)], \ - redshift, *((float *)deltax_unfiltered_original + HII_R_FFT_INDEX(x,y,z))); - // Below sometimes (very rare though) can happen when the density drops too fast and to below T_HI - if (flag_options->USE_TS_FLUCT){ - if (box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] < spin_temp->Tk_box[HII_R_INDEX(x,y,z)]) - box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = spin_temp->Tk_box[HII_R_INDEX(x,y,z)]; - } - else{ - thistk = TK*(1. + cT_ad*perturbed_field->density[HII_R_INDEX(x,y,z)]); - if (box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] < thistk) - box->temp_kinetic_all_gas[HII_R_INDEX(x,y,z)] = thistk; - } - } - } - } + #pragma omp for reduction(+:global_xH) + for (ct = 0; ct < HII_TOT_NUM_PIXELS; ct++) { + global_xH += box->xH_box[ct]; } } + global_xH /= (float) HII_TOT_NUM_PIXELS; +#endif - for (x=0; xHII_DIM; x++){ - for (y=0; yHII_DIM; y++){ - for (z=0; ztemp_kinetic_all_gas[HII_R_INDEX(x,y,z)])==0){ - LOG_ERROR("Tk after fully ioinzation is either infinite or a Nan. Something has gone wrong "\ - "in the temperature calculation: z_re=%.4f, redshift=%.4f, curr_dens=%.4e", box->z_re_box[HII_R_INDEX(x,y,z)], redshift, curr_dens); -// Throw(ParameterError); - Throw(InfinityorNaNError); - } - } - } - } - - // find the neutral fraction - if (LOG_LEVEL >= DEBUG_LEVEL) { - global_xH = 0; - -#pragma omp parallel shared(box) private(ct) num_threads(user_params->N_THREADS) - { -#pragma omp for reduction(+:global_xH) - for (ct = 0; ct < HII_TOT_NUM_PIXELS; ct++) { - global_xH += box->xH_box[ct]; - } - } - global_xH /= (float) HII_TOT_NUM_PIXELS; - } - - if (isfinite(global_xH) == 0) { + if(isfinite(global_xH) == 0) { LOG_ERROR( "Neutral fraction is either infinite or a Nan. Something has gone wrong in the ionisation calculation!"); -// Throw(ParameterError); Throw(InfinityorNaNError); } // update the N_rec field - if (flag_options->INHOMO_RECO) { - -#pragma omp parallel shared(perturbed_field, adjustment_factor, stored_redshift, redshift, box, previous_ionize_box, \ - fabs_dtdz, ZSTEP, something_finite_or_infinite) \ - private(x, y, z, curr_dens, z_eff, dNrec) num_threads(user_params->N_THREADS) - { -#pragma omp for - for (x = 0; x < user_params->HII_DIM; x++) { - for (y = 0; y < user_params->HII_DIM; y++) { - for (z = 0; z < HII_D_PARA; z++) { - - // use the original density and redshift for the snapshot (not the adjusted redshift) - // Only want to use the adjusted redshift for the ionisation field - //NOTE: but the structure field wasn't adjusted, this seems wrong - // curr_dens = 1.0 + (perturbed_field->density[HII_R_INDEX(x, y, z)]) / adjustment_factor; - curr_dens = 1.0 + (perturbed_field->density[HII_R_INDEX(x, y, z)]); - z_eff = pow(curr_dens, 1.0 / 3.0); - - if (flag_options->PHOTON_CONS_TYPE == 1) { - z_eff *= (1 + stored_redshift); - } else { - z_eff *= (1 + redshift); - } - - dNrec = splined_recombination_rate(z_eff - 1., box->Gamma12_box[HII_R_INDEX(x, y, z)]) * - fabs_dtdz * ZSTEP * (1. - box->xH_box[HII_R_INDEX(x, y, z)]); - - if (isfinite(dNrec) == 0) { - something_finite_or_infinite = 1; - } - - box->dNrec_box[HII_R_INDEX(x, y, z)] = - previous_ionize_box->dNrec_box[HII_R_INDEX(x, y, z)] + dNrec; - } - } - } - } - - if (something_finite_or_infinite) { - LOG_ERROR("Recombinations have returned either an infinite or NaN value."); -// Throw(ParameterError); - Throw(InfinityorNaNError); - } + if (flag_options->INHOMO_RECO){ + set_recombination_rates(box,previous_ionize_box,perturbed_field,&ionbox_constants); } fftwf_cleanup_threads(); @@ -1507,68 +1396,17 @@ LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll_MINI: %e", box->mean_f destruct_heat(); - for (i=0; iN_THREADS; i++) { - gsl_rng_free (r[i]); - } - -LOG_DEBUG("global_xH = %e",global_xH); - - fftwf_free(deltax_unfiltered); - fftwf_free(deltax_unfiltered_original); - fftwf_free(deltax_filtered); - if(flag_options->USE_MINI_HALOS){ - fftwf_free(prev_deltax_unfiltered); - fftwf_free(prev_deltax_filtered); - } - if(flag_options->USE_TS_FLUCT) { - fftwf_free(xe_unfiltered); - fftwf_free(xe_filtered); - } - if (recomb_filter_flag){ - fftwf_free(N_rec_unfiltered); - fftwf_free(N_rec_filtered); - } - - if(flag_options->USE_HALO_FIELD) { - fftwf_free(stars_unfiltered); - fftwf_free(stars_filtered); - fftwf_free(sfr_unfiltered); - fftwf_free(sfr_filtered); - } - - - LOG_SUPER_DEBUG("freed fftw boxes"); - - if(flag_options->USE_MASS_DEPENDENT_ZETA) { - if(flag_options->USE_MINI_HALOS){ - fftwf_free(log10_Mturnover_unfiltered); - fftwf_free(log10_Mturnover_filtered); - fftwf_free(log10_Mturnover_MINI_unfiltered); - fftwf_free(log10_Mturnover_MINI_filtered); - } - //fftwf_free(Mcrit_RE_grid); - //fftwf_free(Mcrit_LW_grid); - - } - - if (prev_redshift < 1){ - free(previous_ionize_box->z_re_box); - if (flag_options->USE_MASS_DEPENDENT_ZETA && flag_options->USE_MINI_HALOS){ - free(previous_ionize_box->Gamma12_box); - free(previous_ionize_box->dNrec_box); - free(previous_ionize_box->Fcoll); - free(previous_ionize_box->Fcoll_MINI); - } - } + LOG_DEBUG("global_xH = %e",global_xH); + free_fftw_grids(grid_struct); if(!flag_options->USE_TS_FLUCT && user_params->USE_INTERPOLATION_TABLES) { - freeSigmaMInterpTable(); + freeSigmaMInterpTable(); } - //These functions check for allocation + //This function checks for allocation so don't worry about double-freeing tables free_conditional_tables(); - free(overdense_int_boundexceeded_threaded); + free(radii_spec); LOG_DEBUG("finished!\n"); diff --git a/src/py21cmfast/src/PerturbField.c b/src/py21cmfast/src/PerturbField.c index 21314a96..850b0285 100644 --- a/src/py21cmfast/src/PerturbField.c +++ b/src/py21cmfast/src/PerturbField.c @@ -211,6 +211,7 @@ int ComputePerturbField( // *************** BEGIN INITIALIZATION ************************** // + LOG_DEBUG("Computing Perturbed Field at z=%.3f",redshift); // perform a very rudimentary check to see if we are underresolved and not using the linear approx if ((user_params->BOX_LEN > user_params->DIM) && !(global_params.EVOLVE_DENSITY_LINEARLY)){ LOG_WARNING("Resolution is likely too low for accurate evolved density fields\n \ @@ -247,8 +248,6 @@ int ComputePerturbField( // check if the linear evolution flag was set if (global_params.EVOLVE_DENSITY_LINEARLY){ - LOG_DEBUG("Linearly evolve density field"); - #pragma omp parallel shared(growth_factor,boxes,LOWRES_density_perturb,HIRES_density_perturb,dimension) private(i,j,k) num_threads(user_params->N_THREADS) { #pragma omp for @@ -268,7 +267,6 @@ int ComputePerturbField( } else { // Apply Zel'dovich/2LPT correction - LOG_DEBUG("Apply Zel'dovich"); #pragma omp parallel shared(LOWRES_density_perturb,HIRES_density_perturb,dimension) private(i,j,k) num_threads(user_params->N_THREADS) { @@ -317,7 +315,6 @@ int ComputePerturbField( // * ************************************************************************* * // // reference: reference: Scoccimarro R., 1998, MNRAS, 299, 1097-1118 Appendix D if(user_params->USE_2LPT){ - LOG_DEBUG("Apply 2LPT"); // allocate memory for the velocity boxes and read them in velocity_displacement_factor_2LPT = (displacement_factor_2LPT - init_displacement_factor_2LPT) / user_params->BOX_LEN; @@ -361,7 +358,6 @@ int ComputePerturbField( } // go through the high-res box, mapping the mass onto the low-res (updated) box - LOG_DEBUG("Perturb the density field"); #pragma omp parallel \ shared(init_growth_factor,boxes,f_pixel_factor,resampled_box,dimension) \ private(i,j,k,xi,xf,yi,yf,zi,zf,HII_i,HII_j,HII_k,d_x,d_y,d_z,t_x,t_y,t_z,xp1,yp1,zp1) \ @@ -528,7 +524,6 @@ int ComputePerturbField( } } free(resampled_box); - LOG_DEBUG("Finished perturbing the density field"); LOG_SUPER_DEBUG("density_perturb: "); if(user_params->PERTURB_ON_HIGH_RES){ @@ -581,14 +576,10 @@ int ComputePerturbField( } } } - LOG_DEBUG("Cleanup velocities for perturb"); } // Now, if I still have the high resolution density grid (HIRES_density_perturb) I need to downsample it to the low-resolution grid if(user_params->PERTURB_ON_HIGH_RES) { - - LOG_DEBUG("Downsample the high-res perturbed density"); - // Transform to Fourier space to sample (filter) the box dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, D_PARA, user_params->N_THREADS, HIRES_density_perturb); @@ -705,8 +696,6 @@ int ComputePerturbField( } // **** Convert to velocities ***** // - LOG_DEBUG("Generate velocity fields"); - float dDdt_over_D; dDdt_over_D = dDdt/growth_factor; @@ -767,6 +756,7 @@ int ComputePerturbField( fftwf_free(HIRES_density_perturb_saved); } fftwf_cleanup(); + LOG_DEBUG("Done."); } // End of Try{} Catch(status){ diff --git a/src/py21cmfast/src/SpinTemperatureBox.c b/src/py21cmfast/src/SpinTemperatureBox.c index 0affacb3..efd74765 100644 --- a/src/py21cmfast/src/SpinTemperatureBox.c +++ b/src/py21cmfast/src/SpinTemperatureBox.c @@ -97,9 +97,10 @@ int ComputeTsBox(float redshift, float prev_redshift, UserParams *user_params, C Try{ // This Try{} wraps the whole function. LOG_DEBUG("input values:"); LOG_DEBUG("redshift=%f, prev_redshift=%f perturbed_field_redshift=%f", redshift, prev_redshift, perturbed_field_redshift); - if (LOG_LEVEL >= DEBUG_LEVEL){ - writeAstroParams(flag_options, astro_params); - } + +#if LOG_LEVEL >= SUPER_DEBUG_LEVEL + writeAstroParams(flag_options, astro_params); +#endif // Makes the parameter structs visible to a variety of functions/macros // Do each time to avoid Python garbage collection issues @@ -313,7 +314,6 @@ void setup_z_edges(double zp){ double zpp, prev_zpp, prev_R; double dzpp_for_evolve; int R_ct; - LOG_DEBUG("Starting z edges"); R = L_FACTOR*user_params_global->BOX_LEN/(float)user_params_global->HII_DIM; R_factor = pow(global_params.R_XLy_MAX/R, 1/((float)global_params.NUM_FILTER_STEPS_FOR_Ts)); @@ -780,11 +780,6 @@ void fill_freqint_tables(double zp, double x_e_ave, double filling_factor_of_HI_ freq_int_lya_tbl_diff[x_e_ct-1][R_ct] = freq_int_lya_tbl[x_e_ct][R_ct] - freq_int_lya_tbl[x_e_ct-1][R_ct]; } } - LOG_ULTRA_DEBUG("%d of %d heat: %.3e %.3e %.3e ion: %.3e %.3e %.3e lya: %.3e %.3e %.3e lower %.3e" - ,R_ct,global_params.NUM_FILTER_STEPS_FOR_Ts - ,freq_int_heat_tbl[0][R_ct],freq_int_heat_tbl[x_int_NXHII/2][R_ct],freq_int_heat_tbl[x_int_NXHII-1][R_ct] - ,freq_int_ion_tbl[0][R_ct],freq_int_ion_tbl[x_int_NXHII/2][R_ct],freq_int_ion_tbl[x_int_NXHII-1][R_ct] - ,freq_int_lya_tbl[0][R_ct],freq_int_lya_tbl[x_int_NXHII/2][R_ct],freq_int_lya_tbl[x_int_NXHII-1][R_ct], lower_int_limit); } //separating the inverse diff loop to prevent a race on different R_ct (shouldn't matter) #pragma omp for @@ -908,8 +903,6 @@ int global_reion_properties(double zp, double x_e_ave, double *log10_Mcrit_LW_av fill_freqint_tables(zp,x_e_ave,*Q_HI,log10_Mcrit_LW_ave,0); } - LOG_DEBUG("Done."); - return sum_nion + sum_nion_mini > 1e-15 ? 0 : 1; //NO_LIGHT returned } @@ -997,7 +990,6 @@ struct spintemp_from_sfr_prefactors{ void set_zp_consts(double zp, struct spintemp_from_sfr_prefactors *consts){ //constant prefactors for the R_ct==0 part - LOG_DEBUG("Setting zp constants"); double luminosity_converstion_factor; double gamma_alpha; consts->growth_zp = dicke(zp); @@ -1258,8 +1250,6 @@ void ts_main(float redshift, float prev_redshift, UserParams *user_params, Cosmo for(R_ct=0;R_ct= SUPER_DEBUG_LEVEL +#if LOG_LEVEL >= SUPER_DEBUG_LEVEL if(box_ct==0 && !flag_options->USE_HALO_FIELD){ double integral_db; if(flag_options->USE_MASS_DEPENDENT_ZETA){ @@ -1555,7 +1545,7 @@ void ts_main(float redshift, float prev_redshift, UserParams *user_params, Cosmo if(flag_options->USE_LYA_HEATING) LOG_SUPER_DEBUG("ct %.2e | ij %.2e",dstarlya_cont_dt_box[box_ct],dstarlya_inj_dt_box[box_ct]); } - #endif +#endif } } } diff --git a/src/py21cmfast/src/Stochasticity.c b/src/py21cmfast/src/Stochasticity.c index a635c1df..2441f37d 100644 --- a/src/py21cmfast/src/Stochasticity.c +++ b/src/py21cmfast/src/Stochasticity.c @@ -17,6 +17,7 @@ #include "hmf.h" #include "cosmology.h" #include "InitialConditions.h" +#include "rng.h" #include "Stochasticity.h" //buffer size (per cell of arbitrary size) in the sampling function @@ -257,7 +258,7 @@ void set_prop_rng(gsl_rng *rng, bool from_catalog, double *interp, double * inpu int add_properties_cat(unsigned long long int seed, float redshift, HaloField *halos){ //set up the rng gsl_rng * rng_stoc[user_params_global->N_THREADS]; - seed_rng_threads(rng_stoc,seed); + seed_rng_threads_fast(rng_stoc,seed); LOG_DEBUG("computing rng for %llu halos",halos->n_halos); @@ -299,6 +300,48 @@ int stoc_halo_sample(struct HaloSamplingConstants *hs_constants, gsl_rng * rng, return 0; } +//This performs the number-limited sampling as above, modified to include a mass tolerance +// halo samples are entirely thrown out and re-drawn if they do not meet the tolerance +// This version is only currently used for internal testing against Nikolic et al. 2024, +int stoc_halo_sample_tol(struct HaloSamplingConstants *hs_constants, gsl_rng * rng, int *n_halo_out, float *M_out){ + double exp_N = hs_constants->expected_N; + double exp_M = hs_constants->expected_M; + int ii, nh=0; + int halo_count=0; + double mass_tot=0.; + double m_sample; + int attempts, n_failures; + + double tbl_arg = hs_constants->cond_val; + + //number limits that *can* fall in the tolerance + int min_N = 1; //assuming exp_M*0.9 < M_cond + int max_N = (int)(exp_M*1.1 / hs_constants->M_min); + + for(n_failures=0;n_failures<1000;n_failures++){ + do { + nh = gsl_ran_poisson(rng,exp_N); + } while((nh < min_N) || (nh > max_N)); + for(attempts=0;attempts<10;attempts++){ + halo_count = 0; + mass_tot = 0.; + for(ii=0;ii 0.9*exp_M) && (mass_tot < 1.1*exp_M)){ + *n_halo_out = halo_count; + return 0; + } + } + } + + //if we go through the whole loop raise an error + return TableEvaluationError; +} + double remove_random_halo(gsl_rng * rng, int n_halo, int *idx, double *M_prog, float *M_out){ double last_M_del; int random_idx; @@ -646,6 +689,7 @@ int stoc_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int //We always use Number-Limited sampling for grid-based cases if(user_params_global->SAMPLE_METHOD == 1 || !hs_constants->from_catalog){ err = stoc_halo_sample(hs_constants, rng, n_halo_out, M_out); + // err = stoc_halo_sample_tol(hs_constants, rng, n_halo_out, M_out); } else if(user_params_global->SAMPLE_METHOD == 0){ err = stoc_mass_sample(hs_constants, rng, n_halo_out, M_out); @@ -660,6 +704,10 @@ int stoc_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int LOG_ERROR("Invalid sampling method"); Throw(ValueError); } + if(err){ + LOG_ERROR("Error %d encountered in Halo sampling",err); + Throw(err); + } if(*n_halo_out > MAX_HALO_CELL){ LOG_ERROR("too many halos in condition, buffer overflow"); Throw(ValueError); @@ -953,7 +1001,7 @@ int stochastic_halofield(UserParams *user_params, CosmoParams *cosmo_params, //set up the rng gsl_rng * rng_stoc[user_params->N_THREADS]; - seed_rng_threads(rng_stoc,seed); + seed_rng_threads_fast(rng_stoc,seed); struct HaloSamplingConstants hs_constants; stoc_set_consts_z(&hs_constants,redshift,redshift_desc); @@ -1003,7 +1051,7 @@ int single_test_sample(UserParams *user_params, CosmoParams *cosmo_params, Astro //set up the rng gsl_rng * rng_stoc[user_params->N_THREADS]; - seed_rng_threads(rng_stoc,seed); + seed_rng_threads_fast(rng_stoc,seed); if(z_in > 0 && z_out <= z_in){ LOG_DEBUG("progenitor sampling must go back in time z_out=%.2f z_in=%.2f",z_out,z_in); diff --git a/src/py21cmfast/src/_inputparams_wrapper.h b/src/py21cmfast/src/_inputparams_wrapper.h index 67a2285c..8d3f2eb7 100644 --- a/src/py21cmfast/src/_inputparams_wrapper.h +++ b/src/py21cmfast/src/_inputparams_wrapper.h @@ -116,6 +116,7 @@ typedef struct FlagOptions{ bool CELL_RECOMB; int PHOTON_CONS_TYPE; bool USE_UPPER_STELLAR_TURNOVER; + bool HALO_SCALING_RELATIONS_MEDIAN; } FlagOptions; typedef struct GlobalParams{ diff --git a/src/py21cmfast/src/elec_interp.c b/src/py21cmfast/src/elec_interp.c index e7597ce1..a116988e 100755 --- a/src/py21cmfast/src/elec_interp.c +++ b/src/py21cmfast/src/elec_interp.c @@ -99,7 +99,7 @@ void initialize_interp_arrays() // Skip third line -- first have to get past the second line ending. skipline(input_file, 2); - LOG_SUPER_DEBUG("Reading %s", input_file_name); + LOG_ULTRA_DEBUG("Reading %s", input_file_name); for (i=0;i #include #include diff --git a/src/py21cmfast/src/interp_tables.c b/src/py21cmfast/src/interp_tables.c index a35b90ae..2af85914 100644 --- a/src/py21cmfast/src/interp_tables.c +++ b/src/py21cmfast/src/interp_tables.c @@ -296,18 +296,17 @@ void initialise_Nion_Conditional_spline(float z, float Mcrit_atom, float min_den sigma2 = EvaluateSigma(log(Mcond)); - if(prev){ - table_2d = &Nion_conditional_table_prev; - table_mini = &Nion_conditional_table_MINI_prev; - } - else{ - table_2d = &Nion_conditional_table2D; - table_mini = &Nion_conditional_table_MINI; - } - //If we use minihalos, both tables are 2D (delta,mturn) due to reionisaiton feedback //otherwise, the Nion table is 1D, since reionsaiton feedback is only active with minihalos if (minihalos){ + if(prev){ + table_2d = &Nion_conditional_table_prev; + table_mini = &Nion_conditional_table_MINI_prev; + } + else{ + table_2d = &Nion_conditional_table2D; + table_mini = &Nion_conditional_table_MINI; + } if(!table_2d->allocated) { allocate_RGTable2D_f(NDELTA,NMTURN,table_2d); } diff --git a/src/py21cmfast/src/recombinations.c b/src/py21cmfast/src/recombinations.c index ccf0658e..163d76b1 100755 --- a/src/py21cmfast/src/recombinations.c +++ b/src/py21cmfast/src/recombinations.c @@ -39,6 +39,8 @@ static double RR_table[RR_Z_NPTS][RR_lnGamma_NPTS], lnGamma_values[RR_lnGamma_NP static gsl_interp_accel *RR_acc[RR_Z_NPTS]; static gsl_spline *RR_spline[RR_Z_NPTS]; +static bool oob_warning_printed; + double recombination_rate(double z_eff, double gamma12_bg, double T4, int usecaseB); void free_MHR(); /* deallocates the gsl structures from init_MHR */ double Gamma_SS(double Gamma_bg, double Delta, double T_4, double z);//ionization rate w. self shielding @@ -72,8 +74,13 @@ double splined_recombination_rate(double z_eff, double gamma12_bg){ return 0; } else if (lnGamma >= (RR_lnGamma_min + RR_DEL_lnGamma * ( RR_lnGamma_NPTS - 1 )) ){ - LOG_WARNING("splined_recombination_rate: Gamma12 of %g is outside of interpolation array", gamma12_bg); lnGamma = RR_lnGamma_min + RR_DEL_lnGamma * ( RR_lnGamma_NPTS - 1 ) - FRACT_FLOAT_ERR; + if(!oob_warning_printed){ + LOG_WARNING("splined_recombination_rate: Gamma12 of %g is outside of interpolation array", gamma12_bg); + LOG_WARNING("gamma will be set to the maximum value %.3e", exp(lnGamma)); + LOG_WARNING("THIS WARNING IS PRINTED ONLY ONCE PER RUN"); + oob_warning_printed = true; + } } return gsl_spline_eval(RR_spline[z_ct], lnGamma, RR_acc[z_ct]); @@ -106,6 +113,7 @@ void init_MHR(){ gsl_spline_init(RR_spline[z_ct], lnGamma_values, RR_table[z_ct], RR_lnGamma_NPTS); } // go to next redshift + oob_warning_printed = false; return; } diff --git a/src/py21cmfast/src/rng.c b/src/py21cmfast/src/rng.c new file mode 100644 index 00000000..66bad999 --- /dev/null +++ b/src/py21cmfast/src/rng.c @@ -0,0 +1,161 @@ +#include +#include +#include +#include +#include +#include + +#include "cexcept.h" +#include "exceptions.h" +#include "logger.h" +#include "Constants.h" +#include "InputParameters.h" +#include "OutputStructs.h" + +#include "rng.h" + +//TODO: seed_rng_threads is slow, and isn't great to use every snapshot in Stochasticity.c or IonisationBox.c +// There are two remedies: +// - make it faster (by reducing the number of generated ints, I don't know why it's done that way) +// - have an rng struct which tracks its allocation and is only seeded once +// The second option creates difficulties i.e the same seed initialised at different points +// (for example by stopping/starting a run) can have different results +// seed_rng threads generates random seeds for each thread by choosing N_THREADS ints +// from an array of INT_MAX/16 (100 million?) ints. This array is shuffled, then each int is +// used as a seed with a looping list of 5 different gsl initialisers. +// In order to keep consistency with master I'm leaving it as is for now, but I really want to +// understand why it was done this way. +void seed_rng_threads(gsl_rng * rng_arr[], unsigned long long int seed){ + gsl_rng * rseed = gsl_rng_alloc(gsl_rng_mt19937); // An RNG for generating seeds for multithreading + + gsl_rng_set(rseed, seed); + + unsigned int seeds[user_params_global->N_THREADS]; + + // For multithreading, seeds for the RNGs are generated from an initial RNG (based on the input random_seed) and then shuffled (Author: Fred Davies) + int num_int = INT_MAX/16; + int i, thread_num; + unsigned int *many_ints = (unsigned int *)malloc((size_t)(num_int*sizeof(unsigned int))); // Some large number of possible integers + for (i=0; iN_THREADS, many_ints, num_int, sizeof(unsigned int)); // Populate the seeds array from the large list of integers + gsl_ran_shuffle(rseed, seeds, user_params_global->N_THREADS, sizeof(unsigned int)); // Shuffle the randomly selected integers + + int checker=0; + // seed the random number generators + for (thread_num = 0; thread_num < user_params_global->N_THREADS; thread_num++){ + switch (checker){ + case 0: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); + break; + case 1: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_gfsr4); + gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); + break; + case 2: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_cmrg); + gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); + break; + case 3: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_mrg); + gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); + break; + case 4: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_taus2); + gsl_rng_set(rng_arr[thread_num], seeds[thread_num]); + break; + } // end switch + + checker += 1; + + if(checker==5) { + checker = 0; + } + } + + free(many_ints); + gsl_rng_free(rseed); +} + +//Returns true if all elements of an array are unique, O(N^2) but N is small +bool array_unique(unsigned int array[], int array_size){ + int idx_0, idx_1; + for(idx_0=0;idx_0 1000 threads + int thread_num, checker; + + checker = 0; + // seed the random number generators + sample_n_unique_integers(user_params_global->N_THREADS, seed, st_arr); + + for (thread_num = 0; thread_num < user_params_global->N_THREADS; thread_num++){ + seed_thread = st_arr[thread_num]; + switch (checker){ + case 0: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(rng_arr[thread_num], seed_thread); + break; + case 1: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_gfsr4); + gsl_rng_set(rng_arr[thread_num], seed_thread); + break; + case 2: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_cmrg); + gsl_rng_set(rng_arr[thread_num], seed_thread); + break; + case 3: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_mrg); + gsl_rng_set(rng_arr[thread_num], seed_thread); + break; + case 4: + rng_arr[thread_num] = gsl_rng_alloc(gsl_rng_taus2); + gsl_rng_set(rng_arr[thread_num], seed_thread); + break; + } // end switch + checker += 1; + if(checker==5) { + checker = 0; + } + } +} + +void free_rng_threads(gsl_rng * rng_arr[]){ + int ii; + for(ii=0;iiN_THREADS;ii++){ + gsl_rng_free(rng_arr[ii]); + } +} diff --git a/src/py21cmfast/src/rng.h b/src/py21cmfast/src/rng.h new file mode 100644 index 00000000..0223ee74 --- /dev/null +++ b/src/py21cmfast/src/rng.h @@ -0,0 +1,10 @@ +#ifndef _21CMFRNG_H +#define _21CMFRNG_H + +#include + +void seed_rng_threads(gsl_rng * rng_arr[], unsigned long long int seed); +void seed_rng_threads_fast(gsl_rng * rng_arr[], unsigned long long int seed); +void free_rng_threads(gsl_rng * rng_arr[]); + +#endif diff --git a/src/py21cmfast/wrapper/cfuncs.py b/src/py21cmfast/wrapper/cfuncs.py index f15c481c..d8c36647 100644 --- a/src/py21cmfast/wrapper/cfuncs.py +++ b/src/py21cmfast/wrapper/cfuncs.py @@ -1203,13 +1203,14 @@ def halo_sample_test( if from_cat: z_prev = (1 + redshift) / global_params.ZPRIME_STEP_FACTOR - 1 + buffer_size = int(3e7) # about 500MB total 2e7 * 4 (float) * 4 (mass + 3crd) nhalo_out = np.zeros(1).astype("i4") N_out = np.zeros(n_cond).astype("i4") M_out = np.zeros(n_cond).astype("f8") exp_M = np.zeros(n_cond).astype("f8") exp_N = np.zeros(n_cond).astype("f8") - halomass_out = np.zeros(int(1e7)).astype("f4") - halocrd_out = np.zeros(int(3e7)).astype("i4") + halomass_out = np.zeros(buffer_size).astype("f4") + halocrd_out = np.zeros(int(3 * buffer_size)).astype("i4") lib.single_test_sample( user_params.cstruct, diff --git a/src/py21cmfast/wrapper/globals.py b/src/py21cmfast/wrapper/globals.py index 786940a2..5edb930c 100644 --- a/src/py21cmfast/wrapper/globals.py +++ b/src/py21cmfast/wrapper/globals.py @@ -322,5 +322,12 @@ def use(self, **kwargs): for k, v in prev.items(): setattr(self, k, v) + def validate(self): + """Validate the global parameters to ensure the passed values do not clash.""" + if getattr(self._cobj, "FIND_BUBBLE_ALGORITHM") not in [1, 2]: + raise ValueError( + "FIND_BUBBLE_ALGORITHM MUST BE 1 (entire sphere) or 2 (central pixel)" + ) + global_params = GlobalParams(lib.global_params, ffi) diff --git a/src/py21cmfast/wrapper/inputs.py b/src/py21cmfast/wrapper/inputs.py index f3740427..b4ae5006 100644 --- a/src/py21cmfast/wrapper/inputs.py +++ b/src/py21cmfast/wrapper/inputs.py @@ -454,6 +454,10 @@ class FlagOptions(InputStruct): 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. + HALO_SCALING_RELATIONS_MEDIAN: bool, optional + If True, halo scaling relation parameters (F_STAR10,t_STAR etc...) define the median of their conditional distributions + If False, they describe the mean. + This becomes important when using non-symmetric dristributions such as the log-normal """ _photoncons_models = [ @@ -485,6 +489,7 @@ class FlagOptions(InputStruct): ) USE_UPPER_STELLAR_TURNOVER = field(default=True, converter=bool) M_MIN_in_Mass = field(default=True, converter=bool) + HALO_SCALING_RELATIONS_MEDIAN = field(default=False, converter=bool) @M_MIN_in_Mass.validator def _M_MIN_in_Mass_vld(self, att, val): diff --git a/src/py21cmfast/wrapper/photoncons.py b/src/py21cmfast/wrapper/photoncons.py index bb0eed09..23530fae 100644 --- a/src/py21cmfast/wrapper/photoncons.py +++ b/src/py21cmfast/wrapper/photoncons.py @@ -462,8 +462,7 @@ def photoncons_alpha(cosmo_params, user_params, astro_params, flag_options): ) for i, a in enumerate(alpha_arr): # alter astro params with new alpha - astro_params_photoncons = deepcopy(astro_params) - astro_params_photoncons.ALPHA_ESC = a + astro_params_photoncons = astro_params.clone(ALPHA_ESC=a) # find the analytic curve wth that alpha # TODO: Theres a small memory leak here since global arrays are allocated (for some reason) diff --git a/src/py21cmfast/wrapper/structs.py b/src/py21cmfast/wrapper/structs.py index 28ebf74c..58504835 100644 --- a/src/py21cmfast/wrapper/structs.py +++ b/src/py21cmfast/wrapper/structs.py @@ -1082,26 +1082,33 @@ def ensure_input_computed(self, input_box, load=False) -> bool: def summarize(self, indent=0) -> str: """Generate a string summary of the struct.""" indent = indent * " " - out = f"\n{indent}{self.__class__.__name__}\n" - out += "".join( - f"{indent} {fieldname:>15}: {getattr(self, fieldname, 'non-existent')}\n" - for fieldname in self.struct.primitive_fields + # print array type and column headings + out = ( + f"\n{indent}{self.__class__.__name__:>25} " + + " 1st: End: Min: Max: Mean: \n" ) + # print array extrema and means for fieldname, state in self._array_state.items(): if not state.initialized: - out += f"{indent} {fieldname:>15}: uninitialized\n" + out += f"{indent} {fieldname:>25}: uninitialized\n" elif not state.computed: - out += f"{indent} {fieldname:>15}: initialized\n" + out += f"{indent} {fieldname:>25}: initialized\n" elif not state.computed_in_mem: - out += f"{indent} {fieldname:>15}: computed on disk\n" + out += f"{indent} {fieldname:>25}: computed on disk\n" else: x = getattr(self, fieldname).flatten() if len(x) > 0: - out += f"{indent} {fieldname:>15}: {x[0]:1.4e}, {x[-1]:1.4e}, {x.min():1.4e}, {x.max():1.4e}, {np.mean(x):1.4e}\n" + out += f"{indent} {fieldname:>25}: {x[0]:11.4e}, {x[-1]:11.4e}, {x.min():11.4e}, {x.max():11.4e}, {np.mean(x):11.4e}\n" else: - out += f"{indent} {fieldname:>15}: size zero\n" + out += f"{indent} {fieldname:>25}: size zero\n" + + # print primitive fields + out += "".join( + f"{indent} {fieldname:>25}: {getattr(self, fieldname, 'non-existent')}\n" + for fieldname in self.struct.primitive_fields + ) return out @@ -1242,6 +1249,7 @@ def __setattr__(self, name, value): with contextlib.suppress(AttributeError): setattr(self._cobj, name, value) object.__setattr__(self, name, value) + self.validate() def items(self): """Yield (name, value) pairs for each element of the struct.""" @@ -1277,3 +1285,7 @@ def filtered_repr(self, filter_params): if k not in filter_params ) ) + ")" + + def validate(self): + """Validate members of this structure against each other.""" + pass diff --git a/tests/produce_integration_test_data.py b/tests/produce_integration_test_data.py index fdf27b69..f2ee4a2b 100644 --- a/tests/produce_integration_test_data.py +++ b/tests/produce_integration_test_data.py @@ -49,7 +49,7 @@ "DIM": 150, "BOX_LEN": 100, "NO_RNG": True, - "SAMPLER_MIN_MASS": 5e8, + "SAMPLER_MIN_MASS": 1e9, } DEFAULT_FLAG_OPTIONS = { diff --git a/tests/test_c_interpolation_tables.py b/tests/test_c_interpolation_tables.py index 3ba6dccb..469a087f 100644 --- a/tests/test_c_interpolation_tables.py +++ b/tests/test_c_interpolation_tables.py @@ -1,6 +1,5 @@ import pytest -import attrs import matplotlib as mpl import numpy as np from astropy import constants as c @@ -1035,6 +1034,10 @@ def make_integral_comparison_plot(x1, x2, integral_list, integral_list_second, p axs[1, 1].set_xlabel("delta") axs[1, 0].set_ylabel("Integral") axs[0, 0].set_ylabel("Integral") + axs[0, 0].grid() + axs[0, 1].grid() + axs[1, 0].grid() + axs[1, 1].grid() # copied and expanded from test_integration_features.py diff --git a/tests/test_halo_sampler.py b/tests/test_halo_sampler.py index aaa13bd3..f060f3eb 100644 --- a/tests/test_halo_sampler.py +++ b/tests/test_halo_sampler.py @@ -34,14 +34,14 @@ def test_sampler(name, cond, from_cat, plt): redshift, kwargs = cint.OPTIONS_HMF[name] redshift = 8 opts = prd.get_all_options(redshift, **kwargs) - up = opts["user_params"] + up = opts["user_params"].clone(SAMPLER_MIN_MASS=2e8) cp = opts["cosmo_params"] ap = opts["astro_params"] fo = opts["flag_options"] from_cat = "cat" in from_cat - n_cond = 10000 + n_cond = 15000 if from_cat: mass = 10 ** options_log10mass[cond] cond = mass @@ -70,11 +70,12 @@ def test_sampler(name, cond, from_cat, plt): redshift=redshift, from_cat=from_cat, cond_array=np.full(n_cond, cond), + seed=987, ) # set up histogram l10min = np.log10(up.SAMPLER_MIN_MASS) - l10max = np.log10(1e14) + l10max = np.log10(mass * 1.01) edges = np.logspace(l10min, l10max, num=int(10 * (l10max - l10min))) bin_minima = edges[:-1] bin_maxima = edges[1:] @@ -133,7 +134,7 @@ def test_sampler(name, cond, from_cat, plt): # The histograms get inaccurate when the volume is too small # so only compare when we expect real halos if sample_dict["expected_progenitor_mass"][0] > up.SAMPLER_MIN_MASS: - sel_compare_bins = edges[:-1] < (0.9 * mass) + sel_compare_bins = edges[1:] < (0.9 * mass) print_failure_stats( mf_out[sel_compare_bins], @@ -243,7 +244,6 @@ def plot_sampler_comparison( # mass function axis axs[0].set_title(title) axs[0].set_ylim([1e-6, 1e2]) - axs[0].set_xlim([bin_edges[0], bin_edges[np.argmax(exp_mf == 0)]]) axs[0].set_yscale("log") axs[0].set_ylabel("dn/dlnM") diff --git a/tests/test_segfaults.py b/tests/test_segfaults.py new file mode 100644 index 00000000..f7b20af3 --- /dev/null +++ b/tests/test_segfaults.py @@ -0,0 +1,138 @@ +""" +This file contains tests which run a lightcone under various flag options +to test the C backend for segfaults. +These will not test the outputs of the run past the fact that they are finite, +just that the run completes without error +""" + +import pytest + +import numpy as np + +import py21cmfast as p21c + +from . import produce_integration_test_data as prd + +DEFAULT_USER_PARAMS_CTEST = { + "HII_DIM": 32, + "DIM": 128, + "BOX_LEN": 64, + "NO_RNG": True, + "SAMPLER_MIN_MASS": 1e9, +} + +OPTIONS_CTEST = { + "defaults": [20, {"USE_MASS_DEPENDENT_ZETA": True}], + "no-mdz": [20, {}], + "mini": [ + 20, + { + "USE_MINI_HALOS": True, + "INHOMO_RECO": True, + "USE_TS_FLUCT": True, + "USE_MASS_DEPENDENT_ZETA": True, + }, + ], + "ts": [20, {"USE_TS_FLUCT": True, "USE_MASS_DEPENDENT_ZETA": True}], + "ts_nomdz": [20, {"USE_TS_FLUCT": True}], + "inhomo": [20, {"INHOMO_RECO": True, "USE_MASS_DEPENDENT_ZETA": True}], + "inhomo_ts": [ + 20, + {"INHOMO_RECO": True, "USE_TS_FLUCT": True, "USE_MASS_DEPENDENT_ZETA": True}, + ], + "sampler": [ + 20, + { + "USE_HALO_FIELD": True, + "HALO_STOCHASTICITY": True, + "USE_MASS_DEPENDENT_ZETA": True, + }, + ], + "sampler_mini": [ + 20, + { + "USE_HALO_FIELD": True, + "HALO_STOCHASTICITY": True, + "USE_MINI_HALOS": True, + "USE_TS_FLUCT": True, + "INHOMO_RECO": True, + "USE_MASS_DEPENDENT_ZETA": True, + }, + ], + "sampler_ts": [ + 20, + { + "USE_HALO_FIELD": True, + "HALO_STOCHASTICITY": True, + "USE_TS_FLUCT": True, + "USE_MASS_DEPENDENT_ZETA": True, + }, + ], + "sampler_ir": [ + 20, + { + "USE_HALO_FIELD": True, + "HALO_STOCHASTICITY": True, + "INHOMO_RECO": True, + "USE_MASS_DEPENDENT_ZETA": True, + }, + ], + "sampler_ts_ir": [ + 20, + { + "USE_HALO_FIELD": True, + "HALO_STOCHASTICITY": True, + "USE_TS_FLUCT": True, + "INHOMO_RECO": True, + "USE_MASS_DEPENDENT_ZETA": True, + }, + ], + "photoncons-z": [ + 20, + {"PHOTON_CONS_TYPE": "z-photoncons", "USE_MASS_DEPENDENT_ZETA": True}, + ], + "photoncons-a": [ + 20, + {"PHOTON_CONS_TYPE": "alpha-photoncons", "USE_MASS_DEPENDENT_ZETA": True}, + ], + "photoncons-f": [ + 20, + {"PHOTON_CONS_TYPE": "f-photoncons", "USE_MASS_DEPENDENT_ZETA": True}, + ], + "minimize_mem": [ + 20, + { + "USE_TS_FLUCT": True, + "INHOMO_RECO": True, + "MINIMIZE_MEMORY": True, + "USE_MASS_DEPENDENT_ZETA": True, + }, + ], +} + + +@pytest.mark.parametrize("name", list(OPTIONS_CTEST.keys())) +def test_lc_runs(name, max_redshift): + redshift, kwargs = OPTIONS_CTEST[name] + options = prd.get_all_options(redshift, **kwargs) + + options["user_params"] = p21c.UserParams.new(DEFAULT_USER_PARAMS_CTEST) + + lcn = p21c.RectilinearLightconer.with_equal_cdist_slices( + min_redshift=redshift, + max_redshift=max_redshift, + quantities=[ + "brightness_temp", + ], + resolution=options["user_params"].cell_size, + ) + + with p21c.config.use(ignore_R_BUBBLE_MAX_error=True): + _, _, _, lightcone = p21c.exhaust_lightcone( + lightconer=lcn, + write=False, + **options, + ) + + assert isinstance(lightcone, p21c.LightCone) + assert np.all(np.isfinite(lightcone.brightness_temp))