FML
Typedefs | Functions
FML::NBODY Namespace Reference

This namespace deals with N-body simulations. More...

Typedefs

template<int N>
using FFTWGrid = FML::GRID::FFTWGrid< N >
 
template<class T >
using MPIParticles = FML::PARTICLE::MPIParticles< T >
 

Functions

template<int N, class T >
void DriftParticles (FML::PARTICLE::MPIParticles< T > &part, double delta_time, bool periodic_box=true)
 
template<int N, class T >
void DriftParticles (T *p, size_t NumPart, double delta_time, bool periodic_box)
 This moves the particles according to \( x_{\rm new} = x + v \Delta t \). More...
 
template<int N, class T >
void KickParticles (std::array< FFTWGrid< N >, N > &force_grid, MPIParticles< T > &part, double delta_time, std::string interpolation_method)
 This moves the particle velocities according to \( v_{\rm new} = v + F \Delta t \). More...
 
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)
 This moves the particle velocities according to \( v_{\rm new} = v + F \Delta t \). More...
 
template<int N>
void compute_force_from_density_real (const FFTWGrid< N > &density_grid_real, std::array< FFTWGrid< N >, N > &force_real, std::string density_assignment_method_used, double norm_poisson_equation)
 Take a density grid in real space and returns the force \( \nabla \phi \) where \( \nabla^2 \phi = {\rm norm} \cdot \delta \) Different choices for what kernel to use for \( \nabla / \nabla^2\) are availiable, see the function body (is set too be a compile time option). More...
 
template<int N>
void compute_force_from_density_fourier (const FFTWGrid< N > &density_grid_fourier, std::array< FFTWGrid< N >, N > &force_real, std::string density_assignment_method_used, double norm_poisson_equation)
 Take a density grid in fourier space and returns the force \( \nabla \phi \) where \( \nabla^2 \phi = {\rm norm} \cdot \delta \) Different choices for what kernel to use for \( \nabla / \nabla^2\) are availiable, see the function body (is set too be a compile time option). More...
 
template<int N, class T >
void KickDriftKickNBodyStep (int Nmesh, MPIParticles< T > &part, double delta_time, std::string density_assignment_method, double norm_poisson_equation)
 Take a N-body step with a simple Kick-Drift-Kick method (this method serves mainly as an example for how one can do this). More...
 
template<int N, class T >
void YoshidaNBodyStep (int Nmesh, MPIParticles< T > &part, double delta_time, std::string density_assignment_method, double norm_poisson_equation)
 Take a N-body step with a 4th order symplectic Yoshida method. More...
 
template<int N, class T >
void DriftParticles (MPIParticles< T > &part, double delta_time, bool periodic_box)
 This moves the particles according to \( x_{\rm new} = x + v \Delta t \). More...
 
template<int N, class T >
void NBodyInitialConditions (MPIParticles< T > &part, int Npart_1D, double buffer_factor, const FFTWGrid< N > &delta_fourier, std::vector< FFTWGrid< N > > &phi_nLPT_potentials, int LPT_order, double box, double zini, std::vector< double > velocity_norms)
 Generate particles from a given initial density field using Lagrangian perturbation theory. More...
 
template<int N, class T >
void NBodyInitialConditions (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::vector< double > velocity_norms)
 Generate particles from a given power-spectrum using Lagrangian perturbation theory. More...
 
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)
 This method computes the fifth-force potential for modified gravity models using the linear approximation This computes \( \delta_{\rm MG}(k) \) where the total force in fourier space is \( F(k) \propto \frac{\vec{k}}{k^2}[\delta(k) + \delta_{\rm MG}(k)] \) by solving \( \nabla^2 \phi = m^2 \phi + F^{-1}[g(k) \delta(k)] \) where \( \delta_{\rm MG}(k) = -k^2\phi(k) \). More...
 
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)
 This method computes the fifth-force potential for modified gravity models which has a screening mechanism using the approximate method of Winther & Ferreira 2015. More...
 
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)
 This method computes the fifth-force potential for modified gravity models which has a screening mechanism using the approximate method of Winther & Ferreira 2015. More...
 
double compute_sigma_of_R (std::function< double(double)> Pofk_of_kBox_over_volume, double R_mpch, double boxsize_mpch)
 This computes the standard deviation of linear fluctuations smoothed over a sphere of radius, \( \sigma(R) \), by integrating the power-spectrum convolved with a top-hat windowfunction of radius \( R \). More...
 

Detailed Description

This namespace deals with N-body simulations.

Computing forces, moving particles and generating initial conditions.

Typedef Documentation

◆ FFTWGrid

template<int N>
using FML::NBODY::FFTWGrid = typedef FML::GRID::FFTWGrid<N>

Definition at line 29 of file NBody.h.

◆ MPIParticles

template<class T >
using FML::NBODY::MPIParticles = typedef FML::PARTICLE::MPIParticles<T>

Definition at line 31 of file NBody.h.

Function Documentation

◆ compute_delta_fifth_force()

template<int N>
void FML::NBODY::compute_delta_fifth_force ( const FFTWGrid< N > &  density_fourier,
FFTWGrid< N > &  density_mg_fourier,
std::function< double(double)>  coupling_factor_of_kBox 
)

This method computes the fifth-force potential for modified gravity models using the linear approximation This computes \( \delta_{\rm MG}(k) \) where the total force in fourier space is \( F(k) \propto \frac{\vec{k}}{k^2}[\delta(k) + \delta_{\rm MG}(k)] \) by solving \( \nabla^2 \phi = m^2 \phi + F^{-1}[g(k) \delta(k)] \) where \( \delta_{\rm MG}(k) = -k^2\phi(k) \).

For example in \( f(R) \) gravity we have \( g(k) = \frac{1}{3}\frac{k^2}{k^2 + m^2}\) and in DGP we have \( g(k) = \frac{1}{3\beta} \) (independent of scale).

Template Parameters
NThe dimension we work in.
Parameters
[in]density_fourierThe density contrast in fourier space.
[out]density_mg_fourierThe force potential.
[in]coupling_factor_of_kBoxThe coupling factor \( g(k) \)

Definition at line 1072 of file NBody.h.

1074 {
1075
1076 auto coupling_factor_of_kBox_spline =
1077 density_fourier.make_fourier_spline(coupling_factor_of_kBox, "MG coupling(k)");
1078 const auto Local_nx = density_fourier.get_local_nx();
1079 density_mg_fourier = density_fourier;
1080#ifdef USE_OMP
1081#pragma omp parallel for
1082#endif
1083 for (int islice = 0; islice < Local_nx; islice++) {
1084 [[maybe_unused]] double kmag;
1085 [[maybe_unused]] std::array<double, N> kvec;
1086 for (auto && fourier_index : density_mg_fourier.get_fourier_range(islice, islice + 1)) {
1087 // Get wavevector and magnitude
1088 density_mg_fourier.get_fourier_wavevector_and_norm_by_index(fourier_index, kvec, kmag);
1089
1090 // Compute coupling
1091 auto coupling = coupling_factor_of_kBox_spline(kmag);
1092
1093 // Multiply by coupling
1094 auto value = density_mg_fourier.get_fourier_from_index(fourier_index);
1095 density_mg_fourier.set_fourier_from_index(fourier_index, value * FML::GRID::FloatType(coupling));
1096 }
1097 }
1098 }
float FloatType
Definition: FFTWGlobal.h:29
ptrdiff_t get_local_nx() const
Number of local x-slices in the real grid.
Definition: FFTWGrid.h:1079
void set_fourier_from_index(const IndexIntType ind, const ComplexType value)
Set value of cell in fourier grid using (local) index.
Definition: FFTWGrid.h:1036
ComplexType get_fourier_from_index(const IndexIntType index) const
Fetch value in fourier grid by (local) index.
Definition: FFTWGrid.h:1031
FourierRange get_fourier_range(int islice_begin=0, int islice_end=0) const
Range iterator for going through all active cells in the main fourier grid by index [ e....
Definition: FFTWGrid.h:419
void get_fourier_wavevector_and_norm_by_index(const IndexIntType ind, std::array< double, N > &kvec, double &kmag) const
From index in the grid get the k-vector and the magnitude of it.
Definition: FFTWGrid.h:1103
FML::INTERPOLATION::SPLINE::Spline make_fourier_spline(std::function< double(double)> function_of_kBox, std::string label, int sampling_factor=4) const
std::function can be slow so for looping through a fourier grid and evaluating a function f(k) in eve...
Definition: FFTWGrid.h:1366

◆ compute_delta_fifth_force_density_screening()

template<int N>
void FML::NBODY::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 
)

This method computes the fifth-force potential for modified gravity models which has a screening mechanism using the approximate method of Winther & Ferreira 2015.

This computes \( \delta_{\rm MG}(k) \) where the total force in fourier space is \( F(k) \propto \frac{\vec{k}}{k^2}[\delta(k) + \delta_{\rm MG}(k)] \) by solving \( \nabla^2 \phi = m^2 \phi + f(\rho)F^{-1}[g(k) \delta(k)] \) where \( \delta_{\rm MG}(k) = -k^2\phi(k) \) and \( \rho \) is the density in units of the mean density. For example in DGP gravity we have \( g(k) = \frac{1}{3\beta} \) (independent of scale) and the screening function is \( f(\rho) = 2\frac{\sqrt{1 + C} - 1}{C} \) where \( C = \frac{8\Omega_M(r_cH_0)^2}{9\beta^2}\rho \) If you don't want screening then simply pass the function \( f \equiv 1 \) and the equation reduces to the one in the linear regime.

Template Parameters
NThe dimension we work in.
Parameters
[in]density_fourierThe density contrast in fourier space.
[out]density_mg_fourierThe force potential.
[in]coupling_factor_of_kBoxThe coupling factor \( g(k) \)
[in]screening_factor_of_densityThe screening factor \( f(\rho) \) Should be in \( [0,1] \) and go to 1 for \( \rho \to 0 \) and to 0 for \( \rho \to \infty \).
[in]smoothing_scaleThe smoothing radius in units of the boxsize.
[in]smoothing_methodThe k-space smoothing filter (gaussian, tophat, sharpk).

Definition at line 1224 of file NBody.h.

1229 {
1230
1231 const auto Local_nx = density_fourier.get_local_nx();
1232
1233 // Copy over
1234 density_mg_fourier = density_fourier;
1235
1236 // Smooth density field
1237 if (smoothing_scale > 0.0)
1238 FML::GRID::smoothing_filter_fourier_space(density_mg_fourier, smoothing_scale, smoothing_method);
1239
1240 // To real space density field
1241 density_mg_fourier.fftw_c2r();
1242
1243 // Apply screening function
1244#ifdef USE_OMP
1245#pragma omp parallel for
1246#endif
1247 for (int islice = 0; islice < Local_nx; islice++) {
1248 for (auto && real_index : density_mg_fourier.get_real_range(islice, islice + 1)) {
1249 auto delta = density_mg_fourier.get_real_from_index(real_index);
1250 auto screening_factor = screening_factor_of_density(1.0 + delta);
1251 density_mg_fourier.set_real_from_index(real_index, delta * screening_factor);
1252 }
1253 }
1254
1255 // To fourier spacce
1256 density_mg_fourier.fftw_r2c();
1257
1258 // Apply coupling
1259#ifdef USE_OMP
1260#pragma omp parallel for
1261#endif
1262 for (int islice = 0; islice < Local_nx; islice++) {
1263 [[maybe_unused]] double kmag;
1264 [[maybe_unused]] std::array<double, N> kvec;
1265 for (auto && fourier_index : density_mg_fourier.get_fourier_range(islice, islice + 1)) {
1266 auto value = density_mg_fourier.get_fourier_from_index(fourier_index);
1267
1268 // Get wavevector and magnitude
1269 density_mg_fourier.get_fourier_wavevector_and_norm_by_index(fourier_index, kvec, kmag);
1270
1271 // Compute coupling
1272 auto coupling = coupling_factor_of_kBox(kmag);
1273
1274 // Multiply by coupling
1275 density_mg_fourier.set_fourier_from_index(fourier_index, value * FML::GRID::FloatType(coupling));
1276 }
1277 }
1278 }
void fftw_c2r()
Perform complex-to-real fourier transform.
Definition: FFTWGrid.h:946
void set_real_from_index(const IndexIntType ind, const FloatType value)
Set value in grid from (local) index.
Definition: FFTWGrid.h:1020
FloatType get_real_from_index(const IndexIntType index) const
Fetch value in grid by (local) index.
Definition: FFTWGrid.h:1253
void fftw_r2c()
Perform real-to-complex fourier transform.
Definition: FFTWGrid.h:881
RealRange get_real_range(int islice_begin=0, int islice_end=0) const
Range iterator for going through all active cells in the main real real grid by index [ e....
Definition: FFTWGrid.h:400
void smoothing_filter_fourier_space(FFTWGrid< N > &fourier_grid, double smoothing_scale, std::string smoothing_method)
Low-pass filters (tophat, gaussian, sharpk)

◆ compute_delta_fifth_force_potential_screening()

template<int N>
void FML::NBODY::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 
)

This method computes the fifth-force potential for modified gravity models which has a screening mechanism using the approximate method of Winther & Ferreira 2015.

This computes \( \delta_{\rm MG}(k) \) where the total force in fourier is given by \( F(k) \propto \frac{\vec{k}}{k^2}[\delta(k) + \delta_{\rm MG}(k)] \) by solving \( \nabla^2 \phi = m^2 \phi + f(\Phi)F^{-1}[g(k) \delta(k)] \) where \( \delta_{\rm MG}(k) = -k^2\phi(k) \) For example in \( f(R) \) gravity we have \( g(k) = \frac{1}{3}\frac{k^2}{k^2 + m^2} \) and the screening function is \( f(\Phi) = \min(1, \left|\frac{3f_R}{2\Phi}\right|) \) If you don't want screening then simpy pass the function \( f \equiv 1 \) and the equation reduces to the one in the linear regime

Template Parameters
NThe dimension we work in.
Parameters
[in]density_fourierThe density contrast in fourier space.
[out]density_mg_fourierThe force potential.
[in]coupling_factor_of_kBoxThe coupling factor \(g(k)\)
[in]screening_factor_of_newtonian_potentialThe screening factor \( f(\Phi_N) \) Should be in \( [0,1] \) and go to 1 for \( \Phi_N \to 0 \) and 0 for very large \( \Phi_N \)
[in]poisson_normThe factor \( C \) in \( \nabla^2\Phi = C\delta \) to get the potential in the metric (not the code-potential) so \( C = \frac{3}{2}\Omega_M a \frac{(H_0B)^2}{a^2} \)

Definition at line 1122 of file NBody.h.

1127 {
1128
1129 const auto Local_nx = density_fourier.get_local_nx();
1130 const auto Local_x_start = density_fourier.get_local_x_start();
1131
1132 // Make copy of density grid
1133 density_mg_fourier = density_fourier;
1134
1135 // Transform to Newtonian potential
1136#ifdef USE_OMP
1137#pragma omp parallel for
1138#endif
1139 for (int islice = 0; islice < Local_nx; islice++) {
1140 [[maybe_unused]] double kmag2;
1141 [[maybe_unused]] std::array<double, N> kvec;
1142 for (auto && fourier_index : density_mg_fourier.get_fourier_range(islice, islice + 1)) {
1143
1144 auto value = density_mg_fourier.get_fourier_from_index(fourier_index);
1145 density_mg_fourier.get_fourier_wavevector_and_norm2_by_index(fourier_index, kvec, kmag2);
1146 value *= -poisson_norm / kmag2;
1147
1148 density_mg_fourier.set_fourier_from_index(fourier_index, value);
1149 }
1150 }
1151
1152 // Set DC mode
1153 if (Local_x_start == 0)
1154 density_mg_fourier.set_fourier_from_index(0, 0.0);
1155
1156 // Take another copy of the density field as we need it in real space
1157 auto delta_real = density_fourier;
1158
1159 // Transform to real space: Phi(x) and delta(x)
1160 density_mg_fourier.fftw_c2r();
1161 delta_real.fftw_c2r();
1162
1163 // Apply screening function
1164#ifdef USE_OMP
1165#pragma omp parallel for
1166#endif
1167 for (int islice = 0; islice < Local_nx; islice++) {
1168 for (auto && real_index : density_mg_fourier.get_real_range(islice, islice + 1)) {
1169 auto phi_newton = density_mg_fourier.get_real_from_index(real_index);
1170 auto delta = delta_real.get_real_from_index(real_index);
1171 auto screening_factor = screening_factor_of_newtonian_potential(phi_newton);
1172 density_mg_fourier.set_real_from_index(real_index, delta * screening_factor);
1173 }
1174 }
1175
1176 // Back to fourier space: delta(k)
1177 density_mg_fourier.fftw_r2c();
1178
1179 // Apply coupling
1180#ifdef USE_OMP
1181#pragma omp parallel for
1182#endif
1183 for (int islice = 0; islice < Local_nx; islice++) {
1184 [[maybe_unused]] double kmag;
1185 [[maybe_unused]] std::array<double, N> kvec;
1186 for (auto && fourier_index : density_mg_fourier.get_fourier_range(islice, islice + 1)) {
1187 auto value = density_mg_fourier.get_fourier_from_index(fourier_index);
1188
1189 // Get wavevector and magnitude
1190 density_mg_fourier.get_fourier_wavevector_and_norm_by_index(fourier_index, kvec, kmag);
1191
1192 // Compute coupling
1193 auto coupling = coupling_factor_of_kBox(kmag);
1194
1195 // Multiply by coupling
1196 density_mg_fourier.set_fourier_from_index(fourier_index, value * FML::GRID::FloatType(coupling));
1197 }
1198 }
1199 }
ptrdiff_t get_local_x_start() const
The x-slice the first local slice has in the global grid.
Definition: FFTWGrid.h:1084
void get_fourier_wavevector_and_norm2_by_index(const IndexIntType ind, std::array< double, N > &kvec, double &kmag2) const
From index in the grid get the k-vector and the norm of it (square magnitude)
Definition: FFTWGrid.h:1122

◆ compute_force_from_density_fourier()

template<int N>
void FML::NBODY::compute_force_from_density_fourier ( const FFTWGrid< N > &  density_grid_fourier,
std::array< FFTWGrid< N >, N > &  force_real,
std::string  density_assignment_method_used,
double  norm_poisson_equation 
)

Take a density grid in fourier space and returns the force \( \nabla \phi \) where \( \nabla^2 \phi = {\rm norm} \cdot \delta \) Different choices for what kernel to use for \( \nabla / \nabla^2\) are availiable, see the function body (is set too be a compile time option).

Fiducial choice is the continuous greens function \( 1/k^2\), but we can also choose to also devonvolve the window and discrete kernels (Hamming 1989; same as used in GADGET) and Hockney & Eastwood 1988. See e.g. 1603.00476 for a list and references.

Template Parameters
NThe dimension of the grid
Parameters
[in]density_grid_fourierThe density contrast in fourier space.
[out]force_realThe force in real space.
[in]density_assignment_method_usedThe density assignement we used to compute the density field. Needed only in case kernel_choice (defined in the body of this function) is not CONTINUOUS_GREENS_FUNCTION.
[in]norm_poisson_equationThe prefactor (norm) to the Poisson equation.

Definition at line 253 of file NBody.h.

256 {
257
258 // What fourier space kernel to use for D/D^2
259 enum KernelChoices {
260 // 1/k^2
261 CONTINUOUS_GREENS_FUNCTION,
262 // Divide by square of density assignment window function 1/k^2W^2
263 CONTINUOUS_GREENS_FUNCTION_DECONVOLVE,
264 // Hockney & Eastwood 1988: 1 / [ 4/dx^2 * Sum sin(ki * dx / 2)^2 ] with dx = 1/Ngrid
265 DISCRETE_GREENS_FUNCTION_HOCKNEYEASTWOOD,
266 // Hockney & Eastwood 1988: 1 / [ 4/dx^2 * Sum sin(ki * dx / 2)^2 ] / W^2 with dx = 1/Ngrid
267 DISCRETE_GREENS_FUNCTION_HOCKNEYEASTWOOD_DECONVOLVE,
268 // Hamming: D = D/k^2 where D = (8 sin(k) - sin(2k))/6
269 DISCRETE_GREENS_FUNCTION_HAMMING,
270 // Hamming: D/k^2W^2 where D = (8 sin(k) - sin(2k))/6 (GADGET2 kernel)
271 DISCRETE_GREENS_FUNCTION_HAMMING_DECONVOLVE
272 };
273 constexpr int kernel_choice = CONTINUOUS_GREENS_FUNCTION;
274
275 auto Nmesh = density_grid_fourier.get_nmesh();
276 auto Local_nx = density_grid_fourier.get_local_nx();
277 auto Local_x_start = density_grid_fourier.get_local_x_start();
278
279 // This is needed in case kernel_choice != CONTINUOUS_GREENS_FUNCTION
280 // The order of the density assignment method
281 const int order = FML::INTERPOLATION::interpolation_order_from_name(density_assignment_method_used);
282 // Window function for density assignment
283 const double knyquist = M_PI * Nmesh;
284 [[maybe_unused]] auto window_function = [&](std::array<double, N> & kvec) -> double {
285 double w = 1.0;
286 for (int idim = 0; idim < N; idim++) {
287 const double koverkny = M_PI / 2. * (kvec[idim] / knyquist);
288 w *= koverkny == 0.0 ? 1.0 : std::sin(koverkny) / (koverkny);
289 }
290 // res = pow(w,p);
291 double res = 1;
292 for (int i = 0; i < order; i++)
293 res *= w;
294 return res;
295 };
296
297 // Copy over
298 for (int idim = 0; idim < N; idim++) {
299 force_real[idim] = density_grid_fourier;
300 force_real[idim].add_memory_label("FFTWGrid::compute_force_from_density_fourier::force_real_" +
301 std::to_string(idim));
302 force_real[idim].set_grid_status_real(idim == 0);
303 }
304
305 // Loop over all local fourier grid cells
306#ifdef USE_OMP
307#pragma omp parallel for
308#endif
309 for (int islice = 0; islice < Local_nx; islice++) {
310 [[maybe_unused]] double kmag2;
311 [[maybe_unused]] std::array<double, N> kvec;
312 std::complex<FML::GRID::FloatType> I(0, 1);
313 for (auto && fourier_index : force_real[0].get_fourier_range(islice, islice + 1)) {
314 if (Local_x_start == 0 and fourier_index == 0)
315 continue; // DC mode (k=0)
316
317 force_real[0].get_fourier_wavevector_and_norm2_by_index(fourier_index, kvec, kmag2);
318 auto value = force_real[0].get_fourier_from_index(fourier_index);
319
320 // Divide by k^2 (different kernel choices here, fiducial is just 1/k^2)
321 if constexpr (kernel_choice == CONTINUOUS_GREENS_FUNCTION) {
322 value /= kmag2;
323 } else if constexpr (kernel_choice == CONTINUOUS_GREENS_FUNCTION_DECONVOLVE) {
324 double W = window_function(kvec);
325 value /= (kmag2 * W * W);
326 } else if constexpr (kernel_choice == DISCRETE_GREENS_FUNCTION_HOCKNEYEASTWOOD) {
327 double sum = 0.0;
328 for (int idim = 0; idim < N; idim++) {
329 double s = std::sin(kvec[idim] / (2.0 * double(Nmesh)));
330 sum += s * s;
331 }
332 sum *= 4.0 * double(Nmesh * Nmesh);
333 value /= sum;
334 } else if constexpr (kernel_choice == DISCRETE_GREENS_FUNCTION_HOCKNEYEASTWOOD_DECONVOLVE) {
335 double W = window_function(kvec);
336 double sum = 0.0;
337 for (int idim = 0; idim < N; idim++) {
338 double s = std::sin(kvec[idim] / (2.0 * double(Nmesh)));
339 sum += s * s;
340 }
341 sum *= 4.0 * double(Nmesh * Nmesh) * W * W;
342 value /= sum;
343 } else if constexpr (kernel_choice == DISCRETE_GREENS_FUNCTION_HAMMING) {
344 value *= 1.0 / kmag2;
345 } else if constexpr (kernel_choice == DISCRETE_GREENS_FUNCTION_HAMMING_DECONVOLVE) {
346 double W = window_function(kvec);
347 value *= 1.0 / (kmag2 * W * W);
348 } else {
349 FML::assert_mpi(
350 false,
351 "Unknown kernel_choice in compute_force_from_density_fourier. Method set at the "
352 "head of this function");
353 }
354
355 // Modify F[D] = kvec -> (8*sin(ki dx) - sin(2 ki dx))/6dx
356 if constexpr (kernel_choice == DISCRETE_GREENS_FUNCTION_HAMMING or
357 kernel_choice == DISCRETE_GREENS_FUNCTION_HAMMING_DECONVOLVE) {
358 for (int idim = 0; idim < N; idim++) {
359 kvec[idim] = (8.0 * std::sin(kvec[idim] / double(Nmesh)) - std::sin(2 * double(Nmesh))) /
360 6.0 * double(Nmesh);
361 }
362 }
363
364 // Compute force -ik/k^2 delta(k)
365 for (int idim = 0; idim < N; idim++) {
366 force_real[idim].set_fourier_from_index(
367 fourier_index, -I * value * FML::GRID::FloatType(kvec[idim] * norm_poisson_equation));
368 }
369 }
370 }
371
372 // Deal with DC mode
373 if (Local_x_start == 0)
374 for (int idim = 0; idim < N; idim++)
375 force_real[idim].set_fourier_from_index(0, 0.0);
376
377 // Fourier transform back to real space
378 for (int idim = 0; idim < N; idim++)
379 force_real[idim].fftw_c2r();
380 }
const int Nmesh
Definition: Main.cpp:48
int get_nmesh() const
Number of grid-cells per dimension.
Definition: FFTWGrid.h:1054
void add_memory_label(std::string label)
For memory logging: add a label to the grid.
Definition: FFTWGrid.h:318
void fftw_c2r(FFTWGrid< N > &in_grid, FFTWGrid< N > &out_grid)
Definition: FFTWGrid.h:1225
int interpolation_order_from_name(std::string density_assignment_method)
Get the interpolation order from a string holding the density_assignment_method (NGP,...

◆ compute_force_from_density_real()

template<int N>
void FML::NBODY::compute_force_from_density_real ( const FFTWGrid< N > &  density_grid_real,
std::array< FFTWGrid< N >, N > &  force_real,
std::string  density_assignment_method_used,
double  norm_poisson_equation 
)

Take a density grid in real space and returns the force \( \nabla \phi \) where \( \nabla^2 \phi = {\rm norm} \cdot \delta \) Different choices for what kernel to use for \( \nabla / \nabla^2\) are availiable, see the function body (is set too be a compile time option).

Fiducial choice is the continuous greens function \( 1/k^2\), but we can also choose to also devonvolve the window and discrete kernels (Hamming 1989; same as used in GADGET) and Hockney & Eastwood 1988. See e.g. 1603.00476 for a list.

Template Parameters
NThe dimension of the grid
Parameters
[in]density_grid_realThe density contrast in real space.
[out]force_realThe force in real space.
[in]density_assignment_method_usedThe density assignement we used to compute the density field. Needed only in case kernel_choice (defined in the body of this function) is not CONTINUOUS_GREENS_FUNCTION.
[in]norm_poisson_equationThe prefactor (norm) to the Poisson equation.

Definition at line 222 of file NBody.h.

225 {
226
227 FFTWGrid<N> density_grid_fourier = density_grid_real;
228 density_grid_fourier.add_memory_label("FFTWGrid::compute_force_from_density_real::density_grid_fourier");
229 density_grid_fourier.set_grid_status_real(true);
230 density_grid_fourier.fftw_r2c();
232 density_grid_fourier, force_real, density_assignment_method_used, norm_poisson_equation);
233 }
Class for holding grids and performing real-to-complex and complex-to-real FFTs using FFTW with MPI.
Definition: FFTWGrid.h:88
void set_grid_status_real(bool grid_is_a_real_grid)
Set the status of the grid: is it currently a real grid or a fourier grid?
Definition: FFTWGrid.h:468
void compute_force_from_density_fourier(const FFTWGrid< N > &density_grid_fourier, std::array< FFTWGrid< N >, N > &force_real, std::string density_assignment_method_used, double norm_poisson_equation)
Take a density grid in fourier space and returns the force where Different choices for what kernel ...
Definition: NBody.h:253

◆ compute_sigma_of_R()

double FML::NBODY::compute_sigma_of_R ( std::function< double(double)>  Pofk_of_kBox_over_volume,
double  R_mpch,
double  boxsize_mpch 
)

This computes the standard deviation of linear fluctuations smoothed over a sphere of radius, \( \sigma(R) \), by integrating the power-spectrum convolved with a top-hat windowfunction of radius \( R \).

We do this by solving \( \sigma^2(R) = \int \frac{k^3P(k)}{2\pi^2} |W(kR)|^2 d\log k \)

Parameters
[in]Pofk_of_kBox_over_volumeDimensionless power-spectrum \( P / V \) as function of the dimensionless scale \( kB \) where \(B\) is the boxsize and \(V = B^3\) is the box volume.
[in]R_mpchR in units of Mpc/h
[in]boxsize_mpchBoxsize in units of Mpc/h

Definition at line 1294 of file NBody.h.

1294 {
1298
1299 // We integrate from k = 1e-5 h/Mpc to k = 100 h/Mpc
1300 const double kBoxmin = 1e-5 * boxsize_mpch;
1301 const double kBoxmax = 1e2 * boxsize_mpch;
1302
1303 ODEFunction deriv = [&](double logkBox, [[maybe_unused]] const double * sigma2, double * dsigma2dlogk) {
1304 const double kBox = std::exp(logkBox);
1305 const double dimless_pofk = kBox * kBox * kBox / (2.0 * M_PI * M_PI) * Pofk_of_kBox_over_volume(kBox);
1306 const double kR = kBox * R_mpch / boxsize_mpch;
1307 const double window = kR > 0.0 ? 3.0 * (std::sin(kR) - kR * std::cos(kR)) / (kR * kR * kR) : 1.0;
1308 dsigma2dlogk[0] = dimless_pofk * window * window;
1309 return GSL_SUCCESS;
1310 };
1311
1312 // The initial conditions
1313 DVector sigmaini{0.0};
1314 DVector logkarr{std::log(kBoxmin), std::log(kBoxmax)};
1315
1316 // Solve the ODE
1317 ODESolver ode(1e-3, 1e-10, 1e-10);
1318 ode.solve(deriv, logkarr, sigmaini);
1319 auto sigma2 = ode.get_final_data_by_component(0);
1320
1321 return std::sqrt(sigma2);
1322 }
FML::INTERPOLATION::SPLINE::DVector DVector
Definition: Simulation.h:39
FML::SOLVERS::ODESOLVER::ODESolver ODESolver
Definition: Simulation.h:41
FML::SOLVERS::ODESOLVER::ODEFunction ODEFunction
Definition: Simulation.h:42
This is a wrapper around the GSL library to easily solve ODEs and return the data in whatever format ...
Definition: ODESolver.h:87
std::function< int(double, const double *, double *)> ODEFunction
Definition: ODESolver.h:33
std::vector< double > DVector
Definition: ODESolver.h:29
const double exp
Definition: Global.h:259

◆ DriftParticles() [1/3]

template<int N, class T >
void FML::NBODY::DriftParticles ( FML::PARTICLE::MPIParticles< T > &  part,
double  delta_time,
bool  periodic_box = true 
)

◆ DriftParticles() [2/3]

template<int N, class T >
void FML::NBODY::DriftParticles ( MPIParticles< T > &  part,
double  delta_time,
bool  periodic_box 
)

This moves the particles according to \( x_{\rm new} = x + v \Delta t \).

Note that we assume the velocities are in such units that \( v \Delta t\) is a dimensionless shift in [0,1).

Template Parameters
NThe dimension of the grid
TThe particle class
Parameters
[out]partMPIParticles containing the particles.
[in]delta_timeThe size of the timestep.
[in]periodic_boxIs the box periodic?

Definition at line 395 of file NBody.h.

395 {
396 if (part.get_npart() == 0)
397 return;
398 if (delta_time == 0.0)
399 return;
400
401 // Sanity check on particle
402 assert_mpi(FML::PARTICLE::GetNDIM(T()) == N,
403 "[DriftParticles] NDIM of particles and of grid does not match");
404
405 DriftParticles<N, T>(part.get_particles_ptr(), part.get_npart(), delta_time, periodic_box);
406
407 // Particles might have left the current task
409 }
T * get_particles_ptr()
Get a pointer to the first particle.
Definition: MPIParticles.h:451
void communicate_particles()
Communicate particles across CPU boundaries.
Definition: MPIParticles.h:607
size_t get_npart() const
Number of active particles on the local task.
Definition: MPIParticles.h:436

◆ DriftParticles() [3/3]

template<int N, class T >
void FML::NBODY::DriftParticles ( T *  p,
size_t  NumPart,
double  delta_time,
bool  periodic_box 
)

This moves the particles according to \( x_{\rm new} = x + v \Delta t \).

Note that we assume the velocities are in such units that \( v \Delta t\) is a dimensionless shift in [0,1). NB: after this methods is done the particles might have left the current task and must be communicated (this is done automatically if you use the MPIParticles version of this method).

Template Parameters
NThe dimension of the grid
TThe particle class
Parameters
[out]pPointer to the first particle.
[in]NumPartThe number of local particles.
[in]delta_timeThe size of the timestep.
[in]periodic_boxIs the box periodic?

Definition at line 427 of file NBody.h.

427 {
428 if (NumPart == 0)
429 return;
430 if (delta_time == 0.0)
431 return;
432
433 // Sanity check on particle
434 assert_mpi(FML::PARTICLE::GetNDIM(T()) == N,
435 "[DriftParticles] NDIM of particles and of grid does not match");
436 static_assert(FML::PARTICLE::has_get_pos<T>(),
437 "[DriftParticles] Particle class must have a get_pos method to use this method");
438 static_assert(FML::PARTICLE::has_get_vel<T>(),
439 "[DriftParticles] Particle class must have a get_vel method to use this method");
440
441 double max_disp = 0.0;
442#ifdef USE_OMP
443#pragma omp parallel for reduction(max : max_disp)
444#endif
445 for (size_t i = 0; i < NumPart; i++) {
446 auto * pos = FML::PARTICLE::GetPos(p[i]);
447 auto * vel = FML::PARTICLE::GetVel(p[i]);
448 for (int idim = 0; idim < N; idim++) {
449 double disp = vel[idim] * delta_time;
450 pos[idim] += disp;
451 max_disp = std::max(max_disp, std::abs(disp));
452
453 // Periodic wrap
454 if (periodic_box) {
455 if (pos[idim] >= 1.0)
456 pos[idim] -= 1.0;
457 if (pos[idim] < 0.0)
458 pos[idim] += 1.0;
459 }
460 }
461 }
462 FML::MaxOverTasks(&max_disp);
463
464 if (FML::ThisTask == 0)
465 std::cout << "[Drift] Max displacement: " << max_disp << "\n";
466 }
get_vel constexpr double * GetPos(...)
constexpr double * GetVel(...)
void MaxOverTasks(T *value)
Inplace max over tasks of a single value.
Definition: Global.h:121
int ThisTask
The MPI task number.
Definition: Global.cpp:33

◆ KickDriftKickNBodyStep()

template<int N, class T >
void FML::NBODY::KickDriftKickNBodyStep ( int  Nmesh,
MPIParticles< T > &  part,
double  delta_time,
std::string  density_assignment_method,
double  norm_poisson_equation 
)

Take a N-body step with a simple Kick-Drift-Kick method (this method serves mainly as an example for how one can do this).

  1. Particles to grid to get \( \delta \)
  2. Compute the Newtonian potential via \( \nabla^2 \Phi = {\rm norm} \cdot \delta \)
  3. Compute the force \( F = \nabla \Phi \)
  4. Move the particles using \( x \to x + v \Delta t \) and \( v \to v + F \Delta t \) This method assumes that the velocities are in units of (boxsize / time-step-unit), in other words that \( v\Delta t\) gives rise to a shift in [0,1). For cosmological N-body norm_poisson_equation depends on a and it should be set at the correct time. If one does simple sims with fixed time-step then the last kick of the previous step can be combined with the first kick of the current step to save one force evaluation per step (so basically two times as fast).
Template Parameters
NThe dimension of the grid.
TThe particle class.
Parameters
[in]NmeshThe gridsize to use for computing the density and force.
[out]partThe particles
[in]delta_timeThe time \( \Delta t \) we move forward.
[in]density_assignment_methodThe density assignement method (NGP, CIC, TSC, PCS or PQS).
[in]norm_poisson_equationA possible prefactor to the Poisson equation

Definition at line 88 of file NBody.h.

92 {
93
94 const bool periodic_box = true;
95
96 // Particles -> density field
97 auto nleftright =
99 FFTWGrid<N> density_grid_real(Nmesh, nleftright.first, nleftright.second);
100 density_grid_real.add_memory_label("FFTWGrid::KickDriftKickNBodyStep::density_grid_real");
101 density_grid_real.set_grid_status_real(true);
102 FML::INTERPOLATION::particles_to_grid<N, T>(part.get_particles_ptr(),
103 part.get_npart(),
104 part.get_npart_total(),
105 density_grid_real,
106 density_assignment_method);
107
108 // Density field -> force
109 std::array<FFTWGrid<N>, N> force_real;
111 density_grid_real, force_real, density_assignment_method, norm_poisson_equation);
112
113 // Update velocity of particles
114 KickParticles(force_real, part, delta_time * 0.5, density_assignment_method);
115
116 // Move particles (this does communication)
117 DriftParticles<N, T>(part, delta_time, periodic_box);
118
119 // Particles -> density field
120 FML::INTERPOLATION::particles_to_grid<N, T>(part.get_particles_ptr(),
121 part.get_npart(),
122 part.get_npart_total(),
123 density_grid_real,
124 density_assignment_method);
125
126 // Density field -> force
128 density_grid_real, force_real, density_assignment_method, norm_poisson_equation);
129
130 // Update velocity of particles
131 KickParticles(force_real, part, delta_time * 0.5, density_assignment_method);
132 }
size_t get_npart_total() const
Total number of active particles across all tasks.
Definition: MPIParticles.h:441
std::pair< int, int > get_extra_slices_needed_for_density_assignment(std::string density_assignment_method)
Compute how many extra slices we need in the FFTWGrid for a given density assignement / interpolation...
void compute_force_from_density_real(const FFTWGrid< N > &density_grid_real, std::array< FFTWGrid< N >, N > &force_real, std::string density_assignment_method_used, double norm_poisson_equation)
Take a density grid in real space and returns the force where Different choices for what kernel to ...
Definition: NBody.h:222
void KickParticles(std::array< FFTWGrid< N >, N > &force_grid, T *p, size_t NumPart, double delta_time, std::string interpolation_method)
This moves the particle velocities according to .
Definition: NBody.h:513

◆ KickParticles() [1/2]

template<int N, class T >
void FML::NBODY::KickParticles ( std::array< FFTWGrid< N >, N > &  force_grid,
MPIParticles< T > &  part,
double  delta_time,
std::string  interpolation_method 
)

This moves the particle velocities according to \( v_{\rm new} = v + F \Delta t \).

This method assumes the force is normalized such that \( F \Delta t \) has the same units as your v. If the flag free_force_grids is set in the source then we free up memory of the force grids after we have used them. The defalt is false.

Template Parameters
NThe dimension of the grid
TThe particle class
Parameters
[in]force_gridGrid containing the force.
[out]partMPIParticles containing the particles.
[in]delta_timeThe size of the timestep.
[in]interpolation_methodThe interpolation method for interpolating the force to the particle positions.

Definition at line 485 of file NBody.h.

488 {
489 if (delta_time == 0.0)
490 return;
491 KickParticles<N, T>(
492 force_grid, part.get_particles_ptr(), part.get_npart(), delta_time, interpolation_method);
493 }
const std::string interpolation_method
Definition: Main.cpp:64

◆ KickParticles() [2/2]

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

This moves the particle velocities according to \( v_{\rm new} = v + F \Delta t \).

This method assumes the force is normalized such that \( F \Delta t \) has the same units as your v. If the flag free_force_grids is set in the source then we free up memory of the force grids after we have used them. The defalt is false.

Template Parameters
NThe dimension of the grid
TThe particle class
Parameters
[in]force_gridThe force \( \nabla \Phi \).
[out]pPointer to the first particle.
[in]NumPartThe number of local particles.
[in]delta_timeThe size of the timestep.
[in]interpolation_methodThe interpolation method for interpolating the force to the particle positions.

Definition at line 513 of file NBody.h.

517 {
518
519 // Nothing to do if delta_time = 0.0
520 if (delta_time == 0.0)
521 return;
522
523 // Sanity check on particle
524 assert_mpi(FML::PARTICLE::GetNDIM(T()) == N, "[KickParticles] Dimension of particle and grid do not match");
525 static_assert(FML::PARTICLE::has_get_vel<T>(),
526 "[KickParticles] Particle must have velocity to use this method");
527
528 // Deallocate the force grids (after interpolating to the particles we don't need it here and probably
529 // not elsewhere so lets save some memory)
530 constexpr bool free_force_grids = false;
531
532 // Interpolate force to particle positions
533 for (int idim = 0; idim < N; idim++) {
534 force_grid[idim].communicate_boundaries();
535 }
536 std::array<std::vector<FML::GRID::FloatType>, N> force;
537 FML::INTERPOLATION::interpolate_grid_vector_to_particle_positions<N, T>(
538 force_grid, p, NumPart, force, interpolation_method);
539 if (free_force_grids) {
540 for (int idim = 0; idim < N; idim++) {
541 force_grid[idim].free();
542 }
543 }
544
545 double max_dvel = 0.0;
546#ifdef USE_OMP
547#pragma omp parallel for reduction(max : max_dvel)
548#endif
549 for (size_t i = 0; i < NumPart; i++) {
550 auto * vel = FML::PARTICLE::GetVel(p[i]);
551 for (int idim = 0; idim < N; idim++) {
552 double dvel = -force[idim][i] * delta_time;
553 max_dvel = std::max(max_dvel, std::abs(dvel));
554 vel[idim] += dvel;
555 }
556 }
557
558 FML::MaxOverTasks(&max_dvel);
559
560 if (FML::ThisTask == 0)
561 std::cout << "[Kick] Max delta_vel * delta_time : " << max_dvel * delta_time << "\n";
562 }

◆ NBodyInitialConditions() [1/2]

template<int N, class T >
void FML::NBODY::NBodyInitialConditions ( MPIParticles< T > &  part,
int  Npart_1D,
double  buffer_factor,
const FFTWGrid< N > &  delta_fourier,
std::vector< FFTWGrid< N > > &  phi_nLPT_potentials,
int  LPT_order,
double  box,
double  zini,
std::vector< double >  velocity_norms 
)

Generate particles from a given initial density field using Lagrangian perturbation theory.

We generate particles in [0,1) and velocities are given by \( v_{\rm code} = \frac{a^2 \frac{dx}{dt}}{H_0 B} \)

Template Parameters
NThe dimension we are working in.
TThe particle class. Must have methods get_pos, get_vel, get_D_1LPT and get_D_2LPT. But only get_pos data is required to exist. Return a nullptr if the data does not exist in the particle.
Parameters
[out]partParticle container for particles we are to create.
[in]Npart_1DNumber of particles per dimension (i.e. total is \( {\rm Npart}_{\rm 1D}^N \))
[in]buffer_factorHow many more particles to allocate?
[in]delta_fourierThe initial density field \( \delta(k,z_{\rm ini})\) in fourier space
[in]phi_nLPT_potentialsReturn the LPT potentials: 2LPT, 3LPTa, 3LPTb, ... If the vector has zero size then nothing will be returned.
[in]LPT_orderThe LPT order (1 or 2)
[in]boxThe boxsize (only for prining maximum displacement)
[in]ziniThe initial redshift
[in]velocity_normsA vector of the factors we need to multiply the nLPT displacement fields by to get velocities. E.g. \( 100 {\rm Box_in_Mpch} f_i(z_{\rm ini}) H(z_{\rm ini})/H_0 \cdot a_{\rm ini} \) to get peculiar velocities in km/s and \( f_i(z_{\rm ini}) H(z_{\rm ini})/H_0 \cdot a_{\rm ini}^2 \) to get the velocities we use as the fiducial choice in N-body.

Definition at line 694 of file NBody.h.

704 {
705
706 T tmp{};
707 if (FML::ThisTask == 0) {
708 std::cout << "\n";
709 std::cout << "#=====================================================\n";
710 std::cout << "#\n";
711 std::cout << "# .___ _________ \n";
712 std::cout << "# | | \\_ ___ \\ \n";
713 std::cout << "# | | / \\ \\/ \n";
714 std::cout << "# | | \\ \\____ \n";
715 std::cout << "# |___| \\______ / \n";
716 std::cout << "# \\/ \n";
717 std::cout << "#\n";
718 std::cout << "# Generating initial conditions for N-body\n";
719 std::cout << "# Order in LPT = " << LPT_order << "\n";
720 std::cout << "# The boxsize is " << box << " comoving Mpc/h\n";
721 std::cout << "# The initial redshift zini = " << zini << "\n";
722 if (FML::PARTICLE::has_get_pos<T>())
723 std::cout << "# Particle has [Position] x_code = x / Box ("
724 << sizeof(FML::PARTICLE::GetVel(tmp)[0]) * N << " bytes)\n";
725 if (FML::PARTICLE::has_get_vel<T>())
726 std::cout << "# Particle has [Velocity] v_code = a^2 dxdt / (H0 Box) ("
727 << sizeof(FML::PARTICLE::GetPos(tmp)[0]) * N << " bytes)\n";
728 if (FML::PARTICLE::has_set_mass<T>())
729 std::cout << "# Particle has [Mass] (" << sizeof(FML::PARTICLE::GetMass(tmp)) << " bytes)\n";
730 if (FML::PARTICLE::has_set_id<T>())
731 std::cout << "# Particle has [ID] (" << sizeof(FML::PARTICLE::GetID(tmp)) << " bytes)\n";
732 if (FML::PARTICLE::has_get_D_1LPT<T>())
733 std::cout << "# Particle has [1LPT Displacement field] ("
734 << sizeof(FML::PARTICLE::GetD_1LPT(tmp)[0]) * N << " bytes)\n";
735 if (FML::PARTICLE::has_get_D_2LPT<T>())
736 std::cout << "# Particle has [2LPT Displacement field] ("
737 << sizeof(FML::PARTICLE::GetD_2LPT(tmp)[0]) * N << " bytes)\n";
738 if (FML::PARTICLE::has_get_D_3LPTa<T>())
739 std::cout << "# Particle has [3LPTa Displacement field] ("
740 << sizeof(FML::PARTICLE::GetD_3LPTa(tmp)[0]) * N << " bytes)\n";
741 if (FML::PARTICLE::has_get_D_3LPTb<T>())
742 std::cout << "# Particle has [3LPTb Displacement field] ("
743 << sizeof(FML::PARTICLE::GetD_3LPTb(tmp)[0]) * N << " bytes)\n";
744 if (FML::PARTICLE::has_get_q<T>())
745 std::cout << "# Particle has [Lagrangian position] ("
746 << sizeof(FML::PARTICLE::GetLagrangianPos(tmp)[0]) * N << " bytes)\n";
747 std::cout << "# Total size of particle is " << FML::PARTICLE::GetSize(tmp) << " bytes\n";
748 std::cout << "# We will make " << Npart_1D << "^" << N << " particles\n";
749 std::cout << "# Plus a buffer with room for " << (buffer_factor - 1.0) * 100.0 << "%% more particles\n";
750 std::cout << "# We will allocate ~"
751 << buffer_factor * double(FML::PARTICLE::GetSize(tmp)) * double(FML::power(Npart_1D, N)) /
752 1e6 / double(FML::NTasks)
753 << " MB per task for the particles\n";
754 std::cout << "#\n";
755 std::cout << "#=====================================================\n";
756 std::cout << "\n";
757 }
758
759 // Sanity checks
760 const auto Nmesh = delta_fourier.get_nmesh();
761 assert_mpi(Nmesh > 0, "[NBodyInitialConditions] delta_fourier has to be already allocated");
762 assert_mpi(LPT_order == 1 or LPT_order == 2 or LPT_order == 3,
763 "[NBodyInitialConditions] Only 1LPT, 2LPT and 3LPT implemented so valid choices here are "
764 "LPT_order = 1, 2 or 3");
765 const std::string interpolation_method = "CIC"; // We use n-linear interpolation below
766 const auto nextra_cic =
768 assert_mpi(delta_fourier.get_n_extra_slices_left() >= nextra_cic.first and
769 delta_fourier.get_n_extra_slices_right() >= nextra_cic.second,
770 "[NBodyInitialConditions] We use CIC interpolation in this routine so the grid needs to have "
771 "atleast one extra slice on the right");
772
773 // Sanity check on particle
774 assert_mpi(FML::PARTICLE::GetNDIM(T()) == N,
775 "[NBodyInitialConditions] NDIM of particles and of grid does not match");
776 assert_mpi(FML::PARTICLE::has_get_pos<T>(),
777 "[NBodyInitialConditions] Particle class must have a get_pos method");
778
779 // The scalefactor and log(a) at the initial time
780 const double aini = 1.0 / (1.0 + zini);
781
782 FFTWGrid<N> phi_1LPT;
783 FFTWGrid<N> phi_2LPT;
784 FFTWGrid<N> phi_3LPTa;
785 FFTWGrid<N> phi_3LPTb;
786 if (LPT_order == 1) {
787 // Generate the 1LPT potential phi_1LPT = delta(k)/k^2
789 } else if (LPT_order == 2) {
790 // Generate the 1LPT potential phi_1LPT = delta(k)/k^2
792 // Generate the 2LPT potential phi_2LPT = -1/2k^2 F[phi_ii phi_jj - phi_ij^2]
794 } else if (LPT_order == 3) {
795 // Generate the 3LPT potentials phi_3LPTa, phi_3LPTb plus 3LPT curl term
796 // We ignore the curl term in this implementation for simplicity
797 const bool ignore_3LPT_curl_term = true;
798 std::array<FFTWGrid<N>, N> phi_3LPT_Avec_fourier;
799 FML::COSMOLOGY::LPT::compute_3LPT_potential_fourier<N>(delta_fourier,
800 phi_1LPT,
801 phi_2LPT,
802 phi_3LPTa,
803 phi_3LPTb,
804 phi_3LPT_Avec_fourier,
805 ignore_3LPT_curl_term);
806 }
807
808 //================================================================
809 // Function to compute the displacement from a LPT potential
810 // Frees the memory of phi_nLPT after its used
811 //================================================================
812 auto comp_displacement = [&]([[maybe_unused]] int nLPT,
813 FFTWGrid<N> & phi_nLPT,
814 std::array<std::vector<FML::GRID::FloatType>, N> & displacements_nLPT) {
815 // Generate Psi from phi
816 std::array<FFTWGrid<N>, N> Psi_nLPT_vector;
817 FML::COSMOLOGY::LPT::from_LPT_potential_to_displacement_vector<N>(phi_nLPT, Psi_nLPT_vector);
818 for (int idim = 0; idim < N; idim++) {
819 Psi_nLPT_vector[idim].communicate_boundaries();
820 }
821 phi_nLPT.free();
822
823 // Interpolate it to particle Lagrangian positions
824 FML::INTERPOLATION::interpolate_grid_vector_to_particle_positions<N, T>(Psi_nLPT_vector,
825 part.get_particles_ptr(),
826 part.get_npart(),
827 displacements_nLPT,
829 };
830
831 auto add_displacement = [&]([[maybe_unused]] int nLPT,
832 char type,
833 std::array<std::vector<FML::GRID::FloatType>, N> & displacements_nLPT,
834 double vfac_nLPT) -> void {
835 // Generate Psi from phi
836 // Add displacement to particle position
837 double max_disp_nLPT = 0.0;
838 double max_vel_nLPT = 0.0;
839 auto * part_ptr = part.get_particles_ptr();
840#ifdef USE_OMP
841#pragma omp parallel for reduction(max : max_disp_nLPT, max_vel_nLPT)
842#endif
843 for (size_t ind = 0; ind < part.get_npart(); ind++) {
844
845 // Fetch displacement
846 std::array<double, N> disp;
847 for (int idim = 0; idim < N; idim++) {
848 disp[idim] = displacements_nLPT[idim][ind];
849 }
850
851 // Add to position (particle must have position)
852 auto * pos = FML::PARTICLE::GetPos(part_ptr[ind]);
853 if (ind == 0 and FML::ThisTask == 0)
854 std::cout << "Adding " << std::to_string(nLPT) << "LPT position to particle\n";
855 for (int idim = 0; idim < N; idim++) {
856 const double dpos_nLPT = disp[idim];
857 pos[idim] += dpos_nLPT;
858
859 // Periodic BC
860 if (pos[idim] >= 1.0)
861 pos[idim] -= 1.0;
862 if (pos[idim] < 0.0)
863 pos[idim] += 1.0;
864
865 // Compute maximum displacement
866 if (std::fabs(dpos_nLPT) > max_disp_nLPT)
867 max_disp_nLPT = std::fabs(dpos_nLPT);
868 }
869
870 // Add to velocity (if it exists)
871 if constexpr (FML::PARTICLE::has_get_vel<T>()) {
872 if (ind == 0 and FML::ThisTask == 0)
873 std::cout << "Adding " << std::to_string(nLPT) << "LPT velocity to particle\n";
874 auto * vel = FML::PARTICLE::GetVel(part_ptr[ind]);
875 for (int idim = 0; idim < N; idim++) {
876 vel[idim] += vfac_nLPT * disp[idim];
877
878 if (std::fabs(vfac_nLPT * disp[idim]) > max_vel_nLPT)
879 max_vel_nLPT = std::fabs(vfac_nLPT * disp[idim]);
880 }
881 }
882
883 // Store displacement fields at particle (if it exists)
884 // This is needed if we want to do COLA
885 if (nLPT == 1) {
886 if constexpr (FML::PARTICLE::has_get_D_1LPT<T>()) {
887 if (ind == 0 and FML::ThisTask == 0)
888 std::cout << "Storing 1LPT displacment field in particle\n";
889 auto * D = FML::PARTICLE::GetD_1LPT(part_ptr[ind]);
890 for (int idim = 0; idim < N; idim++) {
891 D[idim] = disp[idim];
892 }
893 }
894 }
895
896 if (nLPT == 2) {
897 if constexpr (FML::PARTICLE::has_get_D_2LPT<T>()) {
898 if (ind == 0 and FML::ThisTask == 0)
899 std::cout << "Storing 2LPT displacment field in particle\n";
900 auto * D2 = FML::PARTICLE::GetD_2LPT(part_ptr[ind]);
901 for (int idim = 0; idim < N; idim++) {
902 D2[idim] = disp[idim];
903 }
904 }
905 }
906
907 if (nLPT == 3 and type == 'a') {
908 if constexpr (FML::PARTICLE::has_get_D_3LPTa<T>()) {
909 if (ind == 0 and FML::ThisTask == 0)
910 std::cout << "Storing 3LPTa displacment field in particle\n";
911 auto * D3a = FML::PARTICLE::GetD_3LPTa(part_ptr[ind]);
912 for (int idim = 0; idim < N; idim++) {
913 D3a[idim] = disp[idim];
914 }
915 }
916 }
917
918 if (nLPT == 3 and type == 'b') {
919 if constexpr (FML::PARTICLE::has_get_D_3LPTb<T>()) {
920 if (ind == 0 and FML::ThisTask == 0)
921 std::cout << "Storing 3LPTb displacment field in particle\n";
922 auto * D3b = FML::PARTICLE::GetD_3LPTb(part_ptr[ind]);
923 for (int idim = 0; idim < N; idim++) {
924 D3b[idim] = disp[idim];
925 }
926 }
927 }
928 }
929
930 // Output the maximum displacment and velocity
931 FML::MaxOverTasks(&max_disp_nLPT);
932 FML::MaxOverTasks(&max_vel_nLPT);
933 if (FML::ThisTask == 0)
934 std::cout << "Maximum " << std::to_string(nLPT) << "LPT displacements: " << max_disp_nLPT * box
935 << " Mpc/h\n";
936 if (FML::ThisTask == 0)
937 std::cout << "Maximum " << std::to_string(nLPT)
938 << "LPT velocity: " << max_vel_nLPT * 100.0 * box / aini << " km/s peculiar\n";
939 };
940
941 // Create particles
943 part.info();
944
945 // Set unique IDs if we have that availiable in the particles
946 if constexpr (FML::PARTICLE::has_set_id<T>()) {
947 if (FML::ThisTask == 0)
948 std::cout << "Storing unique ID in particle\n";
949 long long int npart_local = part.get_npart();
950 auto part_per_task = FML::GatherFromTasks(&npart_local);
951 long long int id_start = 0;
952 for (int i = 0; i < FML::ThisTask; i++)
953 id_start += part_per_task[i];
954 long long int count = 0;
955 for (auto & p : part) {
956 FML::PARTICLE::SetID(p, id_start + count++);
957 }
958 }
959
960 // Set mass if we have that availiable in the particles
961 if constexpr (FML::PARTICLE::has_set_mass<T>()) {
962 if (FML::ThisTask == 0)
963 std::cout << "Storing Mass (1.0) in particle\n";
964 for (auto & p : part) {
966 }
967 }
968
969 // Set Lagrangian position of the particle if we have that availiable
970 if constexpr (FML::PARTICLE::has_get_q<T>()) {
971 if (FML::ThisTask == 0)
972 std::cout << "Storing Lagrangian position q in particle\n";
973 for (auto & p : part) {
974 auto pos = FML::PARTICLE::GetPos(p);
976 for (int idim = 0; idim < N; idim++)
977 q[idim] = pos[idim];
978 }
979 }
980
981 // Compute and add displacements
982 // NB: we must do this in one go as add_displacement changes the position of the particles
983 std::array<std::vector<FML::GRID::FloatType>, N> displacements_1LPT;
984 if (LPT_order >= 1) {
985 const int nLPT = 1;
986 comp_displacement(nLPT, phi_1LPT, displacements_1LPT);
987 }
988
989 std::array<std::vector<FML::GRID::FloatType>, N> displacements_2LPT;
990 if (LPT_order >= 2) {
991 const int nLPT = 2;
992 // Store potential if asked for
993 if (phi_nLPT_potentials.size() > 0)
994 phi_nLPT_potentials[0] = phi_2LPT;
995 comp_displacement(nLPT, phi_2LPT, displacements_2LPT);
996 }
997
998 std::array<std::vector<FML::GRID::FloatType>, N> displacements_3LPTa;
999 if (LPT_order >= 3) {
1000 const int nLPT = 3;
1001 // Store potential if asked for
1002 if (phi_nLPT_potentials.size() > 1)
1003 phi_nLPT_potentials[1] = phi_3LPTa;
1004 comp_displacement(nLPT, phi_3LPTa, displacements_3LPTa);
1005 }
1006
1007 std::array<std::vector<FML::GRID::FloatType>, N> displacements_3LPTb;
1008 if (LPT_order >= 3) {
1009 const int nLPT = 3;
1010 // Store potential if asked for
1011 if (phi_nLPT_potentials.size() > 2)
1012 phi_nLPT_potentials[2] = phi_3LPTb;
1013 comp_displacement(nLPT, phi_3LPTb, displacements_3LPTb);
1014 }
1015
1016 if (LPT_order >= 1) {
1017 const int nLPT = 1;
1018 add_displacement(nLPT, 0, displacements_1LPT, velocity_norms[0]);
1019 for (auto & d : displacements_1LPT) {
1020 d.clear();
1021 d.shrink_to_fit();
1022 }
1023 }
1024
1025 if (LPT_order >= 2) {
1026 const int nLPT = 2;
1027 add_displacement(nLPT, 0, displacements_2LPT, velocity_norms[1]);
1028 for (auto & d : displacements_2LPT) {
1029 d.clear();
1030 d.shrink_to_fit();
1031 }
1032 }
1033
1034 if (LPT_order >= 3) {
1035 const int nLPT = 3;
1036 add_displacement(nLPT, 'a', displacements_3LPTa, velocity_norms[2]);
1037 for (auto & d : displacements_3LPTa) {
1038 d.clear();
1039 d.shrink_to_fit();
1040 }
1041 }
1042
1043 if (LPT_order >= 3) {
1044 const int nLPT = 3;
1045 add_displacement(nLPT, 'b', displacements_3LPTb, velocity_norms[3]);
1046 for (auto & d : displacements_3LPTb) {
1047 d.clear();
1048 d.shrink_to_fit();
1049 }
1050 }
1051
1052 // Communicate particles (they might have left the current task)
1053 part.communicate_particles();
1054 }
const double box
Definition: Main.cpp:45
const int Npart_1D
Definition: Main.cpp:50
const double buffer_factor
Definition: Main.cpp:51
const double zini
Definition: Main.cpp:40
const double aini
Definition: Main.cpp:41
int get_n_extra_slices_right() const
How many extra slices we have allocated to the right.
Definition: FFTWGrid.h:1265
int get_n_extra_slices_left() const
How many extra slices we have allocated to the left.
Definition: FFTWGrid.h:1260
void create_particle_grid(int Npart_1D, double buffer_factor, double xmin_local, double xmax_local)
Create particles in a rectangular grid spread across all tasks buffer_factor is how much extra to al...
Definition: MPIParticles.h:478
void info()
Show some info about the class.
Definition: MPIParticles.h:213
void compute_1LPT_potential_fourier(const FFTWGrid< N > &delta_fourier, FFTWGrid< N > &phi_1LPT_fourier)
Generate the 1LPT potential defined as and .
void compute_2LPT_potential_fourier(const FFTWGrid< N > &delta_fourier, FFTWGrid< N > &phi_2LPT_fourier)
Generate the 2LPT potential defined as and .
constexpr double * GetLagrangianPos(...)
constexpr double * GetD_3LPTa(...)
constexpr void SetMass(...)
set_id constexpr int GetID(...)
constexpr double * GetD_3LPTb(...)
constexpr double * GetD_2LPT(...)
set_mass constexpr double GetMass(Args &... args)
constexpr void SetID(...)
double xmax_domain
The local extent of the domain (global domain goes from 0 to 1)
Definition: Global.cpp:42
constexpr long long int power(int base, int exponent)
Simple integer a^b power-function by squaring.
Definition: Global.h:187
std::vector< T > GatherFromTasks(T *value)
Gather a single value from all tasks. Value from task ThisTask is stored in values[ThisTask].
Definition: Global.h:109
double xmin_domain
The local extent of the domain (global domain goes from 0 to 1)
Definition: Global.cpp:41
int NTasks
Total number of MPI tasks.
Definition: Global.cpp:34

◆ NBodyInitialConditions() [2/2]

template<int N, class T >
void FML::NBODY::NBodyInitialConditions ( 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::vector< double >  velocity_norms 
)

Generate particles from a given power-spectrum using Lagrangian perturbation theory.

We generate particles in [0,1) and velocities are given by \( v_{\rm code} = \frac{a^2 \frac{dx}{dt}}{H_0 B} \)

Template Parameters
NThe dimension we are working in.
TThe particle class. Must have methods get_pos, get_vel, get_D_1LPT and get_D_2LPT. But only get_pos data is required to exist. Return a nullptr if the data does not exist in the particle.
Parameters
[out]partParticle container for particles we are to create.
[in]Npart_1DNumber of particles per dimension (i.e. total is \( {\rm Npart}_{\rm 1D}^N \))
[in]buffer_factorHow many more particles to allocate?
[in]NmeshThe grid to generate the IC on
[in]fix_amplitudeAmplitude fixed? Only random phases if true.
[in]rngRandom number generator
[in]Pofk_of_kBox_over_Pofk_primordalThe ratio of the power-spectrum (for delta) at the time you want the density field to be created at to the primordial one (the function above).
[in]Pofk_of_kBox_over_volume_primordialThe dimensionless function \( P/V\) where \( V = B^N\) is the box volume as function of the dimensionless wavenumber \( kB \) where \( B \) is the boxsize and \( P(k) \) is the primordial power-spectrum for \(\Phi\).
[in]LPT_orderThe LPT order (1 or 2)
[in]type_of_random_fieldWhat random field: gaussian, local, equilateral, orthogonal
[in]fNLIf non-gaussianity the value of fNL
[in]boxThe boxsize (only for prining maximum displacement)
[in]ziniThe initial redshift
[in]velocity_normsA vector of the factors we need to multiply the nLPT displacement fields by to get velocities. E.g. \( 100 {\rm Box_in_Mpch} f_i(z_{\rm ini}) H(z_{\rm ini})/H_0 \cdot a_{\rm ini} \) to get peculiar velocities in km/s and \( f_i(z_{\rm ini}) H(z_{\rm ini})/H_0 \cdot a_{\rm ini}^2 \) to get the velocities we use as the fiducial choice in N-body. The order is: 1LPT, 2LPT, 3LPTa, 3LPTb

Definition at line 609 of file NBody.h.

624 {
625
626 // Some sanity checks
627 assert_mpi(Npart_1D > 0 and Nmesh > 0 and zini >= 0.0 and rng != nullptr and box > 0.0,
628 "[NBodyInitialConditions] Invalid parameters");
629
630 // Generate the random field first
632 FFTWGrid<N> delta_fourier(Nmesh, nextra.first, nextra.second);
633 delta_fourier.add_memory_label("FFTWGrid::NBodyInitialConditions::delta_fourier");
634 delta_fourier.set_grid_status_real(false);
635
636 // Make a gaussian or non-local non-gaussian random field in fourier space
637 if (type_of_random_field == "gaussian") {
638 auto Pofk_of_kBox_over_volume = [&](double kBox) {
639 if (kBox == 0.0)
640 return 0.0;
641 return Pofk_of_kBox_over_Pofk_primordal(kBox) * Pofk_of_kBox_over_volume_primordial(kBox);
642 };
644 delta_fourier, rng, Pofk_of_kBox_over_volume, fix_amplitude);
645 } else {
647 delta_fourier,
648 rng,
649 Pofk_of_kBox_over_Pofk_primordal,
650 Pofk_of_kBox_over_volume_primordial,
652 fNL,
654 }
655
656 // Generate IC from a given fourier grid
657 std::vector<FFTWGrid<N>> phi_nLPT_potentials;
658 NBodyInitialConditions<N, T>(part,
659 Npart_1D,
661 delta_fourier,
662 phi_nLPT_potentials,
663 LPT_order,
664 box,
665 zini,
666 velocity_norms);
667 }
const std::string type_of_random_field
Definition: Main.cpp:45
const bool fix_amplitude
Definition: Main.cpp:44
const double fNL
Definition: Main.cpp:46
void generate_gaussian_random_field_fourier(FFTWGrid< N > &grid, RandomGenerator *rng, std::function< double(double)> Pofk_of_kBox_over_volume, bool fix_amplitude)
void generate_nonlocal_gaussian_random_field_fourier_cosmology(FFTWGrid< N > &delta_fourier, RandomGenerator *rng, std::function< double(double)> Pofk_of_kBox_over_Pofk_primordal, std::function< double(double)> Pofk_of_kBox_over_volume_primordial, bool fix_amplitude, double fNL, std::string type_of_fnl, double u=0.0, std::vector< double > kernel_values={})
This computes a non-gaussian density field at any redshift.

◆ YoshidaNBodyStep()

template<int N, class T >
void FML::NBODY::YoshidaNBodyStep ( int  Nmesh,
MPIParticles< T > &  part,
double  delta_time,
std::string  density_assignment_method,
double  norm_poisson_equation 
)

Take a N-body step with a 4th order symplectic Yoshida method.

This method is mainly an illustration, for using this with cosmology we should take into account that norm_poisson_equation is a function of time

Template Parameters
NThe dimension of the grid.
TThe particle class.
Parameters
[in]NmeshThe gridsize to use for computing the density and force.
[out]partThe particles
[in]delta_timeThe time \( \Delta t \) we move forward.
[in]density_assignment_methodThe density assignement method (NGP, CIC, TSC, PCS or PQS).
[in]norm_poisson_equationA possible prefactor to the Poisson equation

Definition at line 149 of file NBody.h.

153 {
154
155 const bool periodic_box = true;
156
157 // The Yoshida coefficients
158 const double w1 = 1.0 / (2 - std::pow(2.0, 1.0 / 3.0));
159 const double w0 = 1.0 - 2.0 * w1;
160 const double c1 = w1 / 2.0, c4 = c1;
161 const double c2 = (w0 + w1) / 2.0, c3 = c2;
162 const double d1 = w1, d3 = d1;
163 const double d2 = w0;
164
165 // They must sum to unity
166 assert(std::fabs(c1 + c2 + c3 + c4 - 1.0) < 1e-10);
167 assert(std::fabs(d1 + d2 + d3 - 1.0) < 1e-10);
168
169 // Set up a density grid to use
170 auto nleftright =
172 FFTWGrid<N> density_grid_real(Nmesh, nleftright.first, nleftright.second);
173 density_grid_real.add_memory_label("FFTWGrid::YoshidaNBodyStep::density_grid_real");
174 density_grid_real.set_grid_status_real(true);
175
176 // Perform one step: delta_time_pos is the advance for pos positions and delta_time_vel is for velocity
177 auto one_step = [&](double delta_time_pos, double delta_time_vel, double norm_poisson) {
178 // Move particles (this does communication)
179 DriftParticles<N, T>(part, delta_time_pos, periodic_box);
180
181 // Particles -> density field
182 FML::INTERPOLATION::particles_to_grid<N, T>(part.get_particles_ptr(),
183 part.get_npart(),
184 part.get_npart_total(),
185 density_grid_real,
186 density_assignment_method);
187 // Density field -> force
188 std::array<FFTWGrid<N>, N> force_real;
189 compute_force_from_density_real(density_grid_real, force_real, density_assignment_method, norm_poisson);
190
191 // Update velocity of particles
192 KickParticles(force_real, part, delta_time_vel, density_assignment_method);
193 };
194
195 // The norm_poisson_equation in a cosmo sim depends on [aexp] so this should be changed
196 one_step(delta_time * c1, delta_time * d1, norm_poisson_equation);
197 one_step(delta_time * c2, delta_time * d2, norm_poisson_equation);
198 one_step(delta_time * c3, delta_time * d3, norm_poisson_equation);
199
200 // Move particles (this does communication)
201 DriftParticles<N, T>(part, delta_time * c4, periodic_box);
202 }
float w0
Definition: example.py:10