neutrinos_lra.h File Reference

#include <bigfile-mpi.h>
#include "powerspectrum.h"
#include "cosmology.h"

Include dependency graph for neutrinos_lra.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	_delta_tot_table

Typedefs
typedef struct _delta_tot_table	_delta_tot_table

Functions
void	init_neutrinos_lra (const int nk_in, const double TimeTransfer, const double TimeMax, const double Omega0, const _omega_nu *const omnu, const double UnitTime_in_s, const double UnitLength_in_cm)
void	delta_nu_from_power (struct _powerspectrum *PowerSpectrum, Cosmology *CP, const double Time, const double TimeIC)
void	petaio_save_neutrinos (BigFile *bf, int ThisTask)
void	petaio_read_neutrinos (BigFile *bf, int ThisTask)
void	petaio_read_icnutransfer (BigFile *bf, int ThisTask)
void	powerspectrum_nu_save (struct _powerspectrum *PowerSpectrum, const char *OutputDir, const char *filename, const double Time)

Typedef Documentation

_delta_tot_table

typedef struct _delta_tot_table _delta_tot_table

Definition at line 1 of file neutrinos_lra.h.

Function Documentation

delta_nu_from_power()

void delta_nu_from_power	(	struct _powerspectrum *	PowerSpectrum,
		Cosmology *	CP,
		const double	Time,
		const double	TimeIC
	)

Definition at line 134 of file neutrinos_lra.c.

{
    int i;
    /*This is done on the first timestep: we need nk_nonzero for it to work.*/
    if(!delta_tot_table.delta_tot_init_done) {
        if(delta_tot_table.ia == 0) {
            /* Compute delta_nu from the transfer functions if first entry.*/
            delta_tot_first_init(&delta_tot_table, PowerSpectrum->nonzero, PowerSpectrum->kk, PowerSpectrum->Power, TimeIC);
        }
 
        /*Initialise the first delta_nu*/
        get_delta_nu_combined(CP, &delta_tot_table, exp(delta_tot_table.scalefact[delta_tot_table.ia-1]), delta_tot_table.delta_nu_last);
        delta_tot_table.delta_tot_init_done = 1;
    }
    for(i = 0; i < PowerSpectrum->nonzero; i++)
        PowerSpectrum->logknu[i] = log(PowerSpectrum->kk[i]);
 
    double * Power_in = PowerSpectrum->Power;
    /* Rebin the input power if necessary*/
    if(delta_tot_table.nk != PowerSpectrum->nonzero) {
        Power_in = (double *) mymalloc("pkint", delta_tot_table.nk * sizeof(double));
        double * logPower = (double *) mymalloc("logpk", PowerSpectrum->nonzero * sizeof(double));
        for(i = 0; i < PowerSpectrum->nonzero; i++)
            logPower[i] = log(PowerSpectrum->Power[i]);
        gsl_interp * pkint = gsl_interp_alloc(gsl_interp_linear, PowerSpectrum->nonzero);
        gsl_interp_init(pkint, PowerSpectrum->logknu, logPower, PowerSpectrum->nonzero);
        gsl_interp_accel * pkacc = gsl_interp_accel_alloc();
        for(i = 0; i < delta_tot_table.nk; i++) {
            double logk = log(delta_tot_table.wavenum[i]);
            if(pkint->xmax < logk || pkint->xmin > logk)
                Power_in[i] = delta_tot_table.delta_tot[i][delta_tot_table.ia-1];
            else
                Power_in[i] = exp(gsl_interp_eval(pkint, PowerSpectrum->logknu, logPower, logk, pkacc));
        }
        myfree(logPower);
        gsl_interp_accel_free(pkacc);
        gsl_interp_free(pkint);
    }
 
    const double partnu = particle_nu_fraction(&CP->ONu.hybnu, Time, 0);
    /* If we get called twice with the same scale factor, do nothing: delta_nu
     * already stores the neutrino power from the current timestep.*/
    if(1 - partnu > 1e-3 && log(Time)-delta_tot_table.scalefact[delta_tot_table.ia-1] > FLOAT_ACC) {
        /*We need some estimate for delta_tot(current time) to obtain delta_nu(current time).
            Even though delta_tot(current time) is not directly used (the integrand vanishes at a = a(current)),
            it is indeed needed for interpolation */
        /*It was checked that using delta_tot(current time) = delta_cdm(current time) leads to no more than 2%
          error on delta_nu (and moreover for large k). Using the last timestep's delta_nu decreases error even more.
          So we only need one step. */
        /*This increments the number of stored spectra, although the last one is not yet final.*/
        update_delta_tot(&delta_tot_table, Time, Power_in, delta_tot_table.delta_nu_last, 0);
        /*Get the new delta_nu_curr*/
        get_delta_nu_combined(CP, &delta_tot_table, Time, delta_tot_table.delta_nu_last);
        /* Decide whether we save the current time or not */
        if (Time > exp(delta_tot_table.scalefact[delta_tot_table.ia-2]) + 0.009) {
            /* If so update delta_tot(a) correctly, overwriting current power spectrum */
            update_delta_tot(&delta_tot_table, Time, Power_in, delta_tot_table.delta_nu_last, 1);
        }
        /*Otherwise discard the last powerspectrum*/
        else
            delta_tot_table.ia--;
 
        message(0,"Done getting neutrino power: nk = %d, k = %g, delta_nu = %g, delta_cdm = %g,\n", delta_tot_table.nk, delta_tot_table.wavenum[1], delta_tot_table.delta_nu_last[1], Power_in[1]);
        /*kspace_prefac = M_nu (analytic) / M_particles */
        const double OmegaNu_nop = get_omega_nu_nopart(&CP->ONu, Time);
        const double omega_hybrid = get_omega_nu(&CP->ONu, 1) * partnu / pow(Time, 3);
        /* Omega0 - Omega in neutrinos + Omega in particle neutrinos = Omega in particles*/
        PowerSpectrum->nu_prefac = OmegaNu_nop/(delta_tot_table.Omeganonu/pow(Time,3) + omega_hybrid);
    }
    double * delta_nu_ratio = (double *) mymalloc2("dnu_rat", delta_tot_table.nk * sizeof(double));
    double * logwavenum = (double *) mymalloc2("logwavenum", delta_tot_table.nk * sizeof(double));
    gsl_interp * pkint = gsl_interp_alloc(gsl_interp_linear, delta_tot_table.nk);
    gsl_interp_accel * pkacc = gsl_interp_accel_alloc();
    /*We want to interpolate in log space*/
    for(i=0; i < delta_tot_table.nk; i++) {
        if(isnan(delta_tot_table.delta_nu_last[i]))
            endrun(2004,"delta_nu_curr=%g i=%d delta_cdm_curr=%g kk=%g\n",delta_tot_table.delta_nu_last[i],i,Power_in[i],delta_tot_table.wavenum[i]);
        /*Enforce positivity for sanity reasons*/
        if(delta_tot_table.delta_nu_last[i] < 0)
            delta_tot_table.delta_nu_last[i] = 0;
        delta_nu_ratio[i] = delta_tot_table.delta_nu_last[i]/ Power_in[i];
        logwavenum[i] = log(delta_tot_table.wavenum[i]);
    }
    if(delta_tot_table.nk != PowerSpectrum->nonzero)
        myfree(Power_in);
    gsl_interp_init(pkint, logwavenum, delta_nu_ratio, delta_tot_table.nk);
 
    /*We want to interpolate in log space*/
    for(i=0; i < PowerSpectrum->nonzero; i++) {
        if(PowerSpectrum->nonzero == delta_tot_table.nk)
            PowerSpectrum->delta_nu_ratio[i] = delta_nu_ratio[i];
        else {
            double logk = PowerSpectrum->logknu[i];
            if(logk > pkint->xmax)
                logk = pkint->xmax;
            PowerSpectrum->delta_nu_ratio[i] = gsl_interp_eval(pkint, logwavenum, delta_nu_ratio, logk, pkacc);
        }
    }
 
    gsl_interp_accel_free(pkacc);
    gsl_interp_free(pkint);
    myfree(logwavenum);
    myfree(delta_nu_ratio);
}

References CP, _delta_tot_table::delta_nu_last, _powerspectrum::delta_nu_ratio, _delta_tot_table::delta_tot, delta_tot_first_init(), _delta_tot_table::delta_tot_init_done, delta_tot_table, endrun(), FLOAT_ACC, get_delta_nu_combined(), get_omega_nu(), get_omega_nu_nopart(), _omega_nu::hybnu, _delta_tot_table::ia, _powerspectrum::kk, _transfer_init_table::logk, _powerspectrum::logknu, message(), myfree, mymalloc, mymalloc2, _delta_tot_table::nk, _powerspectrum::nonzero, _powerspectrum::nu_prefac, _delta_tot_table::Omeganonu, Cosmology::ONu, particle_nu_fraction(), _powerspectrum::Power, _delta_tot_table::scalefact, update_delta_tot(), and _delta_tot_table::wavenum.

Referenced by compute_neutrino_power().

Here is the call graph for this function:

Here is the caller graph for this function:

init_neutrinos_lra()

void init_neutrinos_lra	(	const int	nk_in,
		const double	TimeTransfer,
		const double	TimeMax,
		const double	Omega0,
		const _omega_nu *const	omnu,
		const double	UnitTime_in_s,
		const double	UnitLength_in_cm
	)

Allocates memory for delta_tot_table.

Parameters

nk_in	Number of bins stored in each power spectrum.
TimeTransfer	Scale factor of the transfer functions.
TimeMax	Final scale factor up to which we will need memory.
Omega0	Matter density at z=0.
omnu	Pointer to structure containing pre-computed tables for evaluating neutrino matter densities.
UnitTime_in_s	Time unit of the simulation in s.
UnitLength_in_cm	Length unit of the simulation in cm

Definition at line 450 of file neutrinos_lra.c.

{
   _delta_tot_table *d_tot = &delta_tot_table;
   int count;
   /*Memory allocations need to be done on all processors*/
   d_tot->nk_allocated=nk_in;
   d_tot->nk=nk_in;
   /*Store starting time*/
   d_tot->TimeTransfer = TimeTransfer;
   /* Allocate memory for delta_tot here, so that we can have further memory allocated and freed
    * before delta_tot_init is called. The number nk here should be larger than the actual value needed.*/
   /*Allocate pointers to each k-vector*/
   d_tot->namax=ceil(100*(TimeMax-TimeTransfer))+2;
   d_tot->ia=0;
   d_tot->delta_tot =(double **) mymalloc("kspace_delta_tot",nk_in*sizeof(double *));
   /*Allocate list of scale factors, and space for delta_tot, in one operation.*/
   d_tot->scalefact = (double *) mymalloc("kspace_scalefact",d_tot->namax*(nk_in+1)*sizeof(double));
   /*Allocate space for delta_nu and wavenumbers*/
   d_tot->delta_nu_last = (double *) mymalloc("kspace_delta_nu", sizeof(double) * 3 * nk_in);
   d_tot->delta_nu_init = d_tot->delta_nu_last + nk_in;
   d_tot->wavenum = d_tot->delta_nu_init + nk_in;
   /*Allocate actual data. Note that this means data can be accessed either as:
    * delta_tot[k][a] OR as
    * delta_tot[0][a+k*namax] */
   d_tot->delta_tot[0] = d_tot->scalefact+d_tot->namax;
   for(count=1; count< nk_in; count++)
        d_tot->delta_tot[count] = d_tot->delta_tot[0] + count*d_tot->namax;
   /*Allocate space for the initial neutrino power spectrum*/
   /*Setup pointer to the matter density*/
   d_tot->omnu = omnu;
   /*Set the prefactor for delta_nu, and the units system*/
   d_tot->light = LIGHTCGS * UnitTime_in_s/UnitLength_in_cm;
   d_tot->delta_nu_prefac = 1.5 *Omega0 * HUBBLE * HUBBLE * pow(UnitTime_in_s,2)/d_tot->light;
   /*Matter fraction excluding neutrinos*/
   d_tot->Omeganonu = Omega0 - get_omega_nu(omnu, 1);
}

References _delta_tot_table::delta_nu_init, _delta_tot_table::delta_nu_last, _delta_tot_table::delta_nu_prefac, _delta_tot_table::delta_tot, delta_tot_table, get_omega_nu(), HUBBLE, _delta_tot_table::ia, _delta_tot_table::light, LIGHTCGS, mymalloc, _delta_tot_table::namax, _delta_tot_table::nk, _delta_tot_table::nk_allocated, _delta_tot_table::Omeganonu, _delta_tot_table::omnu, _delta_tot_table::scalefact, _delta_tot_table::TimeTransfer, UnitLength_in_cm, and _delta_tot_table::wavenum.

Referenced by begrun(), and test_allocate_delta_tot_table().

Here is the call graph for this function:

Here is the caller graph for this function:

petaio_read_icnutransfer()

void petaio_read_icnutransfer	(	BigFile *	bf,
		int	ThisTask
	)

Definition at line 317 of file neutrinos_lra.c.

{
#pragma omp master
    {
    t_init->NPowerTable = 2;
    BigBlock bn;
    /* Read the size of the ICTransfer block.
     * If we can't read it, just set it to zero*/
    if(0 == big_file_mpi_open_block(bf, &bn, "ICTransfers", MPI_COMM_WORLD)) {
        if(0 != big_block_get_attr(&bn, "Nentry", &t_init->NPowerTable, "u8", 1))
            endrun(0, "Failed to read attr: %s\n", big_file_get_error_message());
        if(0 != big_block_mpi_close(&bn, MPI_COMM_WORLD))
            endrun(0, "Failed to close block %s\n",big_file_get_error_message());
    }
    message(0,"Found transfer function, using %d rows.\n", t_init->NPowerTable);
    t_init->logk = (double *) mymalloc2("Transfer_functions", sizeof(double) * 2*t_init->NPowerTable);
    t_init->T_nu=t_init->logk+t_init->NPowerTable;
 
    /*Defaults: a very small value*/
    t_init->logk[0] = -100;
    t_init->logk[t_init->NPowerTable-1] = 100;
 
    t_init->T_nu[0] = 1e-30;
    t_init->T_nu[t_init->NPowerTable-1] = 1e-30;
 
    /*Now read the arrays*/
    BigArray Tnu = {0};
    BigArray logk = {0};
 
    size_t dims[2] = {0, 1};
    ptrdiff_t strides[2] = {sizeof(double), sizeof(double)};
    /*The neutrino state is shared between all processors,
     *so only read on master task and broadcast*/
    if(ThisTask == 0) {
        dims[0] = t_init->NPowerTable;
    }
    big_array_init(&Tnu, t_init->T_nu, "=f8", 2, dims, strides);
    big_array_init(&logk, t_init->logk, "=f8", 2, dims, strides);
    /*This is delta_nu / delta_tot: note technically the most massive eigenstate only.
     * But if there are significant differences the mass is so low this is basically zero.*/
    petaio_read_block(bf, "ICTransfers/DELTA_NU", &Tnu, 0);
    petaio_read_block(bf, "ICTransfers/logk", &logk, 0);
    /*Also want d_{cdm+bar} / d_tot so we can get d_nu/(d_cdm+d_b)*/
    double * T_cb = (double *) mymalloc("tmp1", t_init->NPowerTable* sizeof(double));
    T_cb[0] = 1;
    T_cb[t_init->NPowerTable-1] = 1;
    BigArray Tcb = {0};
    big_array_init(&Tcb, T_cb, "=f8", 2, dims, strides);
    petaio_read_block(bf, "ICTransfers/DELTA_CB", &Tcb, 0);
    int i;
    for(i = 0; i < t_init->NPowerTable; i++)
        t_init->T_nu[i] /= T_cb[i];
    myfree(T_cb);
    /*Broadcast the arrays.*/
    MPI_Bcast(t_init->logk,2*t_init->NPowerTable,MPI_DOUBLE,0,MPI_COMM_WORLD);
    }
}

References endrun(), _transfer_init_table::logk, message(), myfree, mymalloc, mymalloc2, _transfer_init_table::NPowerTable, petaio_read_block(), t_init, _transfer_init_table::T_nu, and ThisTask.

Referenced by petaio_read_snapshot().

Here is the call graph for this function:

Here is the caller graph for this function:

petaio_read_neutrinos()

void petaio_read_neutrinos	(	BigFile *	bf,
		int	ThisTask
	)

Definition at line 376 of file neutrinos_lra.c.

{
#pragma omp master
    {
    size_t nk, ia, ik, i;
    BigBlock bn;
    if(0 != big_file_mpi_open_block(bf, &bn, "Neutrino", MPI_COMM_WORLD)) {
        endrun(0, "Failed to open block at %s:%s\n", "Neutrino",
                    big_file_get_error_message());
    }
    if(
    (0 != big_block_get_attr(&bn, "Nscale", &ia, "u8", 1)) ||
    (0 != big_block_get_attr(&bn, "Nkval", &nk, "u8", 1))) {
        endrun(0, "Failed to read attr: %s\n",
                    big_file_get_error_message());
    }
    double *delta_tot = (double *) mymalloc("tmp_nusave",ia*nk*sizeof(double));
    /*Allocate list of scale factors, and space for delta_tot, in one operation.*/
    if(0 != big_block_get_attr(&bn, "scalefact", delta_tot_table.scalefact, "f8", ia))
        endrun(0, "Failed to read attr: %s\n", big_file_get_error_message());
    if(0 != big_block_mpi_close(&bn, MPI_COMM_WORLD)) {
        endrun(0, "Failed to close block %s\n",
                    big_file_get_error_message());
    }
    BigArray deltas = {0};
    size_t dims[2] = {0, ia};
    ptrdiff_t strides[2] = {(ptrdiff_t) (sizeof(double)*ia), (ptrdiff_t)sizeof(double)};
    /*The neutrino state is shared between all processors,
     *so only read on master task and broadcast*/
    if(ThisTask == 0) {
        dims[0] = nk;
    }
    big_array_init(&deltas, delta_tot, "=f8", 2, dims, strides);
    petaio_read_block(bf, "Neutrino/Deltas", &deltas, 1);
    if(nk > 1Lu*delta_tot_table.nk_allocated || ia > 1Lu*delta_tot_table.namax)
        endrun(5, "Allocated nk %d na %d for neutrino power but need nk %d na %d\n", delta_tot_table.nk_allocated, delta_tot_table.namax, nk, ia);
    /*Save a flat memory block*/
    for(ik=0;ik<nk;ik++)
        for(i=0;i<ia;i++)
            delta_tot_table.delta_tot[ik][i] = delta_tot[ik*ia+i];
    delta_tot_table.nk = nk;
    delta_tot_table.ia = ia;
    myfree(delta_tot);
    /* Read the initial delta_nu. This is basically zero anyway,
     * so for backwards compatibility do not require it*/
    BigArray delta_nu = {0};
    dims[1] = 1;
    strides[0] = sizeof(double);
    memset(delta_tot_table.delta_nu_init, 0, delta_tot_table.nk);
    big_array_init(&delta_nu, delta_tot_table.delta_nu_init, "=f8", 2, dims, strides);
    petaio_read_block(bf, "Neutrino/DeltaNuInit", &delta_nu, 0);
    /* Read the k values*/
    BigArray kvalue = {0};
    memset(delta_tot_table.wavenum, 0, delta_tot_table.nk);
    big_array_init(&kvalue, delta_tot_table.wavenum, "=f8", 2, dims, strides);
    petaio_read_block(bf, "Neutrino/kvalue", &kvalue, 0);
 
    /*Broadcast the arrays.*/
    MPI_Bcast(&(delta_tot_table.ia), 1,MPI_INT,0,MPI_COMM_WORLD);
    MPI_Bcast(&(delta_tot_table.nk), 1,MPI_INT,0,MPI_COMM_WORLD);
    MPI_Bcast(delta_tot_table.delta_nu_init,delta_tot_table.nk,MPI_DOUBLE,0,MPI_COMM_WORLD);
    MPI_Bcast(delta_tot_table.wavenum,delta_tot_table.nk,MPI_DOUBLE,0,MPI_COMM_WORLD);
 
    if(delta_tot_table.ia > 0) {
        /*Broadcast data for scalefact and delta_tot, Delta_tot is allocated as the same block of memory as scalefact.
          Not all this memory will actually have been used, but it is easiest to bcast all of it.*/
        MPI_Bcast(delta_tot_table.scalefact,delta_tot_table.namax*(delta_tot_table.nk+1),MPI_DOUBLE,0,MPI_COMM_WORLD);
    }
    }
}

References _delta_tot_table::delta_nu_init, _delta_tot_table::delta_tot, delta_tot_table, endrun(), _delta_tot_table::ia, myfree, mymalloc, _delta_tot_table::namax, _delta_tot_table::nk, _delta_tot_table::nk_allocated, petaio_read_block(), _delta_tot_table::scalefact, ThisTask, and _delta_tot_table::wavenum.

Referenced by petaio_read_snapshot().

Here is the call graph for this function:

Here is the caller graph for this function:

petaio_save_neutrinos()

void petaio_save_neutrinos	(	BigFile *	bf,
		int	ThisTask
	)

Definition at line 265 of file neutrinos_lra.c.

{
#pragma omp master
    {
    double * scalefact = delta_tot_table.scalefact;
    size_t nk = delta_tot_table.nk, ia = delta_tot_table.ia;
    size_t ik, i;
    double * delta_tot = (double *) mymalloc("tmp_delta",nk * ia * sizeof(double));
    /*Save a flat memory block*/
    for(ik=0;ik< nk;ik++)
        for(i=0;i< ia;i++)
            delta_tot[ik*ia+i] = delta_tot_table.delta_tot[ik][i];
 
    BigBlock bn;
    if(0 != big_file_mpi_create_block(bf, &bn, "Neutrino", NULL, 0, 0, 0, MPI_COMM_WORLD)) {
        endrun(0, "Failed to create block at %s:%s\n", "Neutrino",
                big_file_get_error_message());
    }
    if ( (0 != big_block_set_attr(&bn, "Nscale", &ia, "u8", 1)) ||
       (0 != big_block_set_attr(&bn, "scalefact", scalefact, "f8", ia)) ||
        (0 != big_block_set_attr(&bn, "Nkval", &nk, "u8", 1)) ) {
        endrun(0, "Failed to write neutrino attributes %s\n",
                    big_file_get_error_message());
    }
    if(0 != big_block_mpi_close(&bn, MPI_COMM_WORLD)) {
        endrun(0, "Failed to close block %s\n",
                    big_file_get_error_message());
    }
    BigArray deltas = {0};
    size_t dims[2] = {0 , ia};
    /*The neutrino state is shared between all processors,
     *so only write on master task*/
    if(ThisTask == 0) {
        dims[0] = nk;
    }
    ptrdiff_t strides[2] = {(ptrdiff_t) (sizeof(double) * ia), (ptrdiff_t) sizeof(double)};
    big_array_init(&deltas, delta_tot, "=f8", 2, dims, strides);
    petaio_save_block(bf, "Neutrino/Deltas", &deltas, 1);
    myfree(delta_tot);
    /*Now write the initial neutrino power*/
    BigArray delta_nu = {0};
    dims[1] = 1;
    strides[0] = sizeof(double);
    big_array_init(&delta_nu, delta_tot_table.delta_nu_init, "=f8", 2, dims, strides);
    petaio_save_block(bf, "Neutrino/DeltaNuInit", &delta_nu, 1);
    /*Now write the k values*/
    BigArray kvalue = {0};
    big_array_init(&kvalue, delta_tot_table.wavenum, "=f8", 2, dims, strides);
    petaio_save_block(bf, "Neutrino/kvalue", &kvalue, 1);
    }
}

References _delta_tot_table::delta_nu_init, _delta_tot_table::delta_tot, delta_tot_table, endrun(), _delta_tot_table::ia, myfree, mymalloc, _delta_tot_table::nk, petaio_save_block(), _delta_tot_table::scalefact, ThisTask, and _delta_tot_table::wavenum.

Referenced by petaio_save_snapshot().

Here is the call graph for this function:

Here is the caller graph for this function:

powerspectrum_nu_save()

void powerspectrum_nu_save	(	struct _powerspectrum *	PowerSpectrum,
		const char *	OutputDir,
		const char *	filename,
		const double	Time
	)

Definition at line 240 of file neutrinos_lra.c.

{
    int i;
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask != 0)
        return;
 
    char * fname = fastpm_strdup_printf("%s/%s-%0.4f.txt", OutputDir, filename, Time);
    /* Now save the neutrino power spectrum*/
    FILE * fp = fopen(fname, "w");
    fprintf(fp, "# in Mpc/h Units \n");
    fprintf(fp, "# k P_nu(k) Nmodes\n");
    fprintf(fp, "# a= %g\n", Time);
    fprintf(fp, "# nk = %d\n", PowerSpectrum->nonzero);
    for(i = 0; i < PowerSpectrum->nonzero; i++){
        fprintf(fp, "%g %g %ld\n", PowerSpectrum->kk[i], pow(delta_tot_table.delta_nu_last[i],2), PowerSpectrum->Nmodes[i]);
    }
    fclose(fp);
    myfree(fname);
    /*Clean up the neutrino memory now we saved the power spectrum.*/
    gsl_interp_free(PowerSpectrum->nu_spline);
    gsl_interp_accel_free(PowerSpectrum->nu_acc);
}

References _delta_tot_table::delta_nu_last, delta_tot_table, fastpm_strdup_printf(), _powerspectrum::kk, myfree, _powerspectrum::Nmodes, _powerspectrum::nonzero, _powerspectrum::nu_acc, _powerspectrum::nu_spline, and ThisTask.

Referenced by gravpm_force().

Here is the call graph for this function:

Here is the caller graph for this function: