Table of contents


Back to the main page Back to the documentation 1. Overview: N-body 2. Methods 3. COLA simulations 4. Modified gravity simulations

Overview: N-body


We provide some methods for performing some very simple fixed grid N-body simulations (Particle Mesh). The domain decomposition we use in this library is not suitable for huge high-resolution simuations so we only provide these methods and there are some really good and efficient N-body codes like RAMSES and GADGET that one should use for proper N-body anyway.

Methods


In a N-body simulation we solve the equations $$\frac{d^2 x}{d\tau^2} = - \nabla\Phi$$ $$\nabla^2\Phi = C \delta$$ where, for cosmological simulations, $d\tau = H_0 dt / a^2$ and $C = \frac{3}{2}\Omega_M a$ and for Newtonian ones $C = 4\pi G \rho_{\rm mean} / B_{\rm Boxsize}^2$. This can be solved with the Leapfrog method $$x_{n+1} = x_n + v_{n+1/2}\Delta\tau$$ $$v_{n+1/2} = v_{n-1/2} + (\nabla\Phi)_n\Delta\tau$$ The first step is moving the partices (Drift) and the last step is updating the velocities (Kick). The force is computed by performing ${\rm NDIM}+1$ FFTs: $$\nabla\Phi = C F^{-1}[ -\frac{i\vec{k}}{k^2}F[\delta]]$$ The methods below can perform these steps:

// This is solving the Poisson equation D^2phi = norm_poisson_equation * delta and returns DPhi
template <int N, class T>
void compute_force_from_density_real(
                   const FFTWGrid<N> & density_grid_real,
                   std::array<FFTWGrid<N>, N> & force_real,
                   double norm_poisson_equation = 1.0);

// Same as above, but where the density field is already in fourier space
template <int N, class T>
void compute_force_from_density_fourier(
                   const FFTWGrid<N> & density_grid_fourier,
                   std::array<FFTWGrid<N>, N> & force_real,
                   double norm_poisson_equation = 1.0);

// Advance positions from velocity. Takes care of communication in case particles move past domain boundary
template <int N, class T>
void DriftParticles(
                   FML::PARTICLE::MPIParticles <T> & part, 
                   double delta_time, 
                   bool periodic_box = true);

// Does *not* takes care of communication in case particles move past domain boundary
template <int N, class T>
void DriftParticles(
                   T * p,
                   size_t NumPart, 
                   double delta_time, 
                   bool periodic_box = true);

// Advance velocity from force
template <int N, class T>
void KickParticles(std::array<FFTWGrid<N>, N> & force_grid,
                   FML::PARTICLE::MPIParticles<T> & part,
                   double delta_time,
                   std::string interpolation_method);

template <int N, class T>
void KickParticles(std::array<FFTWGrid<N>, N> & force_grid,
                   T * p,
                   size_t NumPart,
                   double delta_time,
                   std::string interpolation_method);

// Do one kick-drift-kick 2nd order leapfrog N-body step
template <int N, class T>
void KickDriftKickNBodyStep(
                   int Nmesh,
                   FML::PARTICLE::MPIParticles<T> & part,
                   double delta_time,
                   std::string density_assignment_method,
                   double norm_poisson_equation)

// Do one 4th order leapfrog N-body step
template <int N, class T>
void YoshidaNBodyStep(
                   int Nmesh,
                   FML::PARTICLE::MPIParticles<T> & part,
                   double delta_time,
                   std::string density_assignment_method,
                   double norm_poisson_equation);

// From a given density field delta(k,zini) generate N-body IC
// using 1, 2 or 3 LPT. 
template <int N, class T>
void NBodyInitialConditions(FML::PARTICLE::MPIParticles<T> & part,
                   int Npart_1D,
                   double buffer_factor,

                   FFTWGrid<N> & delta_fourier,
                   int LPT_order,

                   double box,
                   double zini,
                   std::function<double(double)> & H_over_H0_of_loga,
                   std::vector<std::function<double(double)>> & growth_rate_f_of_loga);

// Generate a random field (gaussian or non-local non-gaussian) of a given type
// and then call the method above
template <int N, class T>
void NBodyInitialConditions(
                   FML::PARTICLE::MPIParticles<T> & part,
                   int Npart_1D,
                   double buffer_factor,

                   int Nmesh,
                   bool fix_amplitude,
                   FML::RANDOM::RandomGenerator * rng,
                   std::function<double(double)> Pofk_of_kBox_over_Pofk_primordal,
                   std::function<double(double)> Pofk_of_kBox_over_volume_primordial,
                   int LPT_order,
                   std::string type_of_random_field,
                   double fNL,

                   double box,
                   double zini,
                   std::function<double(double)> & H_over_H0_of_loga,
                   std::vector<std::function<double(double)>> & growth_rate_f_of_loga);

COLA simulations


In COLA simulations we solve for displacements on top of that predicted by LPT. With this we can do simulations much faster than usual N-body (by sacrificing accuracy on small scales). Instead of propagating $x$ and $v$ via $$\frac{dx}{d\tau} = v$$ $$\frac{d v}{d\tau} = - \nabla\Phi$$ we instead define $x_{\rm COLA} = x - x_{\rm LPT}$ so that the equation becomes $$\frac{dx}{d\tau} = v_{\rm COLA} + \frac{dx_{\rm LPT}}{d\tau}$$ $$\frac{d v_{\rm COLA}}{d\tau} = - \nabla\Phi - \frac{d^2x_{\rm LPT}}{d\tau^2}$$ and we propagate $x$ and $v_{\rm COLA}$ instead. Note that in this frame the initial velocity is $v_{\rm COLA} = 0$ and this helps us keep in line with linear theory on the largest scales even with large time-steps. Also note that $x_{\rm LPT} \sim D_1(\tau)\Psi_1(q,\tau_{\rm ini}) + D_2(\tau)\Psi_2(q,\tau_{\rm ini}) + \ldots$ so the derivatives of $x_{\rm LPT}$ above just becomes some growth-factor times the initial displacement field at particle positions when we start the simulation. The library is set up such that when we create initial conditions for N-body then we store displacement fields at the particle positions if they have the related get_D_1LPT (and 2LPT, 3LPTa, 3LPTb if you want to use higher order) methods. With this its easy to turn a normal N-body into a COLA simulations - all we need is to compute some ODEs to get the growth-factors and add in the additional terms in the equations of motion (see the N-body example for this).

Modified gravity simulations


The library contains some methods for computing the additional fifth-force in modified gravity models. We can compute the fifth-force in any linear model exactly, we can estimate the fifth-force in models with a screening mechanism using an approximate method and if you want you can always use the multigrid solver to solve the exact equation.

The methods below compute the fifth-force in terms of a fifth-force "density contrast" that is such that the total force (in fourier space) is just $$F[\nabla \Phi_{\rm tot}] = F[\nabla \Phi + \nabla \Phi_{\rm fifth-force}] \propto \frac{i\vec{k}}{k^2}(\delta(k) + \delta_{\rm fifth-force})$$ i.e. after computing this we only need to add this to the density contrast we would use to compute the usual gravitational force in order to get the total-force.

// For models where the fifth-force is given as F = Dphi where 
// F[D^2 phi] = m^2 F[phi] + f(k)F[delta]
// Here the coupling_factor is f(k)k^2/(k^2 + m^2)
// e.g. for linearized f(R) the coupling factor is 1/3 * k^2/(k^2 + m^2)
// and for linearized DGP the coupling factor is 1/(3betaDGP)
template <int N>
  void compute_delta_fifth_force(
                   const FFTWGrid<N> & density_fourier,
                   FFTWGrid<N> & density_mg_fourier,
                   std::function<double(double)> coupling_factor_of_kBox);

// Approximate method of Winther and Ferreira 2015
// For MG models that has a potential dependent screening (f(R) and the alike)
// D^2 phi = m^2 phi + f(Phi)delta
// e.g. for f(R) the screening factor is f(Phi) = 1.5 |f_R / Phi|
// and the coupling factor is 1/3 * k^2 / (k^2 + m^2)
// poisson_norm is used to convert delta to Phi and is defined via D^2 Phi = poisson_norm * delta
template <int N>
  void compute_delta_fifth_force_potential_screening(
                   const FFTWGrid<N> & density_fourier,
                   FFTWGrid<N> & density_mg_fourier,
                   std::function<double(double)> coupling_factor_of_kBox);
                   std::function<double(double)> screening_factor_of_newtonian_potential,
                   double poisson_norm);

// Approximate method of Winther and Ferreira 2015
// For MG models that has a density dependent screening (DGP and the alike)
// D^2 phi = m^2 phi + f(rho)delta
// e.g. for DGP m=0 and the screening factor f(rho) = 2(sqrt(1 + C) - 1)/C where C = const * rho 
// and coupling_factor = 1/(3betaDGP)
template <int N>
  void compute_delta_fifth_force_density_screening(
                   const FFTWGrid<N> & density_fourier,
                   FFTWGrid<N> & density_mg_fourier,
                   std::function<double(double)> coupling_factor_of_kBox);
                   std::function<double(double)> screening_factor_of_density,
                   double smoothing_scale,
                   std::string smoothing_method);