Merge branch '21cmfast:master' into Radio_Excess

.gitignore

@@ -53,3 +53,4 @@ r2c*
 plots/
 .vscode/
 src/21cmfast/_version.py
+src/py21cmfast/_version.py

.pre-commit-config.yaml

@@ -2,7 +2,7 @@ exclude: '(^docs/conf.py|^src/py21cmfast/_data/|changethelog.py)'
 
 repos:
     - repo: https://github.com/pre-commit/pre-commit-hooks
-      rev: v4.2.0
+      rev: v4.3.0
       hooks:
           - id: trailing-whitespace
           - id: check-added-large-files
@@ -17,7 +17,7 @@ repos:
           - id: mixed-line-ending
             args: ['--fix=no']
     - repo: https://github.com/pycqa/flake8
-      rev: 4.0.1  # pick a git hash / tag to point to
+      rev: 5.0.4  # pick a git hash / tag to point to
       hooks:
         - id: flake8
           additional_dependencies:
@@ -34,12 +34,12 @@ repos:
             - flake8-print
 
     - repo: https://github.com/psf/black
-      rev: 22.3.0
+      rev: 22.10.0
       hooks:
           - id: black
 
-    - repo: https://github.com/pre-commit/mirrors-isort
-      rev: v5.10.1
+    - repo: https://github.com/PyCQA/isort
+      rev: 5.10.1
       hooks:
       - id: isort
 
@@ -51,7 +51,7 @@ repos:
       - id: rst-inline-touching-normal
 
     - repo: https://github.com/asottile/pyupgrade
-      rev: v2.32.1
+      rev: v3.1.0
       hooks:
       - id: pyupgrade
         args: [--py36-plus]

CHANGELOG.rst

@@ -21,6 +21,10 @@ Added
 * New ``validate_all_inputs`` function that cross-references the four main input structs
   and ensures all the parameters make sense together. Mostly for internal use.
 * Ability to save/read directly from an open HDF5 File (#170)
+* An implementation of cloud-in-cell to more accurately redistribute the perturbed mass
+  across all neighbouring cells instead of the previous nearest cell approach
+* Changed PhotonConsEndCalibz from z = 5 -> z = 3.5 to handle later reionisation
+  scenarios in line with current observations (#305)
 
 v3.2.1 [13 Sep 2022]
 ----------------------

docs/design.rst

@@ -62,18 +62,12 @@ changing the parameters, and can use the outputs in whatever way the user desire
 To maintain continuity with previous versions, a CLI interface is provided (see below)
 that acts in a similar fashion to previous versions.
 
-When ``21cmFAST`` is installed, it automatically creates a configuration directory in
-the user's home: ``~/.21cmfast``. This houses a number of important configuration
-options; usually default values of parameters. At this stage, the location of this
-directory is not itself configurable. The config directory contains example
-configuration files for the CLI interface (see below), which can also be copied anywhere
-on disk and modified. Importantly, the ``config.yml`` file in this directory specifies
-some of the more important options for how ``21cmFAST`` behaves by default.
-One such option is ``boxdir``, which specifies the directory in which ``21cmFAST`` will
-cache results (see below for details). Finally, the config directory houses several data
-tables which are used to accelerate several calculations. In principle
-these files are over-writeable, but they should only be touched if one knows very well
-what they are doing.
+High-level configuration of ``21cmFAST`` can be set using the ``py21cmfast.config``
+object. It is essentially a dictionary with its key/values the parameters. To make any
+changes in the object permanent, use the ``py21cmfast.config.write()`` method.
+One global configuration option is ``direc``, which specifies the directory in which
+``21cmFAST`` will cache results by default (this can be overriden directly in any
+function, see below for details).
 
 Finally, ``21cmFAST`` contains a more robust cataloguing/caching method. Instead of
 saving data with a selection of the dependent parameters written into the filename --
@@ -102,8 +96,7 @@ To get help on any subcommand, simply use::
     $ 21cmfast <subcommand> --help
 
 Any subcommand which runs some aspect of ``21cmFAST`` will have a ``--config`` option,
-which specifies a configuration file (by default this is
-``~/.21cmfast/runconfig_example.yml``). This config file specifies the parameters of the
+which specifies a configuration file. This config file specifies the parameters of the
 run. Furthermore, any particular parameter that can be specified in the config file can
 be alternatively specified on the command line by appending the command with the
 parameter name, eg.::

src/py21cmfast/src/GenerateICs.c

@@ -571,7 +571,7 @@ int ComputeInitialConditions(
 
                 // generate the phi_1 boxes in Fourier transform
 #pragma omp parallel shared(HIRES_box,phi_1,i,j) private(n_x,n_y,n_z,k_x,k_y,k_z,k_sq,k) num_threads(user_params->N_THREADS)
-                {
+            {
 #pragma omp for
                 for (n_x=0; n_x<user_params->DIM; n_x++){
                     if (n_x>MIDDLE)

src/py21cmfast/src/Globals.h

@@ -116,7 +116,7 @@ extern struct GlobalParams global_params = {
     .PhotonConsStart = 0.995,
     .PhotonConsEnd = 0.3,
     .PhotonConsAsymptoteTo = 0.01,
-    .PhotonConsEndCalibz = 5.0,
+    .PhotonConsEndCalibz = 3.5,
 
     .HEAT_FILTER = 0,
     .CLUMPING_FACTOR = 2.,

src/py21cmfast/src/PerturbField.c

@@ -24,12 +24,17 @@ int ComputePerturbField(
     fftwf_complex *HIRES_density_perturb, *HIRES_density_perturb_saved;
     fftwf_complex *LOWRES_density_perturb, *LOWRES_density_perturb_saved;
 
-    float growth_factor, displacement_factor_2LPT, init_growth_factor, init_displacement_factor_2LPT, xf, yf, zf;
+    float growth_factor, displacement_factor_2LPT, init_growth_factor, init_displacement_factor_2LPT;
+    double xf, yf, zf;
     float mass_factor, dDdt, f_pixel_factor, velocity_displacement_factor, velocity_displacement_factor_2LPT;
     unsigned long long ct, HII_i, HII_j, HII_k;
-    int i,j,k, xi, yi, zi, dimension, switch_mid;
+    int i,j,k,xi, yi, zi, dimension, switch_mid;
     double ave_delta, new_ave_delta;
 
+    // Variables to perform cloud in cell re-distribution of mass for the perturbed field
+    int xp1,yp1,zp1;
+    float d_x,d_y,d_z,t_x,t_y,t_z;
+
     // Function for deciding the dimensions of loops when we could
     // use either the low or high resolution grids.
     switch(user_params->PERTURB_ON_HIGH_RES) {
@@ -196,7 +201,7 @@ 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) num_threads(user_params->N_THREADS)
+                        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) num_threads(user_params->N_THREADS)
         {
 #pragma omp for
             for (i=0; i<user_params->DIM;i++){
@@ -237,15 +242,14 @@ int ComputePerturbField(
                                 zf -= (boxes->lowres_vz_2LPT)[HII_R_INDEX(HII_i,HII_j,HII_k)];
                             }
                         }
-
-                        xf *= (float)(dimension);
-                        yf *= (float)(dimension);
-                        zf *= (float)(dimension);
-                        while (xf >= (float)(dimension)){ xf -= (dimension);}
+                        xf *= (double)(dimension);
+                        yf *= (double)(dimension);
+                        zf *= (double)(dimension);
+                        while (xf >= (double)(dimension)){ xf -= (dimension);}
                         while (xf < 0){ xf += (dimension);}
-                        while (yf >= (float)(dimension)){ yf -= (dimension);}
+                        while (yf >= (double)(dimension)){ yf -= (dimension);}
                         while (yf < 0){ yf += (dimension);}
-                        while (zf >= (float)(dimension)){ zf -= (dimension);}
+                        while (zf >= (double)(dimension)){ zf -= (dimension);}
                         while (zf < 0){ zf += (dimension);}
                         xi = xf;
                         yi = yf;
@@ -257,13 +261,79 @@ int ComputePerturbField(
                         if (zi >= (dimension)){ zi -= (dimension);}
                         if (zi < 0) {zi += (dimension);}
 
+                        // Determine the fraction of the perturbed cell which overlaps with the 8 nearest grid cells,
+                        // based on the grid cell which contains the centre of the perturbed cell
+                        d_x = fabs(xf - (double)(xi+0.5));
+                        d_y = fabs(yf - (double)(yi+0.5));
+                        d_z = fabs(zf - (double)(zi+0.5));
+                        if(xf < (double)(xi+0.5)) {
+                            // If perturbed cell centre is less than the mid-point then update fraction
+                            // of mass in the cell and determine the cell centre of neighbour to be the
+                            // lowest grid point index
+                            d_x = 1. - d_x;
+                            xi -= 1;
+                            if (xi < 0) {xi += (dimension);} // Only this critera is possible as iterate back by one (we cannot exceed DIM)
+                        }
+                        if(yf < (double)(yi+0.5)) {
+                            d_y = 1. - d_y;
+                            yi -= 1;
+                            if (yi < 0) {yi += (dimension);}
+                        }
+                        if(zf < (double)(zi+0.5)) {
+                            d_z = 1. - d_z;
+                            zi -= 1;
+                            if (zi < 0) {zi += (dimension);}
+                        }
+                        t_x = 1. - d_x;
+                        t_y = 1. - d_y;
+                        t_z = 1. - d_z;
+
+                        // Determine the grid coordinates of the 8 neighbouring cells
+                        // Takes into account the offset based on cell centre determined above
+                        xp1 = xi + 1;
+                        if(xp1 >= dimension) { xp1 -= (dimension);}
+                        yp1 = yi + 1;
+                        if(yp1 >= dimension) { yp1 -= (dimension);}
+                        zp1 = zi + 1;
+                        if(zp1 >= dimension) { zp1 -= (dimension);}
+
                         if(user_params->PERTURB_ON_HIGH_RES) {
+                            // Redistribute the mass over the 8 neighbouring cells according to cloud in cell
+#pragma omp atomic
+                                resampled_box[R_INDEX(xi,yi,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*t_y*t_z);
+#pragma omp atomic
+                                resampled_box[R_INDEX(xp1,yi,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*t_y*t_z);
+#pragma omp atomic
+                                resampled_box[R_INDEX(xi,yp1,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*d_y*t_z);
+#pragma omp atomic
+                                resampled_box[R_INDEX(xp1,yp1,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*d_y*t_z);
 #pragma omp atomic
-                            resampled_box[R_INDEX(xi,yi,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)]);
+                                resampled_box[R_INDEX(xi,yi,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*t_y*d_z);
+#pragma omp atomic
+                                resampled_box[R_INDEX(xp1,yi,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*t_y*d_z);
+#pragma omp atomic
+                                resampled_box[R_INDEX(xi,yp1,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*d_y*d_z);
+#pragma omp atomic
+                                resampled_box[R_INDEX(xp1,yp1,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*d_y*d_z);
                         }
                         else {
+                            // Redistribute the mass over the 8 neighbouring cells according to cloud in cell
+#pragma omp atomic
+                                resampled_box[HII_R_INDEX(xi,yi,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*t_y*t_z);
+#pragma omp atomic
+                                resampled_box[HII_R_INDEX(xp1,yi,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*t_y*t_z);
+#pragma omp atomic
+                                resampled_box[HII_R_INDEX(xi,yp1,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*d_y*t_z);
+#pragma omp atomic
+                                resampled_box[HII_R_INDEX(xp1,yp1,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*d_y*t_z);
+#pragma omp atomic
+                                resampled_box[HII_R_INDEX(xi,yi,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*t_y*d_z);
+#pragma omp atomic
+                                resampled_box[HII_R_INDEX(xp1,yi,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*t_y*d_z);
+#pragma omp atomic
+                                resampled_box[HII_R_INDEX(xi,yp1,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(t_x*d_y*d_z);
 #pragma omp atomic
-                            resampled_box[HII_R_INDEX(xi,yi,zi)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)]);
+                                resampled_box[HII_R_INDEX(xp1,yp1,zp1)] += (double)(1. + init_growth_factor*(boxes->hires_density)[R_INDEX(i,j,k)])*(d_x*d_y*d_z);
                         }
                     }
                 }

tests/test_wrapper.py

@@ -405,9 +405,9 @@ def test_coeval_vs_low_level(ic):
         write=False,
     )
 
-    assert np.allclose(coeval.Tk_box, st.Tk_box)
-    assert np.allclose(coeval.Ts_box, st.Ts_box)
-    assert np.allclose(coeval.x_e_box, st.x_e_box)
+    assert np.allclose(coeval.Tk_box, st.Tk_box, rtol=1e-4)
+    assert np.allclose(coeval.Ts_box, st.Ts_box, rtol=1e-4)
+    assert np.allclose(coeval.x_e_box, st.x_e_box, rtol=1e-4)
 
 
 def test_using_cached_halo_field(ic, test_direc):