Merge branch 'lra'

api/fastpm/cosmology.h

@@ -25,7 +25,7 @@ struct FastPMCosmology {
     int N_ncdm;
     int ncdm_freestreaming;  // bool: treat ncdm as free-streaming?
     int ncdm_matterlike;     // bool: treat ncdm as matter-like?
-
+    int ncdm_linearresponse; // bool: Enable the linear response module for neutrinos.
     FastPMGrowthMode growth_mode;
     FastPMFDInterp * FDinterp;
 };

api/fastpm/gravity.h

@@ -14,7 +14,8 @@ fastpm_solver_compute_force(FastPMSolver * fastpm,
     FastPMPainter * painter,
     FastPMSofteningType dealias,
     FastPMKernelType kernel,
-    FastPMFloat * delta_k);
+    FastPMFloat * delta_k,
+    double Time);
 
 void
 gravity_apply_kernel_transfer(FastPMKernelType kernel, PM * pm, FastPMFloat * delta_k, FastPMFloat * canvas, FastPMFieldDescr field);

api/fastpm/neutrinos_lra.h

@@ -0,0 +1,67 @@
+#ifndef NEUTRINOS_LRA_H
+#define NEUTRINOS_LRA_H
+
+#include <fastpm/libfastpm.h>
+#include <bigfile.h>
+
+/** Now we want to define a static object to store all previous delta_tot.
+ * This object needs a constructor, a few private data members, and a way to be read and written from disk.
+ * nk is fixed, delta_tot, scalefact and ia are updated in get_delta_nu_update*/
+typedef struct _delta_tot_table {
+    /** Number of actually non-zero k values stored in each power spectrum*/
+    int nk;
+    /** Size of arrays allocated to store power spectra*/
+    int nk_allocated;
+    /** Maximum number of redshifts to store. Redshifts are stored every delta a = 0.01 */
+    int namax;
+    /** Number of already "recorded" time steps, i.e. scalefact[0...ia-1] is recorded.
+    * Current time corresponds to index ia (but is only recorded if sufficiently far from previous time).
+    * Caution: ia here is different from Na in get_delta_nu (Na = ia+1).*/
+    int ia;
+    /** Prefactor for use in get_delta_nu. Should be 3/2 Omega_m H^2 /c. Units are h^2 / sec / Mpc  */
+    double delta_nu_prefac;
+    /** Set to unity once the init routine has run.*/
+    int delta_tot_init_done;
+    /** Pointer to nk arrays of length namax containing the total power spectrum.*/
+    double **delta_tot;
+    /** Array of length namax containing scale factors at which the power spectrum is stored*/
+    double * scalefact;
+    /** Pointer to array of length nk storing initial neutrino power spectrum*/
+    double * delta_nu_init;
+    /** Last-seen neutrino power spectrum*/
+    double * delta_nu_last;
+    /**Pointer to array storing the effective wavenumbers for the above power spectra*/
+    double * wavenum;
+    /** Matter density excluding neutrinos*/
+    double Omeganonu;
+    /** Light speed in internal units, Mpc/s.*/
+    double light;
+    /** The time at which the simulation starts*/
+    double TimeTransfer;
+    /* Cosmology factors*/
+    FastPMCosmology * cosmo;
+} _delta_tot_table;
+
+/* Structure for the computed neutrino data.*/
+typedef struct nu_lra_power
+{
+    double * logknu;
+    double * delta_nu_ratio;
+    int size;
+    double nu_prefac;
+    gsl_interp *nu_spline;
+} nu_lra_power;
+
+/*Computes delta_nu from a CDM power spectrum.*/
+void delta_nu_from_power(nu_lra_power * nupow, FastPMFuncK* ps, FastPMCosmology * CP, const double Time);
+
+/* Load the neutrino transfer function*/
+void load_transfer_data(const double TimeTransfer, FastPMFuncK t_init_in[]);
+
+/*These functions save and load neutrino related data from the snapshots*/
+void ncdm_lr_save_neutrinos(BigFile * bf, int ThisTask);
+int ncdm_lr_read_neutrinos(BigFile * bf, int ThisTask);
+
+/*Save the neutrino power spectrum to a file*/
+void powerspectrum_nu_save(FastPMPowerSpectrum * ps, char powerspectrum_file[], double MtotbyMcdm);
+#endif

libfastpm/Makefile

@@ -29,6 +29,7 @@ LIBSOURCES = libfastpm.c \
     Ftable.c \
     FDinterp.c \
     cosmology.c \
+    neutrinos_lra.c \
     horizon.c \
     powerspectrum.c \
     logging.c \

libfastpm/gravity.c

@@ -5,6 +5,7 @@
 #include <fastpm/prof.h>
 #include <fastpm/transfer.h>
 #include <fastpm/logging.h>
+#include <fastpm/neutrinos_lra.h>
 
 #include "pmpfft.h"
 #include "pmghosts.h"
@@ -30,7 +31,7 @@ apply_grad_transfer(PM * pm, FastPMFloat * from, FastPMFloat * to, int dir, int
             fastpm_info("fourier space kernel[%d] = %g\n", i, k_finite);
         }
     }
-#pragma omp parallel 
+#pragma omp parallel
     {
         PMKIter kiter;
         ptrdiff_t * Nmesh = pm_nmesh(pm);
@@ -68,7 +69,7 @@ apply_gaussian_softening(PM * pm, FastPMFloat * from, FastPMFloat * to, double N
     /* N is rms in mesh size */
     double r0 = N * pm->BoxSize[0] / pm->Nmesh[0];
 
-#pragma omp parallel 
+#pragma omp parallel
     {
         PMKIter kiter;
         int d;
@@ -296,7 +297,7 @@ _fastpm_solver_compute_delta_k(FastPMSolver * fastpm, PM * pm, FastPMPainter * p
     pm_clear(pm, canvas);
     /* Watch out: paint paints the mass per cell;
      * divide by mean mass per cell to convert to matter overdensity, which
-     * goes into Poisson's equation. 
+     * goes into Poisson's equation.
      *
      * In this perspective, we are still operating with the dimension-less
      * Poisson's equation where the critical density factors canceled
@@ -413,13 +414,40 @@ _fastpm_solver_compute_force(FastPMSolver * fastpm,
 
 }
 
+static double
+lra_neutrinos(double k, nu_lra_power * nulra)
+{
+    /*Add neutrino power if desired*/
+    if(k <= 0)
+        return 1;
+    /* Change the units of k to match those of logkk*/
+    double logk = log(k);
+    /* Floating point roundoff and the binning means there may be a mode just beyond the box size.*/
+    if(logk < nulra->logknu[0] && logk > nulra->logknu[0]-log(2) )
+        logk = nulra->logknu[0];
+    else if( logk > nulra->logknu[nulra->size-1])
+        logk = nulra->logknu[nulra->size-1];
+    /* Note get_neutrino_powerspec returns Omega_nu / (Omega0 -OmegaNu) * delta_nu / P_cdm^1/2, which is dimensionless.
+        * So below is: M_cdm * delta_cdm (1 + Omega_nu/(Omega0-OmegaNu) (delta_nu / delta_cdm))
+        *            = M_cdm * (delta_cdm (Omega0 - OmegaNu)/Omega0 + Omega_nu/Omega0 delta_nu) * Omega0 / (Omega0-OmegaNu)
+        *            = M_cdm * Omega0 / (Omega0-OmegaNu) * (delta_cdm (1 - f_nu)  + f_nu delta_nu) )
+        *            = M_cdm * Omega0 / (Omega0-OmegaNu) * delta_t
+        *            = (M_cdm + M_nu) * delta_t
+        * This is correct for the forces, and gives the right power spectrum,
+        * once we multiply PowerSpectrum.Norm by (Omega0 / (Omega0 - OmegaNu))**2 */
+    double delta_nu = gsl_interp_eval(nulra->nu_spline, nulra->logknu, nulra->delta_nu_ratio, logk, NULL);
+    const double nufac = 1 + nulra->nu_prefac * delta_nu;
+    return nufac;
+}
+
 void
 fastpm_solver_compute_force(FastPMSolver * fastpm,
     PM * pm,
     FastPMPainter * painter,
     FastPMSofteningType dealias,
     FastPMKernelType kernel,
-    FastPMFloat * delta_k)
+    FastPMFloat * delta_k,
+    double Time)
 {
     PMGhostData * pgd[FASTPM_SOLVER_NSPECIES];
 
@@ -449,6 +477,36 @@ fastpm_solver_compute_force(FastPMSolver * fastpm,
         nacc = 3;
     }
 
+    FastPMCosmology * cosmo = fastpm->cosmology;
+    /* Transfer delta_k to delta_m with adding delta_nu*/
+    /*Computes delta_nu from a CDM power spectrum.*/
+    if(cosmo->ncdm_linearresponse) {
+        FastPMPowerSpectrum dmps[1] ={0};
+        /* calculate the neutrino power spectrum. pm_alloc for delta_k uses pm. */
+        fastpm_powerspectrum_init_from_delta(dmps, pm, delta_k, delta_k);
+        /* Need to allocate memory for this!*/
+        nu_lra_power nulra[1] = {0};
+        nulra->size = dmps->base.size;
+        nulra->logknu = malloc(sizeof(nulra->logknu[0]) * nulra->size);
+        nulra->delta_nu_ratio = malloc(sizeof(nulra->delta_nu_ratio[0]) * nulra->size);
+
+        /* Get delta_cdm_curr , which is P(k)^1/2.
+         * The actual power spectrum data is stored in base.[k,f]*/
+        int i;
+        for(i=0; i<dmps->base.size; i++)
+            dmps->base.f[i] = sqrt(dmps->base.f[i]);
+        delta_nu_from_power(nulra, &dmps->base, cosmo, Time);
+        fastpm_powerspectrum_destroy(dmps);
+        /*Initialize the interpolation for the neutrinos*/
+        nulra->nu_spline = gsl_interp_alloc(gsl_interp_linear, nulra->size);
+        // nulra->nu_acc = gsl_interp_accel_alloc();
+        gsl_interp_init(nulra->nu_spline, nulra->logknu, nulra->delta_nu_ratio, nulra->size);
+        /* Now apply the neutrino transfer function to the field stored in delta_k.*/
+        fastpm_apply_any_transfer(pm, delta_k, canvas, (fastpm_fkfunc) lra_neutrinos, &nulra);
+        free(nulra->delta_nu_ratio);
+        free(nulra->logknu);
+    }
+
     _fastpm_solver_compute_force(fastpm, pm, painter, kernel, pgd, canvas, delta_k, ACC, nacc);
 
     _fastpm_solver_destroy_ghosts(fastpm, pgd);