This namespace deals with various cosmology specific things: background evolution, perturbation theory, recombination history etc. More...

Namespaces
namespace	HALOMODEL

namespace	LPT
	This namespace contains things related to Lagrangian Perturbation Theory (displacement fields, reconstruction, initial conditions etc.)

namespace	SPHERICALCOLLAPSE

Data Structures
struct	_PerturbationSystemInfo
	Struct for keeping track over what quantities we have at what index when solving the linear Einstein-Boltzmann system. More...

class	BackgroundCosmology
	Computing the background evolution of our Universe (LCDM). More...

class	Perturbations
	Class for solving the linear perturbations (LCDM) in temperature, baryon, CDM, massless neutrinos etc. More...

class	PowerSpectrum
	Computing power-spectra (matter, CMB, correlation functions etc.) in linear perturbation theory. More...

class	RecombinationHistory
	For computing the recombination history of our Universe. More...

Typedefs
using	ParameterMap = FML::UTILS::ParameterMap

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

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

using	Spline = FML::INTERPOLATION::SPLINE::Spline

using	DVector = std::vector< double >

using	BackgroundCosmology = FML::COSMOLOGY::BackgroundCosmology

using	RecombinationHistory = FML::COSMOLOGY::RecombinationHistory

using	ODEFunctionJacobian = FML::SOLVERS::ODESOLVER::ODEFunctionJacobian

using	Spline2D = FML::INTERPOLATION::SPLINE::Spline2D

using	DVector2D = std::vector< DVector >

typedef struct FML::COSMOLOGY::_PerturbationSystemInfo	PerturbationSystemInfo
	Struct for keeping track over what quantities we have at what index when solving the linear Einstein-Boltzmann system. More...

using	Perturbations = FML::COSMOLOGY::Perturbations

Functions
double	correlation_function_single (double r, std::function< double(double)> &Delta_P, double kmin, double kmax)

std::pair< DVector, DVector >	correlation_function_single_fftlog (std::function< double(double)> &Delta_P, double rmin, double rmax, int ngrid)

std::pair< DVector, DVector >	correlation_function_single_fftw (std::function< double(double)> &Delta_P, double rmin, double rmax, int ngrid_min)

template<class T >
void	particles_to_redshiftspace (FML::PARTICLE::MPIParticles< T > &part, std::vector< double > line_of_sight_direction, double velocity_to_displacement)
	This takes a set of particles and displace them from realspace to redshiftspace Using a fixed line of sight direction DeltaX = (v * r)r * velocity_to_displacement If velocities are peculiar then velocity_to_displacement = 1/(100 * aH/H0 * box) with box in Mpc/h You can revert back to real-space by calling the method again with velocity_to_displacement being negative. More...

Variables
FML::UTILS::ConstantsAndUnits	Constants ("SI")

Detailed Description

This namespace deals with various cosmology specific things: background evolution, perturbation theory, recombination history etc.

Typedef Documentation

◆ BackgroundCosmology

typedef FML::COSMOLOGY::BackgroundCosmology FML::COSMOLOGY::BackgroundCosmology

Definition at line 31 of file Perturbations.h.

◆ DVector

typedef std::vector< double > FML::COSMOLOGY::DVector

Definition at line 37 of file BackgroundCosmology.h.

◆ DVector2D

typedef std::vector< DVector > FML::COSMOLOGY::DVector2D

Definition at line 39 of file Perturbations.h.

◆ ODEFunction

typedef FML::SOLVERS::ODESOLVER::ODEFunction FML::COSMOLOGY::ODEFunction

Definition at line 34 of file BackgroundCosmology.h.

◆ ODEFunctionJacobian

typedef FML::SOLVERS::ODESOLVER::ODEFunctionJacobian FML::COSMOLOGY::ODEFunctionJacobian

Definition at line 35 of file Perturbations.h.

◆ ODESolver

typedef FML::SOLVERS::ODESOLVER::ODESolver FML::COSMOLOGY::ODESolver

Definition at line 35 of file BackgroundCosmology.h.

◆ ParameterMap

typedef FML::UTILS::ParameterMap FML::COSMOLOGY::ParameterMap

Definition at line 33 of file BackgroundCosmology.h.

◆ Perturbations

using FML::COSMOLOGY::Perturbations = typedef FML::COSMOLOGY::Perturbations

Definition at line 43 of file PowerSpectrum.h.

◆ PerturbationSystemInfo

typedef struct FML::COSMOLOGY::_PerturbationSystemInfo FML::COSMOLOGY::PerturbationSystemInfo

Struct for keeping track over what quantities we have at what index when solving the linear Einstein-Boltzmann system.

◆ RecombinationHistory

typedef FML::COSMOLOGY::RecombinationHistory FML::COSMOLOGY::RecombinationHistory

Definition at line 32 of file Perturbations.h.

◆ Spline

typedef FML::INTERPOLATION::SPLINE::Spline FML::COSMOLOGY::Spline

Definition at line 36 of file BackgroundCosmology.h.

◆ Spline2D

typedef FML::INTERPOLATION::SPLINE::Spline2D FML::COSMOLOGY::Spline2D

Definition at line 37 of file Perturbations.h.

Function Documentation

◆ correlation_function_single()

double FML::COSMOLOGY::correlation_function_single	(	double	r,
		std::function< double(double)> &	Delta_P,
		double	kmin,
		double	kmax
	)

Definition at line 920 of file PowerSpectrum.cpp.

                                                                                                               {
 
            ODEFunction dxidlogk_func = [&](double logkr, [[maybe_unused]] const double * y, double * dxidlogkr) {
                double kr = std::exp(logkr);
                dxidlogkr[0] = std::sin(kr) / (kr)*Delta_P(kr / r);
                return GSL_SUCCESS;
            };
 
            DVector xi_ini{0.0};
            DVector range{std::log(kmin * r), std::log(kmax * r)};
            ODESolver xi_ode(1e-5, 1e-9, 1e-9);
            xi_ode.solve(dxidlogk_func, range, xi_ini);
 
            return xi_ode.get_final_data_by_component(0);
        }

◆ correlation_function_single_fftlog()

std::pair< DVector, DVector > FML::COSMOLOGY::correlation_function_single_fftlog	(	std::function< double(double)> &	Delta_P,
		double	rmin,
		double	rmax,
		int	ngrid
	)

Definition at line 938 of file PowerSpectrum.cpp.

                                                                                  {
            // We want xi(r) for rmin ... rmax. We solve for it in rmin/padding ... rmax * padding
            // to avoid ringing/aliasing close to the edges. A factor of 15 seems to be enough
            // The number of points sets the accuracy. The fiducial choice below is to ensure
            // 1% accuracy around the BAO peak
            const double r0 = std::sqrt(rmin * rmax);
            const double paddingfactor = 15.0;
            const double k0 = 1.0 / r0;
 
            // Set ranges and make a log-spaced k array
            const double kmin_fft = k0 * r0 / rmax / paddingfactor;
            const double kmax_fft = k0 * r0 / rmin * paddingfactor;
            DVector k_array = FML::MATH::linspace(std::log(kmin_fft), std::log(kmax_fft), ngrid);
            for (auto & k : k_array)
                k = std::exp(k);
 
            // Fill P(k) array
            DVector Pk_array(ngrid);
            for (int i = 0; i < ngrid; i++) {
                Pk_array[i] = Delta_P(k_array[i]) / std::pow(k_array[i], 3) * (2.0 * M_PI * M_PI);
            }
 
            // FFTLog algorithm
            auto res = FML::SOLVERS::FFTLog::ComputeCorrelationFunction(k_array, Pk_array);
 
            return res;
        }

◆ correlation_function_single_fftw()

std::pair< DVector, DVector > FML::COSMOLOGY::correlation_function_single_fftw	(	std::function< double(double)> &	Delta_P,
		double	rmin,
		double	rmax,
		int	ngrid_min
	)

Definition at line 970 of file PowerSpectrum.cpp.

                                                                                    {
 
            // Set the boxsize and the grid resolution
            double Box = 4.0 * rmax;
            int ngrid = std::max(int(4.0 * (Box / rmin)), ngrid_min);
 
            // Find closest power of two for optimal FFTW
            for (int N = 2;; N *= 2) {
                if (N >= ngrid or N > 1e9) {
                    ngrid = N;
                    break;
                }
            }
 
            // Set up grid in k-space and fill it with the source
            fftw_complex * source = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * ngrid);
            for (int i = 0; i < ngrid; i++) {
                double im = 0.0;
                if (i == 0 or i == ngrid / 2) {
                    // Zero mode and largest mode
                    im = 0.0;
                } else if (i < ngrid / 2) {
                    // Positive modes
                    double k = 2.0 * M_PI / Box * i;
                    im = -Delta_P(k) / (k * k) / 2.0;
                } else {
                    // Negative modes
                    double k = 2.0 * M_PI / Box * (i - ngrid);
                    im = +Delta_P(-k) / (k * k) / 2.0;
                }
                source[i][0] = 0.0;
                source[i][1] = im;
            }
 
            // Transform to real space to get r*xi(r)
            fftw_plan plan = fftw_plan_dft_1d(ngrid, source, source, FFTW_BACKWARD, FFTW_MEASURE);
            fftw_execute(plan);
            fftw_destroy_plan(plan);
 
            // Normalize to get xi (i = 0 omitted to avoid division by 0)
            for (int i = 0; i < ngrid; i++) {
                double r = i * Box / double(ngrid);
                source[i][0] *= 2.0 * M_PI / Box;
                if (i > 0)
                    source[i][0] /= r;
            }
 
            // Make a spline (we only keep the
            DVector r_arr;
            DVector xi_arr;
            r_arr.reserve(ngrid);
            xi_arr.reserve(ngrid);
            for (int i = 0; i < ngrid; i++) {
                double r = i * Box / double(ngrid);
                if (r >= rmin and r <= rmax) {
                    r_arr.push_back(r);
                    xi_arr.push_back(source[i][0]);
                }
            }
 
            // Clean up FFTW
            fftw_free(source);
 
            return {r_arr, xi_arr};
        }

◆ particles_to_redshiftspace()

template<class T >

void FML::COSMOLOGY::particles_to_redshiftspace	(	FML::PARTICLE::MPIParticles< T > &	part,
		std::vector< double >	line_of_sight_direction,
		double	velocity_to_displacement
	)

This takes a set of particles and displace them from realspace to redshiftspace Using a fixed line of sight direction DeltaX = (v * r)r * velocity_to_displacement If velocities are peculiar then velocity_to_displacement = 1/(100 * aH/H0 * box) with box in Mpc/h You can revert back to real-space by calling the method again with velocity_to_displacement being negative.

Template Parameters

T	The particle class

Parameters

[out]	part	MPIParticle container
[in]	line_of_sight_direction	The fixed line of sight vector, e.g. (0,0,1) for the z-axis
[in]	velocity_to_displacement	Convert from user velocity to RSD shift.

Definition at line 244 of file Reconstruction.h.

                                                                         {
 
            // Fetch how many dimensjons we are working in
            const int N = FML::PARTICLE::GetNDIM(T());
 
            // Check that velocities really exists (i.e. get_vel is not just set to return a nullptr)
            static_assert(FML::PARTICLE::has_get_pos<T>(),
                          "[particles_to_redshiftspace] Partices must have positions to use this method");
            static_assert(FML::PARTICLE::has_get_vel<T>(),
                          "[particles_to_redshiftspace] Partices must have velocity to use this method");
 
            // Periodic box? Yes, this is only meant to be used with simulation boxes
            const bool periodic_box = true;
 
            // Make sure line_of_sight_direction is a unit vector
            double norm = 0.0;
            for (int idim = 0; idim < N; idim++) {
                norm += line_of_sight_direction[idim] * line_of_sight_direction[idim];
            }
            norm = std::sqrt(norm);
            assert_mpi(norm > 0.0, "[particles_to_redshiftspace] Line of sight vector cannot be the zero vector\n");
            for (int idim = 0; idim < N; idim++) {
                line_of_sight_direction[idim] /= norm;
            }
 
            double max_disp = 0.0;
            auto NumPart = part.get_npart();
            auto * p = part.get_particles_ptr();
#ifdef USE_OMP
#pragma omp parallel for reduction(max : max_disp)
#endif
            for (size_t i = 0; i < NumPart; i++) {
                auto * pos = FML::PARTICLE::GetPos(p[i]);
                auto * vel = FML::PARTICLE::GetVel(p[i]);
                double vdotr = 0.0;
                for (int idim = 0; idim < N; idim++) {
                    vdotr += vel[idim] * line_of_sight_direction[idim];
                }
                max_disp = std::max(max_disp, std::fabs(vdotr));
                for (int idim = 0; idim < N; idim++) {
                    pos[idim] += vdotr * line_of_sight_direction[idim] * velocity_to_displacement;
                    // Periodic boundary conditions
                    if (periodic_box) {
                        if (pos[idim] < 0.0)
                            pos[idim] += 1.0;
                        if (pos[idim] >= 1.0)
                            pos[idim] -= 1.0;
                    }
                }
            }
 
            // Only need to communicate if we have shifted along the x-axis
            if (line_of_sight_direction[0] != 0.0)
                part.communicate_particles();
 
            // If velocity_to_displacement < 0 we are going the other way around
            FML::MaxOverTasks(&max_disp);
            max_disp *= velocity_to_displacement;
            if (FML::ThisTask == 0)
                std::cout << "[particles_to_redshiftspace] Maximum displacement "
                          << (velocity_to_displacement < 0 ? "(reversing) : " : "            :  ") << max_disp << "\n";
        }

Variable Documentation

◆ Constants

FML::UTILS::ConstantsAndUnits FML::COSMOLOGY::Constants ( "SI" )

Definition at line 40 of file BackgroundCosmology.h.

Namespaces

Data Structures

Typedefs

Functions

Variables