Ionbox refactor

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){

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)],

src/py21cmfast/src/InitialConditions.c

@@ -22,17 +22,26 @@
 
 #include "InitialConditions.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){
-    // 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 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; i<num_int; i++) {
@@ -42,9 +51,7 @@ void seed_rng_threads(gsl_rng * rng_arr[], unsigned long long int seed){
     gsl_ran_choose(rseed, seeds, user_params_global->N_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;
+    int checker=0;
     // seed the random number generators
     for (thread_num = 0; thread_num < user_params_global->N_THREADS; thread_num++){
         switch (checker){
@@ -77,8 +84,56 @@ void seed_rng_threads(gsl_rng * rng_arr[], unsigned long long int seed){
         }
     }
 
-    gsl_rng_free(rseed);
     free(many_ints);
+    gsl_rng_free(rseed);
+}
+
+//I'm putting a faster version of the RNG seeding here to use when we need one seeding per snapshot.
+//  This is exactly the same as above except instead of generating a huge list of ints to select from,
+//  it generates seeds using gsl_rng_uniform_int(), which means there *can* be a repetition.
+//  This should be very rare, and even when it occurs, it will usually use a different initialiser,
+//  as well as be applied to different cells/halos, so I can't forsee any issues.
+//Just in case I'm not using it for the initial conditions, which is okay since the slower version is only used once
+void seed_rng_threads_fast(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 seed_thread;
+    int num_int = INT_MAX/16;
+    int thread_num, checker;
+
+    checker = 0;
+    // seed the random number generators
+    for (thread_num = 0; thread_num < user_params_global->N_THREADS; thread_num++){
+        seed_thread = gsl_rng_uniform_int(rseed,num_int);
+        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;
+        }
+    }
+    gsl_rng_free(rseed);
 }
 
 void free_rng_threads(gsl_rng * rng_arr[]){

src/py21cmfast/src/InitialConditions.h

@@ -11,7 +11,10 @@ int ComputeInitialConditions(
     CosmoParams *cosmo_params, InitialConditions *boxes
 );
 
+//TODO: these seeding functions should probably be somewhere else
+//  Possibly make an rng.c/h
 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