Shift the first 128 k-bins in the power spectrum output to take into account...

Shift the first 128 k-bins in the power spectrum output to take into account the weights of the different contributing modes.

doc/RTD/source/ParameterFiles/parameter_description.rst

@@ -1280,6 +1280,7 @@ The options are:
  * The number of grid foldings to use: ``num_folds``.
  * The factor by which to fold at each iteration: ``fold_factor`` (default: 4)
  * The order of the window function: ``window_order`` (default: 3)
+ * Whether or not to correct the placement of the centre of the k-bins for small k values: ``shift_centre_small_k_bins`` (default: 1)
 
 The window order sets the way the particle properties get assigned to the mesh.
 Order 1 corresponds to the nearest-grid-point (NGP), order 2 to cloud-in-cell
@@ -1319,6 +1320,20 @@ are mutually exclusive. The particles selected in each half are different in
 each output. Note that neutrino PS can only be computed when neutrinos are
 simulated using particles.
 
+SWIFT uses bins of integer :math:`k`, with bins :math:`[0.5,1.5]`, :math:`[1.5,2.5]` etc.  The
+representative :math:`k` values used to be assigned to the bin centres (so k=1, 2, etc), which are
+then transformed to physical :math:`k` by a factor :math:`kL/(2*pi)`. For the first few bins, only
+few modes contribute to each bin. It is then advantageous to move the "centre" of the bin to the
+actual location correponding to the mean of the contributing modes. The :math:`k` label of the bin
+is thus shifted by a small amount. The way to calculate these shifts is to consider a 3D cube of
+:math:`(kx,ky,kz)` cells and check which cells fall inside a spherical shell with boundaries
+:math:`(i+0.5,i+1.5)`, then calculate the average :math:`k=\sqrt{kx^2+ky^2+kz^2}`. So for
+:math:`i=0` there cells :math:`k=1` and 12 cells :math:`k=\sqrt(2)`, so the weighted k becomes
+:math:`(6 * 1 + 12 * \sqrt(2)) / 18 = 1.2761424`. Note that only the first 7 (22) bins require a
+correction larger than 1 (0.1) percent. We apply a correction to the first 128 terms. This
+correction is activated when ``shift_centre_small_k_bins`` is switched on (that's the default
+behaviour).
+
 An example of a valid power-spectrum section of the parameter file looks like:
 
 .. code:: YAML

examples/parameter_example.yml

@@ -932,6 +932,7 @@ PowerSpectrum:
   num_folds:         6                    # Number of foldings (1 means no foldings), determines the max k
   fold_factor:       4                    # (Optional) factor by which to reduce the box along each side each folding (default: 4)
   window_order:      3                    # (Optional) order of the mass assignment scheme (default: 3, TSC)
+  shift_centre_small_k_bins: 1            # (Optional) Correct the centre of the bins with a small k to account for the small number of modes entering the bin.
   output_list_on:    0                    # (Optional) Enable the output list
   output_list:       ./output_list_ps.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
   requested_spectra: ["matter-matter","cdm-cdm","starBH-starBH","gas-matter","pressure-pressure","matter-pressure", "neutrino0-neutrino1"] # Array of strings indicating which components should be correlated for power spectra

src/power_spectrum.c

@@ -51,6 +51,37 @@
 #define power_data_default_fold_factor 4
 #define power_data_default_window_order 3
 
+/* The way to calculate these shifts is to consider a 3D cube of (kx,ky,kz)
+ * cells and check which cells fall inside a spherical shell with boundaries
+ * (i+0.5,i+1.5), then calculate the average k=sqrt(kx^2+ky^2+kz^2). So for i=0
+ * you'd find 6 cells k=1 and 12 cells k=sqrt(2), so the weighted k becomes
+ * (6 * 1 + 12 * sqrt(2)) / 18 = 1.2761424 – etc.
+ * Note that beyond the 7th term, the correction is < 1%. */
+#define number_of_corrected_bins 128
+static const float correction_shift_k_values[number_of_corrected_bins] = {
+    1.2761424f, 1.1154015f, 1.0447197f, 1.0151449f, 1.0195166f, 1.0203214f,
+    1.0102490f, 1.0031348f, 1.0063766f, 1.0093355f, 1.0055681f, 1.0024279f,
+    1.0034435f, 1.0038386f, 1.0011069f, 1.0002888f, 1.0018693f, 1.0029172f,
+    1.0019128f, 1.0009282f, 1.0015312f, 1.0016361f, 1.0009436f, 1.0003777f,
+    1.0005931f, 1.0010948f, 1.0010581f, 1.0009779f, 1.0010282f, 1.0008224f,
+    1.0006637f, 1.0004002f, 1.0002419f, 1.0005172f, 1.0005523f, 1.0004342f,
+    1.0005183f, 1.0005357f, 1.0003162f, 1.0001836f, 1.0003737f, 1.0004792f,
+    1.0004169f, 1.0003660f, 1.0004468f, 1.0004218f, 1.0001436f, 1.0000479f,
+    1.0002012f, 1.0003710f, 1.0003234f, 1.0002661f, 1.0003446f, 1.0003313f,
+    1.0001844f, 1.0000630f, 1.0001714f, 1.0002382f, 1.0001507f, 1.0001663f,
+    1.0002199f, 1.0002403f, 1.0000911f, 0.9999714f, 1.0001136f, 1.0001907f,
+    1.0001917f, 1.0001684f, 1.0001875f, 1.0002158f, 1.0000941f, 1.0000646f,
+    1.0000930f, 1.0001497f, 1.0001589f, 1.0001215f, 1.0001563f, 1.0001254f,
+    1.0000557f, 1.0000220f, 1.0000517f, 1.0001039f, 1.0001185f, 1.0000778f,
+    1.0000848f, 1.0001415f, 1.0001108f, 1.0000709f, 1.0000724f, 1.0001201f,
+    1.0001480f, 1.0001204f, 1.0001185f, 1.0000844f, 1.0000224f, 0.9999752f,
+    0.9999997f, 1.0000969f, 1.0001076f, 1.0000756f, 1.0000700f, 1.0000854f,
+    1.0001067f, 1.0000390f, 1.0000443f, 1.0000863f, 1.0000585f, 1.0000352f,
+    1.0000677f, 1.0001081f, 1.0000537f, 1.0000199f, 1.0000308f, 1.0000585f,
+    1.0000479f, 1.0000304f, 1.0000751f, 1.0000710f, 1.0000152f, 1.0000083f,
+    1.0000342f, 1.0000530f, 1.0000543f, 1.0000442f, 1.0000680f, 1.0000753f,
+    1.0000369f, 1.0000117f};
+
 #ifdef HAVE_FFTW
 
 /**
@@ -1154,7 +1185,9 @@ void power_spectrum(const enum power_type type1, const enum power_type type2,
       else
         sprintf(powunits, "Mpc^3 eV cm^(-3)");
 
-      fprintf(outputfile, "# Folding %d, all lengths/volumes are comoving\n",
+      fprintf(outputfile,
+              "# Folding %d, all lengths/volumes are comoving. k-bin centres "
+              "are not corrected for the weights of the modes.\n",
               i);
       fprintf(outputfile, "# Shotnoise [%s]\n", powunits);
       fprintf(outputfile, "%g\n", shot);
@@ -1208,8 +1241,17 @@ void power_spectrum(const enum power_type type1, const enum power_type type2,
     power_init_output_file(outputfile, type1, type2, us, phys_const);
 
     for (int j = 0; j < numtot; ++j) {
+
+      float k = kcomb[j];
+
+      /* Shall we correct the position of the k-space bin
+       * to account for the different weights of the modes entering the bin? */
+      if (pow_data->shift_centre_small_k_bins && j < number_of_corrected_bins) {
+        k *= correction_shift_k_values[j];
+      }
+
       fprintf(outputfile, "%15.8f %15.8e %15.8e %15.8e\n", s->e->cosmology->z,
-              kcomb[j], (pcomb[j] - shot), shot);
+              k, (pcomb[j] - shot), shot);
     }
     fclose(outputfile);
   }
@@ -1273,6 +1315,9 @@ void power_init(struct power_spectrum_data* p, struct swift_params* params,
         "WARNING: fold factor is recommended not to exceed 6 for a "
         "mass assignment order of 3 (TSC) or below.");
 
+  p->shift_centre_small_k_bins = parser_get_opt_param_int(
+      params, "PowerSpectrum:shift_centre_small_k_bins", 1);
+
   /* Make sensible choices for the k-cuts */
   const int kcutn = (p->windoworder >= 3) ? 90 : 70;
   const int kcutleft = (int)(p->Ngrid / 256.0 * kcutn);

src/power_spectrum.h

@@ -71,6 +71,9 @@ struct power_spectrum_data {
   /*! The order of the mass assignment window */
   int windoworder;
 
+  /* Shall we correct the position of the k-space bin? */
+  int shift_centre_small_k_bins;
+
   /*! Array of component types to correlate on the "left" side */
   enum power_type* types1;