Main.cpp File Reference

Go to the source code of this file.

Typedefs
template<int N>
using	FFTWGrid = FML::GRID::FFTWGrid< N >
using	Particle = COLAParticle< NDIM >
using	DVector = FML::SOLVERS::ODESOLVER::DVector
using	ODEFunction = FML::SOLVERS::ODESOLVER::ODEFunction
using	ODESolver = FML::SOLVERS::ODESOLVER::ODESolver
using	FoFHalo = FML::FOF::FoFHalo< Particle, NDIM >

Functions
double	HoverH0 (double a)
int	main ()

Variables
const int	NDIM = 3
const double	box = 1000.0
const double	zini = 20.0
const double	aini = 1.0 / (1.0 + zini)
const double	aend = 1.0
const int	nsteps = 20
const bool	fix_amplitude = true
const std::string	type_of_random_field = "gaussian"
const double	fNL = 100.0
const unsigned int	random_seed = 1234
const int	Nmesh = 256
const int	Nmesh_IC = 128
const int	Npart_1D = 128
const double	buffer_factor = 1.5
const bool	COLA = true
const int	LPT_order_IC = 3
const int	n_min_FoF_group = 20
const int	Ngrid_max_FoF = 256
const double	FoFBoundarySizeMpch
const double	linking_length = 0.2
const int	nbin_bispectrum = 16
const double	apofk = 1.0
const std::string	filename = "pofk.txt"
const std::string	interpolation_method = "CIC"
const double	OmegaM = 1.0
const double	h = 0.7
const double	A_s = 2e-9
const double	n_s = 0.96
const double	kpivot = 0.05 / h
FML::INTERPOLATION::SPLINE::Spline	D_1LPT_of_loga
FML::INTERPOLATION::SPLINE::Spline	D_2LPT_of_loga
FML::INTERPOLATION::SPLINE::Spline	D_3LPTa_of_loga
FML::INTERPOLATION::SPLINE::Spline	D_3LPTb_of_loga

Typedef Documentation

DVector

using DVector = FML::SOLVERS::ODESOLVER::DVector

Definition at line 82 of file Main.cpp.

FFTWGrid

template<int N>

using FFTWGrid = FML::GRID::FFTWGrid<N>

Definition at line 80 of file Main.cpp.

FoFHalo

using FoFHalo = FML::FOF::FoFHalo<Particle, NDIM>

Definition at line 85 of file Main.cpp.

ODEFunction

using ODEFunction = FML::SOLVERS::ODESOLVER::ODEFunction

Definition at line 83 of file Main.cpp.

ODESolver

using ODESolver = FML::SOLVERS::ODESOLVER::ODESolver

Definition at line 84 of file Main.cpp.

Particle

using Particle = COLAParticle<NDIM>

Definition at line 81 of file Main.cpp.

Function Documentation

HoverH0()

double HoverH0

(

double

a

)

Definition at line 76 of file Main.cpp.

{ return std::sqrt(OmegaM / (a * a * a) + 1.0 - OmegaM); };

main()

int main

(

)

Definition at line 87 of file Main.cpp.

           {
 
    //========================================================
    // If COLA then the particles must have get_D_nLPT method(s)
    // This determines the LPT order we use with COLA
    // The IC method will store the relevant stuff in the particles
    //========================================================
    int LPT_order_COLA = 0;
    if (COLA) {
        if constexpr (FML::PARTICLE::has_get_D_1LPT<Particle>())
            LPT_order_COLA = 1;
        if constexpr (FML::PARTICLE::has_get_D_2LPT<Particle>())
            LPT_order_COLA = 2;
        if constexpr (FML::PARTICLE::has_get_D_3LPTa<Particle>() or FML::PARTICLE::has_get_D_3LPTb<Particle>())
            LPT_order_COLA = 3;
        LPT_order_COLA = std::min(LPT_order_COLA, LPT_order_IC);
    }
 
    //========================================================
    // We need a spline of the LPT growth factor(s)
    // Normalized such that D1LPT(zini) = 1
    //========================================================
    FML::COSMOLOGY::LPT::compute_LPT_growth_factors_LCDM(
        OmegaM, zini, HoverH0, D_1LPT_of_loga, D_2LPT_of_loga, D_3LPTa_of_loga, D_3LPTb_of_loga);
 
    //============================================================
    // Read k, P(k) (File has units of h/Mpc, (Mpc/h)^3 but what matters is that it has the same
    // units as the boxsize)
    //============================================================
    std::ifstream fp(filename.c_str());
    if (!fp.is_open())
        throw std::runtime_error("Cannot open file: " + filename);
    std::vector<double> logk;
    std::vector<double> logpofk;
    for (;;) {
        double kin, pofkin;
        fp >> kin;
        if (fp.eof())
            break;
        fp >> pofkin;
        logk.push_back(std::log(kin));
        logpofk.push_back(std::log(pofkin));
    }
    // Make a spline and the P(k) function needed below
    FML::INTERPOLATION::SPLINE::Spline logpofk_spline(logk, logpofk);
 
    //============================================================
    // Power spectrum scaled back to the initial time using growth-factors
    // For non-gaussian initial conditions we also need the primordial
    // power-spectrum
    //============================================================
    const double DiniOverD = D_1LPT_of_loga(std::log(aini)) / D_1LPT_of_loga(std::log(apofk));
    auto Pofk_of_kBox_over_volume = [&](double kBox) {
        const double logk = std::log(kBox / box);
        const double volume = std::pow(box, NDIM);
        return DiniOverD * DiniOverD * std::exp(logpofk_spline(logk)) / volume;
    };
    auto Pofk_of_kBox_over_volume_primordial = [&](double kBox) {
        const double k = kBox / box;
        return 2.0 * M_PI * M_PI / (kBox * kBox * kBox) * A_s * std::pow(k / kpivot, n_s - 1.0);
    };
    auto Pofk_of_kBox_over_Pofk_primordial = [&](double kBox) {
        return Pofk_of_kBox_over_volume(kBox) / Pofk_of_kBox_over_volume_primordial(kBox);
    };
 
    //============================================================
    // Hubble function and growth-rate(s) needed to compute IC
    //============================================================
    auto H_over_H0_of_loga = [&](double loga) { return HoverH0(std::exp(loga)); };
    std::function<double(double)> growth_rate_f_1LPT_of_loga = [&](double loga) {
        return D_1LPT_of_loga.deriv_x(loga) / D_1LPT_of_loga(loga);
    };
    std::function<double(double)> growth_rate_f_2LPT_of_loga = [&](double loga) {
        return D_2LPT_of_loga.deriv_x(loga) / D_2LPT_of_loga(loga);
    };
    std::function<double(double)> growth_rate_f_3LPTa_of_loga = [&](double loga) {
        return D_3LPTa_of_loga.deriv_x(loga) / D_3LPTa_of_loga(loga);
    };
    std::function<double(double)> growth_rate_f_3LPTb_of_loga = [&](double loga) {
        return D_3LPTb_of_loga.deriv_x(loga) / D_3LPTb_of_loga(loga);
    };
    std::vector<std::function<double(double)>> growth_rate_f_of_loga;
    growth_rate_f_of_loga.push_back(growth_rate_f_1LPT_of_loga);
    growth_rate_f_of_loga.push_back(growth_rate_f_2LPT_of_loga);
    growth_rate_f_of_loga.push_back(growth_rate_f_3LPTa_of_loga);
    growth_rate_f_of_loga.push_back(growth_rate_f_3LPTb_of_loga);
 
    //============================================================
    // Set up random number generator
    //============================================================
    std::shared_ptr<FML::RANDOM::RandomGenerator> rng = std::make_shared<FML::RANDOM::RandomGenerator>();
    rng->set_seed(random_seed);
 
    //============================================================
    // Generate initial conditions
    //============================================================
 
    std::vector<double> velocity_norms(4);
    for (int i = 0; i < 4; i++) {
        velocity_norms[i] =
            100.0 * box * H_over_H0_of_loga(std::log(aini)) * aini * growth_rate_f_of_loga[i](std::log(aini));
    }
 
    FML::PARTICLE::MPIParticles<Particle> part;
    FML::NBODY::NBodyInitialConditions<NDIM, Particle>(part,
                                                       Npart_1D,
                                                       buffer_factor,
                                                       Nmesh_IC,
                                                       fix_amplitude,
                                                       rng.get(),
                                                       Pofk_of_kBox_over_Pofk_primordial,
                                                       Pofk_of_kBox_over_volume_primordial,
                                                       LPT_order_IC,
                                                       type_of_random_field,
                                                       fNL,
                                                       box,
                                                       zini,
                                                       velocity_norms);
 
    //============================================================
    // In the COLA frame v=0 so reset velocities
    //============================================================
    if (COLA) {
        // Loop over all active particles
        for (auto & p : part) {
            auto * vel = FML::PARTICLE::GetVel(p);
            for (int idim = 0; idim < NDIM; idim++) {
                vel[idim] = 0.0;
            }
        }
    }
 
    //============================================================
    // Instead of stepping in dt_code = dt H0/a^2 we have unit time-steps in a. Here we
    // compute the integral of dt_code/da over one time-step
    //============================================================
 
    // The prefactor to delta in the Poisson equation
    auto poisson_factor = [&](double a) { return 1.5 * OmegaM * a; };
 
    auto compute_vel_timestep = [&](double alow, double ahigh) {
        // For the kick step we have an extra 'a' coming from the force via the Poisson equation
        const double amid = (ahigh + alow) / 2.;
        ODEFunction deriv = [&](double a, [[maybe_unused]] const double * t, double * dtda) {
            dtda[0] = 1.0 / (a * a * a * HoverH0(a)) * poisson_factor(a) / poisson_factor(amid);
            return GSL_SUCCESS;
        };
 
        // Solve the integral
        DVector tini{0.0};
        DVector avec{alow, ahigh};
        ODESolver ode;
        ode.solve(deriv, avec, tini);
        return ode.get_final_data()[0];
    };
 
    auto compute_pos_timestep = [&](double alow, double ahigh) {
        ODEFunction deriv = [&](double a, [[maybe_unused]] const double * t, double * dtda) {
            dtda[0] = 1.0 / (a * a * a * HoverH0(a));
            return GSL_SUCCESS;
        };
        // Solve the integral
        DVector tini{0.0};
        DVector avec{alow, ahigh};
        ODESolver ode;
        ode.solve(deriv, avec, tini);
        return ode.get_final_data()[0];
    };
 
    //============================================================
    // Main timestepping loop
    //============================================================
    for (int i = 0; i <= nsteps; i++) {
 
        // Scale factor for velocity and positions. Equidistance steps in 'a'
        const double delta_aexp = (aend - aini) / double(nsteps);
        const double apos_old = aini + delta_aexp * i;
        const double apos_new = (i == nsteps) ? aend : apos_old + delta_aexp;
        const double avel_old = (i == 0) ? aini : apos_old - 0.5 * delta_aexp;
        const double avel_new = (i == 0) ? aini + 0.5 * delta_aexp : (i == nsteps ? aend : apos_old + 0.5 * delta_aexp);
 
        // Timesteps: first step we displace pos and vel by 1/2 timestep to do leapfrog
        // Last time-step we synchronize back by only moving vel
        const double delta_time_pos = compute_pos_timestep(apos_old, apos_new);
        const double delta_time_vel = compute_vel_timestep(avel_old, avel_new);
        const double norm_poisson_equation = poisson_factor(apos_old);
 
        // Particles -> density field(x) -> density fielk(k)
        auto nleftright = FML::INTERPOLATION::get_extra_slices_needed_for_density_assignment(interpolation_method);
        FFTWGrid<NDIM> density_grid_fourier(Nmesh, nleftright.first, nleftright.second);
        FML::INTERPOLATION::particles_to_grid(part.get_particles_ptr(),
                                              part.get_npart(),
                                              part.get_npart_total(),
                                              density_grid_fourier,
                                              interpolation_method);
        density_grid_fourier.fftw_r2c();
 
        // Here we can bin up a power-spectrum if we want (its basically free)
        // Here we can also add contribution from radiation on large scales and linear massless neutrinos etc...
 
        // Density field -> force
        std::array<FFTWGrid<NDIM>, NDIM> force_real;
        FML::NBODY::compute_force_from_density_fourier<NDIM>(
            density_grid_fourier, force_real, interpolation_method, norm_poisson_equation);
 
        // Update velocity of particles (frees force)
        FML::NBODY::KickParticles<NDIM>(force_real, part, delta_time_vel, interpolation_method);
 
        // Add on COLA displacement and velocity
        if (COLA) {
            const double D1new = D_1LPT_of_loga(std::log(apos_new));
            const double D1old = D_1LPT_of_loga(std::log(apos_old));
            const double D1ini = D_1LPT_of_loga(std::log(aini));
            const double fac1_pos = (D1new - D1old) / D1ini;
            const double fac1_vel = -norm_poisson_equation * D1old / D1ini * delta_time_vel;
 
            const double D2new = D_2LPT_of_loga(std::log(apos_new));
            const double D2old = D_2LPT_of_loga(std::log(apos_old));
            const double D2ini = D_2LPT_of_loga(std::log(aini));
            const double fac2_pos = (D2new - D2old) / D2ini;
            const double fac2_vel = -norm_poisson_equation * (D2old - D1old * D1old) / D2ini * delta_time_vel;
 
            const double D3anew = D_3LPTa_of_loga(std::log(apos_new));
            const double D3aold = D_3LPTa_of_loga(std::log(apos_old));
            const double D3aini = D_3LPTa_of_loga(std::log(aini));
            const double fac3a_pos = (D3anew - D3aold) / D3aini;
            const double fac3a_vel =
                -norm_poisson_equation * (D3aold - 2.0 * D1old * D1old * D1old) / D3aini * delta_time_vel;
 
            const double D3bnew = D_3LPTb_of_loga(std::log(apos_new));
            const double D3bold = D_3LPTb_of_loga(std::log(apos_old));
            const double D3bini = D_3LPTb_of_loga(std::log(aini));
            const double fac3b_pos = (D3bnew - D3bold) / D3bini;
            const double fac3b_vel =
                -norm_poisson_equation * (D3bold + D1old * D1old * D1old - D1old * D2old) / D3bini * delta_time_vel;
 
            // Loop over all active particles
            for (auto & p : part) {
                auto * pos = FML::PARTICLE::GetPos(p);
                auto * vel = FML::PARTICLE::GetVel(p);
 
                if (LPT_order_COLA >= 1) {
                    auto * D1 = FML::PARTICLE::GetD_1LPT(p);
                    for (int idim = 0; idim < NDIM; idim++) {
                        pos[idim] += D1[idim] * fac1_pos;
                        vel[idim] += D1[idim] * fac1_vel;
                    }
                }
 
                if (LPT_order_COLA >= 2) {
                    auto * D2 = FML::PARTICLE::GetD_2LPT(p);
                    for (int idim = 0; idim < NDIM; idim++) {
                        pos[idim] += D2[idim] * fac2_pos;
                        vel[idim] += D2[idim] * fac2_vel;
                    }
                }
 
                if constexpr (FML::PARTICLE::has_get_D_3LPTa<Particle>())
                    if (LPT_order_COLA >= 3) {
                        auto * D3a = FML::PARTICLE::GetD_3LPTa(p);
                        for (int idim = 0; idim < NDIM; idim++) {
                            pos[idim] += D3a[idim] * fac3a_pos;
                            vel[idim] += D3a[idim] * fac3a_vel;
                        }
                    }
 
                if constexpr (FML::PARTICLE::has_get_D_3LPTb<Particle>())
                    if (LPT_order_COLA >= 3) {
                        auto * D3b = FML::PARTICLE::GetD_3LPTb(p);
                        for (int idim = 0; idim < NDIM; idim++) {
                            pos[idim] += D3b[idim] * fac3b_pos;
                            vel[idim] += D3b[idim] * fac3b_vel;
                        }
                    }
            }
        }
 
        // Move particles (this does communication)
        FML::NBODY::DriftParticles<NDIM, Particle>(part, delta_time_pos);
 
        if (FML::ThisTask == 0)
            std::cout << i << " / " << nsteps << " Positions at a: " << apos_new << " Velocities at a: " << avel_new
                      << "\n";
 
        // Print info
        part.info();
    }
 
    //============================================================
    // For COLA we must add on the LPT velocity as we solve in
    // the COLA frame
    //============================================================
 
    if (COLA) {
        const double D1end = D_1LPT_of_loga(std::log(aend));
        const double D1ini = D_1LPT_of_loga(std::log(aini));
        const double D2end = D_2LPT_of_loga(std::log(aend));
        const double D2ini = D_2LPT_of_loga(std::log(aini));
        const double D3aend = D_3LPTa_of_loga(std::log(aend));
        const double D3aini = D_3LPTa_of_loga(std::log(aini));
        const double D3bend = D_3LPTb_of_loga(std::log(aend));
        const double D3bini = D_3LPTb_of_loga(std::log(aini));
        const double vfac_1LPT =
            D1end / D1ini * growth_rate_f_1LPT_of_loga(std::log(aend)) * aend * aend * HoverH0(aend);
        const double vfac_2LPT =
            D2end / D2ini * growth_rate_f_2LPT_of_loga(std::log(aend)) * aend * aend * HoverH0(aend);
        const double vfac_3LPTa =
            D3aend / D3aini * growth_rate_f_3LPTa_of_loga(std::log(aend)) * aend * aend * HoverH0(aend);
        const double vfac_3LPTb =
            D3bend / D3bini * growth_rate_f_3LPTb_of_loga(std::log(aend)) * aend * aend * HoverH0(aend);
 
        // Here we add on the COLA velocity kick
        for (auto & p : part) {
            auto * vel = FML::PARTICLE::GetVel(p);
 
            if constexpr (FML::PARTICLE::has_get_D_1LPT<Particle>())
                if (LPT_order_COLA >= 1) {
                    auto * D1 = FML::PARTICLE::GetD_1LPT(p);
                    for (int idim = 0; idim < NDIM; idim++) {
                        vel[idim] += D1[idim] * vfac_1LPT;
                    }
                }
 
            if constexpr (FML::PARTICLE::has_get_D_2LPT<Particle>())
                if (LPT_order_COLA >= 2) {
                    auto * D2 = FML::PARTICLE::GetD_2LPT(p);
                    for (int idim = 0; idim < NDIM; idim++) {
                        vel[idim] += D2[idim] * vfac_2LPT;
                    }
                }
 
            if constexpr (FML::PARTICLE::has_get_D_3LPTa<Particle>())
                if (LPT_order_COLA >= 3) {
                    auto * D3a = FML::PARTICLE::GetD_3LPTa(p);
                    for (int idim = 0; idim < NDIM; idim++) {
                        vel[idim] += D3a[idim] * vfac_3LPTa;
                    }
                }
 
            if constexpr (FML::PARTICLE::has_get_D_3LPTb<Particle>())
                if (LPT_order_COLA >= 3) {
                    auto * D3b = FML::PARTICLE::GetD_3LPTb(p);
                    for (int idim = 0; idim < NDIM; idim++) {
                        vel[idim] += D3b[idim] * vfac_3LPTb;
                    }
                }
        }
    }
 
    //============================================================
    // Compute power-spectrum at final redshift
    //============================================================
 
    if (FML::ThisTask == 0) {
        std::cout << "Compute power spectrum\n";
    }
 
    const bool interlacing = true;
    const std::string density_assignment = "CIC";
    FML::CORRELATIONFUNCTIONS::PowerSpectrumBinning<NDIM> pofk(Nmesh / 2);
    FML::CORRELATIONFUNCTIONS::compute_power_spectrum<NDIM>(Nmesh,
                                                            part.get_particles_ptr(),
                                                            part.get_npart(),
                                                            part.get_npart_total(),
                                                            pofk,
                                                            density_assignment,
                                                            interlacing);
    pofk.scale(box);
    if (FML::ThisTask == 0) {
        std::ofstream fp("pofk_out.txt");
        fp << "\n# k   (h/Mpc)      P(k)  (Mpc/h)^3\n";
        for (int i = 0; i < pofk.n; i++) {
            fp << pofk.kbin[i] << " " << pofk.pofk[i] << "\n";
        }
    }
 
    //============================================================
    // Compute power-spectrum multipoles at final redshift
    //============================================================
 
    // Factor to convert from velocity v_code = a^2 dxdt / (H0 Box) to the RSD shift (dxdt / H) in units of the box
    const double velocity_to_displacement = 1.0 / (aend * aend * HoverH0(aend));
    const int ellmax = 4;
    std::vector<FML::CORRELATIONFUNCTIONS::PowerSpectrumBinning<NDIM>> Pells(1 + ellmax, Nmesh / 2);
    FML::CORRELATIONFUNCTIONS::compute_power_spectrum_multipoles<NDIM>(
        Nmesh, part, velocity_to_displacement, Pells, density_assignment, interlacing);
 
    // To physical units
    for (size_t ell = 0; ell < Pells.size(); ell++) {
        Pells[ell].scale(box);
    }
 
    // Output: k, P0, P0_kaiser P2, P2_kaiser, P4, P4_kaiser
    if (FML::ThisTask == 0) {
        FML::assert_mpi(Pells[0].n == pofk.n, "We require P(k) and Pell(k) to have the same bins");
        std::ofstream fp("pofk_multipoles_out.txt");
        const double f = growth_rate_f_1LPT_of_loga(std::log(aend));
        for (int i = 0; i < Pells[0].n; i++) {
            fp << Pells[0].k[i] << " ";
            for (size_t ell = 0; ell < Pells.size(); ell += 2) {
                fp << Pells[ell].pofk[i] << " ";
 
                // Output the Kaiser limit relation between P(k) and P_ell(k)
                if (ell == 0)
                    fp << pofk.pofk[i] * (1.0 + 2. * f / 3. + f * f / 5.) << " ";
                if (ell == 2)
                    fp << pofk.pofk[i] * (4. * f / 3. + 4. * f * f / 7.) << " ";
                if (ell == 4)
                    fp << pofk.pofk[i] * (8. * f * f / 35.) << " ";
            }
            fp << "\n";
        }
        fp << "\n";
    }
 
    //============================================================
    // Compute bispectrum
    //============================================================
    FML::CORRELATIONFUNCTIONS::BispectrumBinning<NDIM> bofk(
        0.0,
        2.0 * M_PI * Nmesh / 2,
        nbin_bispectrum,
        FML::CORRELATIONFUNCTIONS::BispectrumBinning<NDIM>::LINEAR_SPACING);
    FML::CORRELATIONFUNCTIONS::compute_bispectrum(Nmesh,
                                                  part.get_particles_ptr(),
                                                  part.get_npart(),
                                                  part.get_npart_total(),
                                                  bofk,
                                                  density_assignment,
                                                  interlacing);
    bofk.scale(box);
 
    if (FML::ThisTask == 0) {
        std::ofstream fp("bofk_out.txt");
        fp << "# k1         k2         k3         B123         Q123         N123\n";
        for (int i = 0; i < nbin_bispectrum; i++) {
            for (int j = 0; j < nbin_bispectrum; j++) {
                for (int k = 0; k < nbin_bispectrum; k++) {
                    double k1 = bofk.kmean[i];
                    double k2 = bofk.kmean[j];
                    double k3 = bofk.kmean[k];
                    double B123 = bofk.get_spectrum(i, j, k);
                    double Q123 = bofk.get_reduced_spectrum(i, j, k);
                    double N123 = bofk.get_bincount(i, j, k);
                    fp << k1 << " " << k2 << " " << k3 << " " << B123 << " " << Q123 << " " << N123 << "\n";
                }
            }
        }
    }
 
    //============================================================
    // Locate halos
    //============================================================
 
    part.communicate_particles();
    std::vector<FoFHalo> FoFGroups;
    double bufferlength = FoFBoundarySizeMpch / box;
    FML::FOF::FriendsOfFriends<Particle, NDIM>(part.get_particles_ptr(),
                                               part.get_npart(),
                                               linking_length,
                                               n_min_FoF_group,
                                               true,
                                               bufferlength,
                                               FoFGroups,
                                               Ngrid_max_FoF);
 
    // Output halos task by task
    for (int i = 0; i < FML::NTasks; i++) {
        if (i == FML::ThisTask) {
            std::cout << "\n# Found " << FoFGroups.size() << " halos on task " << FML::ThisTask << "\n";
            std::ofstream fp("halo_out.txt");
            for (auto & g : FoFGroups) {
                if (g.np > 0) {
                    fp << g.np << " ";
                    for (int idim = 0; idim < NDIM; idim++)
                        fp << g.pos[idim] * box << " ";
                    fp << "\n";
                }
            }
        }
#ifdef USE_MPI
        MPI_Barrier(MPI_COMM_WORLD);
#endif
    }
 
#ifdef MEMORY_LOGGING
    // Show full info about allocations we as function of time
    FML::MemoryLog::get()->print();
#endif
}

Variable Documentation

A_s

const double A_s = 2e-9

Definition at line 67 of file Main.cpp.

aend

const double aend = 1.0

Definition at line 42 of file Main.cpp.

aini

const double aini = 1.0 / (1.0 + zini)

Definition at line 41 of file Main.cpp.

apofk

const double apofk = 1.0

Definition at line 62 of file Main.cpp.

box

const double box = 1000.0

Definition at line 39 of file Main.cpp.

buffer_factor

const double buffer_factor = 1.5

Definition at line 51 of file Main.cpp.

COLA

const bool COLA = true

Definition at line 52 of file Main.cpp.

D_1LPT_of_loga

FML::INTERPOLATION::SPLINE::Spline D_1LPT_of_loga

Definition at line 71 of file Main.cpp.

D_2LPT_of_loga

FML::INTERPOLATION::SPLINE::Spline D_2LPT_of_loga

Definition at line 72 of file Main.cpp.

D_3LPTa_of_loga

FML::INTERPOLATION::SPLINE::Spline D_3LPTa_of_loga

Definition at line 73 of file Main.cpp.

D_3LPTb_of_loga

FML::INTERPOLATION::SPLINE::Spline D_3LPTb_of_loga

Definition at line 74 of file Main.cpp.

filename

const std::string filename = "pofk.txt"

Definition at line 63 of file Main.cpp.

fix_amplitude

const bool fix_amplitude = true

Definition at line 44 of file Main.cpp.

fNL

const double fNL = 100.0

Definition at line 46 of file Main.cpp.

FoFBoundarySizeMpch

const double FoFBoundarySizeMpch

Initial value:

=
    3.0

Definition at line 57 of file Main.cpp.

h

const double h = 0.7

Definition at line 66 of file Main.cpp.

interpolation_method

const std::string interpolation_method = "CIC"

Definition at line 64 of file Main.cpp.

kpivot

const double kpivot = 0.05 / h

Definition at line 69 of file Main.cpp.

linking_length

const double linking_length = 0.2

Definition at line 59 of file Main.cpp.

LPT_order_IC

const int LPT_order_IC = 3

Definition at line 53 of file Main.cpp.

n_min_FoF_group

const int n_min_FoF_group = 20

Definition at line 55 of file Main.cpp.

n_s

const double n_s = 0.96

Definition at line 68 of file Main.cpp.

nbin_bispectrum

const int nbin_bispectrum = 16

Definition at line 61 of file Main.cpp.

NDIM

const int NDIM = 3

Definition at line 38 of file Main.cpp.

Ngrid_max_FoF

const int Ngrid_max_FoF = 256

Definition at line 56 of file Main.cpp.

Nmesh

const int Nmesh = 256

Definition at line 48 of file Main.cpp.

Nmesh_IC

const int Nmesh_IC = 128

Definition at line 49 of file Main.cpp.

Npart_1D

const int Npart_1D = 128

Definition at line 50 of file Main.cpp.

nsteps

const int nsteps = 20

Definition at line 43 of file Main.cpp.

OmegaM

const double OmegaM = 1.0

Definition at line 65 of file Main.cpp.

random_seed

const unsigned int random_seed = 1234

Definition at line 47 of file Main.cpp.

type_of_random_field

const std::string type_of_random_field = "gaussian"

Definition at line 45 of file Main.cpp.

zini

const double zini = 20.0

Definition at line 40 of file Main.cpp.