FML
Data Structures | Typedefs | Functions
FML::GRID Namespace Reference

This namespace contains various grids for holding data shared across different tasks. More...

Data Structures

class  FFTWGrid
 Class for holding grids and performing real-to-complex and complex-to-real FFTs using FFTW with MPI. More...
 
class  FourierRange
 For range based for-loops over the fourier-grid. Loop over all local cells. More...
 
class  LoopIteratorFourier
 An iterator that deal with looping through the fourier grid. More...
 
class  LoopIteratorReal
 An iterator that deal with looping through a real grid. More...
 
class  MPIGrid
 A simple multidimensional grid-class for any type that works over MPI. More...
 
class  MPIMultiGrid
 A stack of _Nlevel MPIGrids with \( N^{\rm NDIM} / 2^{\rm Level} \) cells in each level. More...
 
class  RealRange
 For range based for-loops over the real-grid. Loop over all local cells. More...
 

Typedefs

using IndexInt = long long int
 

Functions

template<int N>
void fftw_r2c (FFTWGrid< N > &in_grid, FFTWGrid< N > &out_grid)
 
template<int N>
void fftw_c2r (FFTWGrid< N > &in_grid, FFTWGrid< N > &out_grid)
 
void init_fftw (int *argc, char ***argv)
 
void finalize_fftw ()
 
void set_fftw_nthreads (int nthreads)
 
template<class T >
double AbsoluteValue (T &x)
 
 OPS (+=)
 
 OPS (-=)
 
OPS * OPS (/=)#undef OPS#define OPS(OP) \ template< int NDIM, class T > \ template< class U > \ auto MPIGrid< NDIM, T >::operator OP(const U &rhs) ->MPIGrid< NDIM, T > &{ \ for(IndexInt i=0;i< _NtotLocal;i++) \ this->_y[_NtotLocalLeft+i] OP rhs;\ return *this;\ } OPS(+=
 
template<int N>
void whitening_fourier_space (FFTWGrid< N > &fourier_grid)
 Take a fourier grid and divide each mode by its norm: \( f(k) \to f(k) / |f(k)| \). More...
 
template<int N>
void smoothing_filter_fourier_space (FFTWGrid< N > &fourier_grid, double smoothing_scale, std::string smoothing_method)
 Low-pass filters (tophat, gaussian, sharpk) More...
 
template<int N>
void convolution_fourier_space (const FFTWGrid< N > &fourier_grid_f, const FFTWGrid< N > &fourier_grid_g, FFTWGrid< N > &fourier_grid_result)
 From two fourier grids, f and g, compute the convolution \( f(k) * g(k) = \int d^{\rm N}q f(q) g(k-q) \) This is done via multuplication in reals-space. More...
 
template<int N>
void convolution_real_space (const FFTWGrid< N > &real_grid_f, const FFTWGrid< N > &real_grid_g, FFTWGrid< N > &real_grid_result)
 From two real grids, f and g, compute the convolution \( f(x) * g(x) = \int d^{\rm N}y f(y) g(x-y) \) This is done via multiplication in fourier-space. More...
 
template<int N>
void compute_grid_PDF (const FFTWGrid< N > &real_grid, int nbins, std::vector< double > &x, std::vector< double > &pdf)
 This computes the PDF of whatever quantity is in the grid (e.g. More...
 

Detailed Description

This namespace contains various grids for holding data shared across different tasks.

Typedef Documentation

◆ IndexInt

using FML::GRID::IndexInt = typedef long long int

Definition at line 24 of file MPIGrid.h.

Function Documentation

◆ AbsoluteValue()

template<class T >
double FML::GRID::AbsoluteValue ( T &  x)

Definition at line 28 of file MPIGrid.h.

28 {
29 return std::abs(x);
30 }
test.x
x
Definition: test.py:18

◆ compute_grid_PDF()

template<int N>
void FML::GRID::compute_grid_PDF ( const FFTWGrid< N > &  real_grid,
int  nbins,
std::vector< double > &  x,
std::vector< double > &  pdf 
)

This computes the PDF of whatever quantity is in the grid (e.g.

the density if its a density grid) The binning is set to be linear. The range is set by the values we find in the grid This realy belongs in the namespace CORRELATIONFUNCTIONS so should move it there

Template Parameters
NThe dimension of the grid
Parameters
[in]real_grid
[in]nbinsNumber of bins
[out]xBins for the quantity we are computing the PDF of
[out]pdfThe binned PDF normalized such that \( \int_{-\infty}^\infty p(x)dx = 1\).

Definition at line 245 of file SmoothingFourier.h.

245 {
246
247 // Multiply the two grids in real space
248 auto Local_nx = real_grid.get_local_nx();
249
250 // Find minimum and maximum value in the grid
251 double grid_min = std::numeric_limits<double>::max();
252 double grid_max = -grid_min;
253#ifdef USE_OMP
254#pragma omp parallel for reduction(max : grid_max) reduction(min : grid_min)
255#endif
256 for (int islice = 0; islice < Local_nx; islice++) {
257 for (auto && real_index : real_grid.get_real_range(islice, islice + 1)) {
258 auto value = real_grid.get_real_from_index(real_index);
259 grid_min = std::min(grid_min, value);
260 grid_max = std::max(grid_max, value);
261 }
262 }
263 FML::MinOverTasks(&grid_min);
264 FML::MaxOverTasks(&grid_max);
265
266 // Set up binning
267 x.resize(nbins);
268 pdf.resize(nbins, 0.0);
269 for (int i = 0; i < nbins; i++) {
270 x[i] = grid_min + (grid_max - grid_min) / double(nbins) * (i + 0.5);
271 }
272
273 // For binning over threads
274 std::vector<std::vector<double>> pdfthreads(FML::NThreads, std::vector<double>(nbins, 0.0));
275
276#ifdef USE_OMP
277#pragma omp parallel for
278#endif
279 for (int islice = 0; islice < Local_nx; islice++) {
280 int id = 0;
281#ifdef USE_OMP
282 id = omp_get_thread_num();
283#endif
284 for (auto && real_index : real_grid.get_real_range(islice, islice + 1)) {
285 auto value = real_grid.get_real_from_index(real_index);
286 int ibin = int((value - grid_min) / (grid_max - grid_min) * nbins);
287 if (ibin >= 0 and ibin < nbins)
288 pdfthreads[id][ibin] += 1;
289 }
290 }
291
292 // Sum up over threads
293 for (int i = 0; i < FML::NThreads; i++) {
294 for (int j = 0; j < nbins; j++) {
295 pdf[j] += pdfthreads[i][j];
296 }
297 }
298
299 // Sum over tasks
300#ifdef USE_MPI
301 MPI_Allreduce(MPI_IN_PLACE, pdf.data(), nbins, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
302#endif
303
304 // Normalize so that the PDF integrates to unity
305 const double dx = (grid_max - grid_min) / double(nbins);
306 double integral = 0.0;
307 for (int i = 0; i < nbins; i++) {
308 integral += pdf[i] * dx;
309 }
310 for (int i = 0; i < nbins; i++) {
311 pdf[i] /= integral;
312 }
313 }
FML::GRID::FFTWGrid::get_real_from_index
FloatType get_real_from_index(const IndexIntType index) const
Fetch value in grid by (local) index.
Definition: FFTWGrid.h:1253
FML::GRID::FFTWGrid::get_local_nx
ptrdiff_t get_local_nx() const
Number of local x-slices in the real grid.
Definition: FFTWGrid.h:1079
FML::GRID::FFTWGrid::get_real_range
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
FML::MinOverTasks
void MinOverTasks(T *value)
Inplace min over tasks of a single value.
Definition: Global.h:134
FML::NThreads
int NThreads
Total number of threads availiable.
Definition: Global.cpp:35
FML::MaxOverTasks
void MaxOverTasks(T *value)
Inplace max over tasks of a single value.
Definition: Global.h:121

◆ convolution_fourier_space()

template<int N>
void FML::GRID::convolution_fourier_space ( const FFTWGrid< N > &  fourier_grid_f,
const FFTWGrid< N > &  fourier_grid_g,
FFTWGrid< N > &  fourier_grid_result 
)

From two fourier grids, f and g, compute the convolution \( f(k) * g(k) = \int d^{\rm N}q f(q) g(k-q) \) This is done via multuplication in reals-space.

We allocate one temporary grid and perform 3 fourier tranforms.

Template Parameters
Ndimension of the grid
Parameters
[in]fourier_grid_fThe fourier grid f
[in]fourier_grid_gThe fourier grid g
[out]fourier_grid_resultFourier grid containing the convolution of the two gridsf

Definition at line 134 of file SmoothingFourier.h.

136 {
137
138 bool f_and_g_are_the_same_grid = (&fourier_grid_f == &fourier_grid_g);
139
140 // Make copies: tmp contains f and result contains g
141 // unless f = g for which we can don't need the copy
142 // and save doing 1 FFT
143 FFTWGrid<N> tmp;
144 fourier_grid_result = fourier_grid_g;
145 fourier_grid_result.add_memory_label("FFTWGrid::convolution_fourier_space::fourier_grid_result");
146
147 // Fourier transform to real space
148 fourier_grid_result.fftw_c2r();
149 if (not f_and_g_are_the_same_grid) {
150 tmp = fourier_grid_f;
151 tmp.add_memory_label("FFTWGrid::convolution_fourier_space::tmp");
152 tmp.fftw_c2r();
153 }
154
155 // Multiply the two grids in real space
156 auto Local_nx = fourier_grid_result.get_local_nx();
157#ifdef USE_OMP
158#pragma omp parallel for
159#endif
160 for (int islice = 0; islice < Local_nx; islice++) {
161 for (auto && real_index : fourier_grid_result.get_real_range(islice, islice + 1)) {
162 auto g_real = fourier_grid_result.get_real_from_index(real_index);
163 auto f_real = g_real;
164 if (not f_and_g_are_the_same_grid)
165 f_real = tmp.get_real_from_index(real_index);
166 fourier_grid_result.set_real_from_index(real_index, f_real * g_real);
167 }
168 }
169
170 // Transform back to obtain the desired convolution
171 fourier_grid_result.fftw_r2c();
172 }
FML::GRID::FFTWGrid
Class for holding grids and performing real-to-complex and complex-to-real FFTs using FFTW with MPI.
Definition: FFTWGrid.h:88
FML::GRID::FFTWGrid::fftw_c2r
void fftw_c2r()
Perform complex-to-real fourier transform.
Definition: FFTWGrid.h:946
FML::GRID::FFTWGrid::set_real_from_index
void set_real_from_index(const IndexIntType ind, const FloatType value)
Set value in grid from (local) index.
Definition: FFTWGrid.h:1020
FML::GRID::FFTWGrid::add_memory_label
void add_memory_label(std::string label)
For memory logging: add a label to the grid.
Definition: FFTWGrid.h:318
FML::GRID::FFTWGrid::fftw_r2c
void fftw_r2c()
Perform real-to-complex fourier transform.
Definition: FFTWGrid.h:881

◆ convolution_real_space()

template<int N>
void FML::GRID::convolution_real_space ( const FFTWGrid< N > &  real_grid_f,
const FFTWGrid< N > &  real_grid_g,
FFTWGrid< N > &  real_grid_result 
)

From two real grids, f and g, compute the convolution \( f(x) * g(x) = \int d^{\rm N}y f(y) g(x-y) \) This is done via multiplication in fourier-space.

We allocate one temporary grid and perform 3 fourier tranforms. We can merge this with convolution_fourier_space and just have one method and get do the real or fourier space convolution depening on the status of the grids, but its easy to forget to set the status so we have two methods for this.

Template Parameters
Ndimension of the grid
Parameters
[in]real_grid_fThe real grid f
[in]real_grid_gThe real grid g
[out]real_grid_resultReal grid containing the convolution of the two grids.

Definition at line 190 of file SmoothingFourier.h.

192 {
193
194 bool f_and_g_are_the_same_grid = (&real_grid_f == &real_grid_g);
195
196 // Make copies: tmp contains f and result contains g
197 // unless f = g for which we can don't need the copy
198 // and save doing 1 FFT
199 FFTWGrid<N> tmp;
200 real_grid_result = real_grid_g;
201 real_grid_result.add_memory_label("FFTWGrid::convolution_real_space::real_grid_result");
202
203 // Fourier transform to fourier space
204 real_grid_result.fftw_r2c();
205 if (not f_and_g_are_the_same_grid) {
206 tmp = real_grid_f;
207 tmp.add_memory_label("FFTWGrid::convolution_real_space::tmp");
208 tmp.fftw_r2c();
209 }
210
211 // Multiply the two grids in fourier space
212 auto Local_nx = real_grid_result.get_local_nx();
213#ifdef USE_OMP
214#pragma omp parallel for
215#endif
216 for (int islice = 0; islice < Local_nx; islice++) {
217 for (auto && fourier_index : real_grid_result.get_fourier_range(islice, islice + 1)) {
218 auto g_fourier = real_grid_result.get_fourier_from_index(fourier_index);
219 auto f_fourier = g_fourier;
220 if (not f_and_g_are_the_same_grid)
221 f_fourier = tmp.get_fourier_from_index(fourier_index);
222 real_grid_result.set_fourier_from_index(fourier_index, f_fourier * g_fourier);
223 }
224 }
225
226 // Transform back to obtain the desired convolution
227 real_grid_result.fftw_c2r();
228 }
FML::GRID::FFTWGrid::set_fourier_from_index
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
FML::GRID::FFTWGrid::get_fourier_from_index
ComplexType get_fourier_from_index(const IndexIntType index) const
Fetch value in fourier grid by (local) index.
Definition: FFTWGrid.h:1031
FML::GRID::FFTWGrid::get_fourier_range
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

◆ fftw_c2r()

template<int N>
void FML::GRID::fftw_c2r ( FFTWGrid< N > &  in_grid,
FFTWGrid< N > &  out_grid 
)

Definition at line 1225 of file FFTWGrid.h.

1225 {
1226#ifdef DEBUG_FFTWGRID
1227 if (FML::ThisTask == 0) {
1228 std::cout << "[fftw_c2r] Transforming grid to real space\n";
1229 }
1230#endif
1231 out_grid = in_grid;
1232 out_grid.fftw_c2r();
1233 }
FML::ThisTask
int ThisTask
The MPI task number.
Definition: Global.cpp:33

◆ fftw_r2c()

template<int N>
void FML::GRID::fftw_r2c ( FFTWGrid< N > &  in_grid,
FFTWGrid< N > &  out_grid 
)

Definition at line 1236 of file FFTWGrid.h.

1236 {
1237#ifdef DEBUG_FFTWGRID
1238 if (FML::ThisTask == 0) {
1239 std::cout << "[fftw_r2c] Transforming grid to real space\n";
1240 }
1241#endif
1242 out_grid = in_grid;
1243 out_grid.fftw_r2c();
1244 }

◆ finalize_fftw()

void FML::GRID::finalize_fftw ( )

Definition at line 293 of file Global.cpp.

293 {
294#if defined(USE_FFTW) && defined(USE_MPI)
295 CLEANUP_FFTW_MPI();
296 MPI_Finalize();
297#endif
298 }

◆ init_fftw()

void FML::GRID::init_fftw ( int *  argc,
char ***  argv 
)

Definition at line 279 of file Global.cpp.

279 {
280#if defined(USE_FFTW) && defined(USE_FFTW_THREADS)
281 if (FML::MPIThreadsOK) {
282 FML::FFTWThreadsOK = INIT_FFTW_THREADS();
283 if (FML::FFTWThreadsOK) {
284 SET_FFTW_NTHREADS(FML::NThreads);
285 }
286 }
287#endif
288#if defined(USE_FFTW) && defined(USE_MPI)
289 INIT_FFTW_MPI();
290#endif
291 }
FML::MPIThreadsOK
bool MPIThreadsOK
If MPI and threads can work together with FFTW.
Definition: Global.cpp:38
FML::FFTWThreadsOK
bool FFTWThreadsOK
Definition: Global.cpp:37

◆ OPS() [1/3]

FML::GRID::OPS ( +  )

◆ OPS() [2/3]

FML::GRID::OPS ( -  )

◆ OPS() [3/3]

OPS * FML::GRID::OPS ( /  )

Definition at line 258 of file MPIGrid.h.

263 {
264 assert_mpi(n > 0, "[intlog2] Can't take a log of argument <= 0\n");
265 if (n == 1)
266 return 0;
267 int res = 0;
268 while (n > 1) {
269 n /= 2;
270 res++;
271 }
272 return res;
273 }

◆ set_fftw_nthreads()

void FML::GRID::set_fftw_nthreads ( int  nthreads)

Definition at line 300 of file Global.cpp.

300 {
301#if defined(USE_FFTW) && defined(USE_FFTW_THREADS)
302 if (FML::FFTWThreadsOK)
303 SET_FFTW_NTHREADS(nthreads);
304#endif
305 }

◆ smoothing_filter_fourier_space()

template<int N>
void FML::GRID::smoothing_filter_fourier_space ( FFTWGrid< N > &  fourier_grid,
double  smoothing_scale,
std::string  smoothing_method 
)

Low-pass filters (tophat, gaussian, sharpk)

Template Parameters
NThe dimension of the grid
Parameters
[out]fourier_gridThe fourier grid we do the smoothing of
[in]smoothing_scaleThe smoothing radius of the filter (in units of the boxsize)
[in]smoothing_methodThe smoothing filter (tophat, gaussian, sharpk)

Definition at line 55 of file SmoothingFourier.h.

57 {
58
59 // Sharp cut off kR = 1
60 std::function<double(double)> filter_sharpk = [=](double k2) -> double {
61 double kR2 = k2 * smoothing_scale * smoothing_scale;
62 if (kR2 < 1.0)
63 return 1.0;
64 return 0.0;
65 };
66 // Gaussian exp(-kR^2/2)
67 std::function<double(double)> filter_gaussian = [=](double k2) -> double {
68 double kR2 = k2 * smoothing_scale * smoothing_scale;
69 return std::exp(-0.5 * kR2);
70 };
71 // Top-hat F[ (|x| < R) ]. Implemented only for 2D and 3D
72 std::function<double(double)> filter_tophat_2D = [=](double k2) -> double {
73 double kR2 = k2 * smoothing_scale * smoothing_scale;
74 double kR = std::sqrt(kR2);
75 if (kR2 < 1e-8)
76 return 1.0;
77 return 2.0 / (kR2) * (1.0 - std::cos(kR));
78 };
79 std::function<double(double)> filter_tophat_3D = [=](double k2) -> double {
80 double kR2 = k2 * smoothing_scale * smoothing_scale;
81 double kR = std::sqrt(kR2);
82 if (kR2 < 1e-8)
83 return 1.0;
84 return 3.0 * (std::sin(kR) - kR * std::cos(kR)) / (kR2 * kR);
85 };
86
87 // Select the filter
88 std::function<double(double)> filter;
89 if (smoothing_method == "sharpk") {
90 filter = filter_sharpk;
91 } else if (smoothing_method == "gaussian") {
92 filter = filter_gaussian;
93 } else if (smoothing_method == "tophat") {
94 assert_mpi(N == 2 or N == 3,
95 "[smoothing_filter_fourier_space] Tophat filter only implemented in 2D and 3D");
96 if (N == 2)
97 filter = filter_tophat_2D;
98 if (N == 3)
99 filter = filter_tophat_3D;
100 } else {
101 throw std::runtime_error("Unknown filter " + smoothing_method + " Options: sharpk, gaussian, tophat");
102 }
103
104 // Do the smoothing
105 auto Local_nx = fourier_grid.get_local_nx();
106#ifdef USE_OMP
107#pragma omp parallel for
108#endif
109 for (int islice = 0; islice < Local_nx; islice++) {
110 [[maybe_unused]] double kmag2;
111 [[maybe_unused]] std::array<double, N> kvec;
112 for (auto && fourier_index : fourier_grid.get_fourier_range(islice, islice + 1)) {
113 fourier_grid.get_fourier_wavevector_and_norm2_by_index(fourier_index, kvec, kmag2);
114 auto value = fourier_grid.get_fourier_from_index(fourier_index);
115 value *= filter(kmag2);
116 fourier_grid.set_fourier_from_index(fourier_index, value);
117 }
118 }
119 }
FML::GRID::FFTWGrid::get_fourier_wavevector_and_norm2_by_index
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
FML::exp
const double exp
Definition: Global.h:259

◆ whitening_fourier_space()

template<int N>
void FML::GRID::whitening_fourier_space ( FFTWGrid< N > &  fourier_grid)

Take a fourier grid and divide each mode by its norm: \( f(k) \to f(k) / |f(k)| \).

Template Parameters
NThe dimension of the grid
Parameters
[out]fourier_gridThe fourier grid we do the whitening on

Definition at line 27 of file SmoothingFourier.h.

27 {
28 auto Local_nx = fourier_grid.get_local_nx();
29#ifdef USE_OMP
30#pragma omp parallel for
31#endif
32 for (int islice = 0; islice < Local_nx; islice++) {
33 [[maybe_unused]] double kmag2;
34 [[maybe_unused]] std::array<double, N> kvec;
35 for (auto && fourier_index : fourier_grid.get_fourier_range(islice, islice + 1)) {
36 auto value = fourier_grid.get_fourier_from_index(fourier_index);
37 double norm = std::sqrt(std::norm(value));
38 norm = norm == 0.0 ? 0.0 : 1.0 / norm;
39 fourier_grid.set_fourier_from_index(fourier_index, value * norm);
40 }
41 }
42 }
Generated by doxygen 1.9.3