CoLoRe-2LPT: 2LPT velocities added for particles, tracers and beams

Makefile

@@ -9,7 +9,7 @@ OPTIONS = -Wall -Wno-format-overflow -O3 -std=c99
 #Use double precision integer (enable in general)
 DEFINEFLAGS += -D_LONGIDS
 #Use normalized bias model
-DEFINEFLAGS += -D_BIAS_MODEL_2
+DEFINEFLAGS += -D_BIAS_MODEL_3
 #Use linear bias model
 #DEFINEFLAGS += -D_BIAS_MODEL_3
 #Use new lensing method
@@ -33,8 +33,8 @@ USE_MPI = yes
 #GSL
 #GSL_INC = -I/add/path
 #GSL_LIB = -L/add/path
-GSL_INC = -I/home/alonso/include
-GSL_LIB = -L/home/alonso/lib
+GSL_INC = 
+GSL_LIB = 
 #FFTW
 FFTW_INC =
 FFTW_LIB =
@@ -45,8 +45,8 @@ FITS_LIB =
 HDF5_INC =
 HDF5_LIB =
 #libconfig
-CONF_INC =
-CONF_LIB =
+CONF_INC = -I/home/mruizh/software_libconf/include
+CONF_LIB = -L/home/mruizh/software_libconf/lib
 #healpix
 HPIX_INC =
 HPIX_LIB =

param_colore_2lpt.cfg

@@ -0,0 +1,91 @@
+global:
+{
+  #Output prefix. Output will be in prefix_<node ID>.<fits/txt>
+  prefix_out= "/lustre/home/mruizh/2lpt_paper/codes/CoLoRe-master/new_out/out";
+  #Output format. Select HDF5, FITS or ASCII
+  output_format= "FITS";
+  #Output Gaussian overdensity field at z=0?
+  output_density= false
+  #Path to power spectrum at z=0. Power spectrum file must
+  #be in CAMB format: k (h/Mpc), P(k) (Mpc/h)^3.
+  pk_filename= "/lustre/home/mruizh/CoLoRe_2LPT/CoLoRe_2LPT/input_files/Pk_Abacus_Cosmology_z0.txt"
+  #This redshift range also defines the size of the box
+  z_min= 1.6
+  z_max= 3.79
+  #RNG seed note that output will depend on number of nodes, etc not only
+  #on the RNG seed
+  seed= 1
+  #Set to true if you want to generate the theory prediction for the 3D power spectrum
+  #of the different tracers.
+  write_pred=false
+  #Intervals of redshift at which the prediction will be produced.
+  pred_dz=0.1
+  #If write_pred is true, set this to true if you just want to generate the prediction
+  #and then exit. This is also useful if you want to inspect the memory requirements
+  #before a run.
+  just_write_pred= false
+}
+
+field_par:
+{
+  #Extra Gaussian smoothing scale [Mpc/h] (set to a
+  #negative value if you don't want any smoothing)
+  r_smooth= -1
+  #Do you want to smooth the Newtonian potential as well?
+  smooth_potential= false
+  #Will use a Cartesian grid with n_grid^3 cells
+  n_grid= 4096
+  #Density field type
+  # 0-lognormal
+  # 1-1LPT
+  # 2-2LPT
+  dens_type= 2
+  #If dens_type==1 or 2, buffer size (fraction per particle)
+  lpt_buffer_fraction= 0.6
+  #If dens_type==1 or 2, scheme to interpolate particle
+  #positions into a grid
+  # 0-NGP
+  # 1-CIC
+  # 2-TSC
+  lpt_interp_type= 1
+  lpt_vzty= 1
+  #Set to 1 if you want to output the LPT particle positions
+  output_lpt= 0
+}
+
+cosmo_par:
+{
+  #Non-relativistic matter
+  omega_M= 0.3153
+  #Dark energy
+  omega_L= 0.6847
+  #Baryons
+  omega_B= 0.049301692328524445
+  #Hubble parameter (in units of 100 km/s/Mpc)
+  h= 0.6736
+  #Dark energy equation of state
+  w= -1.0
+  #Primordial scalar spectral index, used only to extrapolate
+  #P(k) at low k end (-3 used at high k end)
+  ns= 0.9649
+  #Power spectrum normalization. The input power spectrum will be
+  #renormalized to this sigma8
+  sigma_8= 0.8111
+}
+
+#For each galaxy population, create a section called srcsX, starting with X=1
+srcs1:
+{
+  #Path to N(z) file. Should contain two columns
+  # 1-> z, 2-> dN(z)/dz*dOmega
+  # with dN/dzdOmega in units of deg^-2
+  nz_filename= "/lustre/home/mruizh/CoLoRe_2LPT/CoLoRe_2LPT/input_files/Nz_qso_130618_2_colore1_hZs.txt"
+  #Path to bias file. Should contain two columns
+  # 1-> z, 2-> b(z)
+  bias_filename= "/lustre/home/mruizh/CoLoRe_2LPT/CoLoRe_2LPT/input_files/Bz_clustering.txt"
+  threshold_filename= "/lustre/home/mruizh/CoLoRe_2LPT/CoLoRe_2LPT/input_files/Tz_clustering.txt"
+  #Do you want to include shear ellipticities?
+  include_lensing= false
+  #Do you want to store line-of-sight skewers for each object?
+  store_skewers= true
+}

param_example.cfg

@@ -48,7 +48,10 @@ field_par:
   # 1-CIC
   # 2-TSC
   lpt_interp_type= 1
-  #Set to 1 if you want to output the LPT particle positions
+
+  lpt_vzty= 0 #Set to 1 if you want to use 2LPT velocities with 2LPT densities. If used, only CIC available (for now) for this case!
+
+  #Set to 1 if you want to output the LPT particle positions (if using 2LPT, also velocities will be in the output!)
   output_lpt= 0
 }
 

src/beaming.c

@@ -67,14 +67,20 @@ static void get_element(ParamCoLoRe *par,long ix,long iy,long iz,
 
   //Get velocity
   if(flag_return & RETURN_VEL) {
-    v[0]=par->grid_npot[ix_hi+iy_0+iz_0]-par->grid_npot[ix_lo+iy_0+iz_0];
-    v[1]=par->grid_npot[ix_0+iy_hi+iz_0]-par->grid_npot[ix_0+iy_lo+iz_0];
-    if(iz==0)
-      v[2]=par->grid_npot[ix_0+iy_0+iz_hi]-par->slice_left[ix_0+iy_0];
-    else if(iz==par->nz_here-1)
-      v[2]=par->slice_right[ix_0+iy_0]-par->grid_npot[ix_0+iy_0+iz_lo];
-    else
-      v[2]=par->grid_npot[ix_0+iy_0+iz_hi]-par->grid_npot[ix_0+iy_0+iz_lo];
+    if (par->lpt_vzty){
+       v[0]=par->grid_velx[ix_0+iy_0+iz_0];
+       v[1]=par->grid_vely[ix_0+iy_0+iz_0];
+       v[2]=par->grid_velz[ix_0+iy_0+iz_0];
+    } else {
+       v[0]=par->grid_npot[ix_hi+iy_0+iz_0]-par->grid_npot[ix_lo+iy_0+iz_0];
+       v[1]=par->grid_npot[ix_0+iy_hi+iz_0]-par->grid_npot[ix_0+iy_lo+iz_0];
+       if(iz==0)
+          v[2]=par->grid_npot[ix_0+iy_0+iz_hi]-par->slice_left[ix_0+iy_0];
+       else if(iz==par->nz_here-1)
+         v[2]=par->slice_right[ix_0+iy_0]-par->grid_npot[ix_0+iy_0+iz_lo];
+       else
+         v[2]=par->grid_npot[ix_0+iy_0+iz_hi]-par->grid_npot[ix_0+iy_0+iz_lo];
+    }
   }
 
   //Get ISW
@@ -319,6 +325,9 @@ void get_beam_properties(ParamCoLoRe *par)
 #ifdef _HAVE_MPI
   long size_slice_npot=(par->nz_max+2)*((long)(2*par->n_grid*(par->n_grid/2+1)));
   long size_slice_dens=par->nz_max*((long)(2*par->n_grid*(par->n_grid/2+1)));
+  long size_slice_velx=par->nz_max*((long)(2*par->n_grid*(par->n_grid/2+1)));
+  long size_slice_vely=par->nz_max*((long)(2*par->n_grid*(par->n_grid/2+1)));
+  long size_slice_velz=par->nz_max*((long)(2*par->n_grid*(par->n_grid/2+1)));
   flouble *buffer_sr=my_malloc(size_slice_npot*sizeof(flouble));
 #endif //_HAVE_MPI
 
@@ -330,6 +339,9 @@ void get_beam_properties(ParamCoLoRe *par)
 #endif //_DEBUG
     mpi_sendrecv_wrap(par->grid_npot,buffer_sr,size_slice_npot,i);
     mpi_sendrecv_wrap(par->grid_dens,buffer_sr,size_slice_dens,i);
+    mpi_sendrecv_wrap(par->grid_velx,buffer_sr,size_slice_velx,i);
+    mpi_sendrecv_wrap(par->grid_vely,buffer_sr,size_slice_vely,i);
+    mpi_sendrecv_wrap(par->grid_velz,buffer_sr,size_slice_velz,i);
 
 #endif //_HAVE_MPI
     par->nz_here=par->nz_all[node_i_am_now];

src/common.h

@@ -110,6 +110,8 @@
 #define BG_KZ_CSTM 1012
 #define BG_BZ_CSTM 1013
 #define BG_NORM_CSTM 1014
+#define BG_F1 1015
+#define BG_F2 1016
 
 //Density field parameters
 #define DZ_SIGMA 0.05
@@ -248,6 +250,8 @@ typedef struct {
   double growth_d2_arr[NA]; //Array of density growth factors used to compute D_d(r)
   double growth_v_arr[NA]; //Array of velocity growth factors used to compute D_v(r)
   double growth_pd_arr[NA]; //Array of potential derivative factors used to compute \dot{\phi}
+  double growth_f1_arr[NA]; // Array of 1 order growth rates used to compute f_d(r)
+  double growth_f2_arr[NA]; // Array of 2 order growth rates used to compute f_d(r)
   double ihub_arr[NA]; //Array of 1/H(z)
   double glob_idr; //1/dr, where dr is the radial comoving distance interval used in the arrays above
   // Power spectra
@@ -267,6 +271,7 @@ typedef struct {
   int smooth_potential; //Do we smooth the newtonian potential as well?
   int dens_type; //Method to produce the density field
   int lpt_interp_type;
+  int lpt_vzty; // Use 2LPT velocities (only makes sense if dens_type is 2LPT)
   double lpt_buffer_fraction; //Fraction of memory saved for buffer particles
   int output_lpt;
   unsigned int seed_rng; //RNG seed
@@ -284,6 +289,12 @@ typedef struct {
   dftw_complex *grid_dens_f; //Fourier-space grid for the density field
   flouble *grid_dens; //Real-space grid for the density field
   dftw_complex *grid_npot_f; //Fourier-space grid for the Newtonian potential
+  flouble *grid_velx; //Real-space grid for the x component of the 2lpt velocity field
+  flouble *grid_vely; //Real-space grid for the y component of the 2lpt velocity field
+  flouble *grid_velz; //Real-space grid for the z component of the 2lpt velocity field
+  dftw_complex *grid_velx_f; //Fourier-space grid for the x component of the 2lpt velocity field (not used right?)
+  dftw_complex *grid_vely_f; //Fourier-space grid for the y component of the 2lpt velocity field (not used right?)
+  dftw_complex *grid_velz_f; //Fourier-space grid for the z component of the 2lpt velocity field (not used right?)
   flouble *grid_npot; //Real-space grid for the Newtonian potential
   flouble *slice_left; //Dummy array to store grid cells coming from the left node
   flouble *slice_right; //Dummy array to store grid cells coming from the right node
@@ -448,6 +459,7 @@ flouble *compute_lensing_spacing(ParamCoLoRe *par);
 ParamCoLoRe *read_run_params(char *fname,int test_memory);
 void write_density_grid(ParamCoLoRe *par,char *prefix_dens);
 void write_lpt(ParamCoLoRe *par,unsigned long long npart,flouble *x,flouble *y,flouble *z);
+void write_lpt_wvel(ParamCoLoRe *par,unsigned long long npart,flouble *x,flouble *y,flouble *z, flouble *vx, flouble *vy, flouble *vz); //Same as above but including 2lpt velocities (merge)
 void write_srcs(ParamCoLoRe *par);
 void write_imap(ParamCoLoRe *par);
 void write_cstm(ParamCoLoRe *par);

src/cosmo.c

@@ -47,6 +47,11 @@ double get_bg(ParamCoLoRe *par,double r,int tag,int ipop)
     return f_of_r_linear(par,r,par->growth_d2_arr,par->growth_d2_arr[0],
 			 par->growth_d2_arr[NA-1]);
   }
+  else if(tag==BG_F1)
+    return f_of_r_linear(par,r,par->growth_f1_arr,par->growth_f1_arr[0],par->growth_f1_arr[NA-1]);
+  else if(tag==BG_F2) {
+    return f_of_r_linear(par,r,par->growth_f2_arr,par->growth_f2_arr[0],par->growth_f2_arr[NA-1]);
+  }
   else if(tag==BG_V1)
     return f_of_r_linear(par,r,par->growth_v_arr,par->growth_v_arr[0],par->growth_v_arr[NA-1]);
   else if(tag==BG_PD) {
@@ -717,11 +722,14 @@ void cosmo_set(ParamCoLoRe *par)
     double gz=csm_growth_factor(pars,a)/growth0;
     double om=csm_omega_m(pars,a);
     double fz=csm_f_growth(pars,a);
+    double f2z=2*pow(om, 0.5714285714285714); // Using 2LPTic Paper aproximation for f2(z) consistent with the one for D2(z)
     double hhz=csm_hubble(pars,a);
     par->z_arr_r2z[ii]=z;
     par->r_arr_r2z[ii]=r;
     par->growth_d_arr[ii]=gz;
-    par->growth_d2_arr[ii]=-0.42857142857*gz*gz*pow(om,-0.00699300699);
+    par->growth_d2_arr[ii]=-0.42857142857*gz*gz*pow(om,-0.00699300699); // Using Paper aproximation for D2(z) = -3/7*D1**2*Om**eps; consistent with the one for f2(z)
+    par->growth_f1_arr[ii]=fz; // Initialize the f1 and f2 arrays
+    par->growth_f2_arr[ii]=f2z;
     par->growth_v_arr[ii]=(gz*hhz*fz)/(par->fgrowth_0*par->hubble_0); //This is for the comoving velocity
     par->growth_pd_arr[ii]=gz*hhz*(fz-1);
     par->ihub_arr[ii]=1./hhz;