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);