GADGET-4
pm_periodic.h
Go to the documentation of this file.
1/*******************************************************************************
2 * \copyright This file is part of the GADGET4 N-body/SPH code developed
3 * \copyright by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
4 * \copyright (vspringel@mpa-garching.mpg.de) and all contributing authors.
5 *******************************************************************************/
6
12#ifndef PM_PERIODIC_H
13#define PM_PERIODIC_H
14
15#include "gadgetconfig.h"
16
17#include <gsl/gsl_integration.h>
18
19#if defined(PMGRID) || defined(NGENIC)
20
21#include <gsl/gsl_integration.h>
22#include <gsl/gsl_rng.h>
23#include <math.h>
24
25#include "../data/allvars.h"
26#include "../data/dtypes.h"
27#include "../data/intposconvert.h"
28#include "../data/mymalloc.h"
29#include "../data/simparticles.h"
30#include "../domain/domain.h"
31#include "../logs/timer.h"
32#include "../mpi_utils/mpi_utils.h"
33#include "../pm/pm_mpi_fft.h"
34#include "../system/system.h"
35#include "../time_integration/timestep.h"
36
37class pm_periodic : public pm_mpi_fft
38{
39 public:
40 pm_periodic(MPI_Comm comm) : setcomm(comm), pm_mpi_fft(comm) {}
41
42#if defined(PMGRID) && defined(PERIODIC)
43
44 private:
45#ifdef LONG_X_BITS
46#if PMGRID != ((PMGRID / LONG_X) * LONG_X)
47#error "PMGRID must be a multiple of the stretch factor in the x-direction"
48#endif
49#endif
50
51#ifdef LONG_Y_BITS
52#if PMGRID != ((PMGRID / LONG_Y) * LONG_Y)
53#error "PMGRID must be a multiple of the stretch factor in the y-direction"
54#endif
55#endif
56
57#ifdef LONG_Z_BITS
58#if PMGRID != ((PMGRID / LONG_Z) * LONG_Z)
59#error "PMGRID must be a multiple of the stretch factor in the x-direction"
60#endif
61#endif
62
63#define GRIDX ((PMGRID / LONG_X) * DBX + DBX_EXTRA)
64#define GRIDY ((PMGRID / LONG_Y) * DBY + DBY_EXTRA)
65#define GRIDZ ((PMGRID / LONG_Z) * DBZ + DBZ_EXTRA)
66
67#define INTCELL ((~((MyIntPosType)0)) / PMGRID + 1)
68
69#if(GRIDX > 1024) || (GRIDY > 1024) || (GRIDZ > 1024)
70 typedef long long large_array_offset; /* use a larger data type in this case so that we can always address all cells of the 3D grid
71 with a single index */
72#else
73 typedef int large_array_offset;
74#endif
75
76#ifdef NUMPART_PER_TASK_LARGE
77 typedef long long large_numpart_type; /* if there is a risk that the local particle number times 8 overflows a 32-bit integer, this
78 data type should be used */
79#else
80 typedef int large_numpart_type;
81#endif
82
83 /* variables for power spectrum estimation */
84#ifndef BINS_PS
85#define BINS_PS 4000 /* number of bins for power spectrum computation */
86#endif
87#ifndef POWERSPEC_FOLDFAC
88#define POWERSPEC_FOLDFAC 16 /* folding factor to obtain an estimate of the power spectrum on very small scales */
89#endif
90
91 char power_spec_fname[MAXLEN_PATH_EXTRA];
92
93 int NSource;
94
95 void pmforce_measure_powerspec(int flag, int *typeflag);
96 void pmforce_do_powerspec(int *typeflag);
97#if defined(GRAVITY_TALLBOX)
98 void pmforce_setup_tallbox_kernel(void);
99 double pmperiodic_tallbox_long_range_potential(double x, double y, double z);
100#endif
101
102 fft_plan myplan;
107 long long maxfftsize;
108
119 fft_real *rhogrid, *forcegrid, *workspace;
120
126 fft_complex *fft_of_rhogrid;
127
128#if defined(GRAVITY_TALLBOX)
129 fft_real *kernel;
131 fft_complex *fft_of_kernel;
132#endif
133
134#ifdef PM_ZOOM_OPTIMIZED
135
143 public:
144 struct part_slab_data
145 {
146 large_array_offset globalindex;
147 large_numpart_type partindex;
149 large_array_offset localindex;
151 };
152 part_slab_data *part;
154 /* realize the comparison function as a functor, so that it can have an internal state (here the data array for which we sort indices
155 */
156 struct pm_periodic_sortindex_comparator
157 {
158 private:
159 const part_slab_data *data;
160
161 public:
162 pm_periodic_sortindex_comparator(const part_slab_data *data_) : data(data_) {}
163
164 bool operator()(const large_numpart_type &a, const large_numpart_type &b) const
165 {
166 return data[a].globalindex < data[b].globalindex;
167 }
168 };
169
170 private:
171 size_t *localfield_sendcount, *localfield_first, *localfield_offset, *localfield_recvcount;
172 large_array_offset *localfield_globalindex, *import_globalindex;
173 fft_real *localfield_data, *import_data;
174
175 void pmforce_zoom_optimized_prepare_density(int mode, int *typelist);
176 void pmforce_zoom_optimized_readout_forces_or_potential(fft_real *grid, int dim);
177
178#else
179
180 struct partbuf
181 {
182#ifndef LEAN
183 MyFloat Mass;
184#endif
185 MyIntPosType IntPos[3];
186 };
187 partbuf *partin, *partout;
188
189 size_t nimport, nexport;
190
191 size_t *Sndpm_count, *Sndpm_offset;
192 size_t *Rcvpm_count, *Rcvpm_offset;
193
194 void pmforce_uniform_optimized_prepare_density(int mode, int *typelist);
195
196 void pmforce_uniform_optimized_readout_forces_or_potential_xy(fft_real *grid, int dim);
197
198 void pmforce_uniform_optimized_readout_forces_or_potential_xz(fft_real *grid, int dim);
199 void pmforce_uniform_optimized_readout_forces_or_potential_zy(fft_real *grid, int dim);
200#endif
201
202 public:
203 simparticles *Sp;
204
205 void pm_init_periodic(simparticles *Sp_ptr);
206 void pmforce_periodic(int mode, int *typelist);
207
208 void calculate_power_spectra(int num);
209
210 static double growthfactor_integrand(double a, void *param)
211 {
212 return pow(a / (All.Omega0 + (1 - All.Omega0 - All.OmegaLambda) * a + All.OmegaLambda * a * a * a), 1.5);
213 }
214
215 double linear_growth_factor(double astart, double aend) { return linear_growth(aend) / linear_growth(astart); }
216
217 double linear_growth(double a)
218 {
219 double hubble_a = sqrt(All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda);
220
221 const int worksize = 100000;
222
223 double result, abserr;
224 gsl_function F;
225
226 gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(worksize);
227 F.function = &growthfactor_integrand;
228
229 gsl_integration_qag(&F, 0, a, 0, 1.0e-8, worksize, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
230
231 gsl_integration_workspace_free(workspace);
232
233 return hubble_a * result;
234 }
235#endif
236};
237
238#endif
239
240#endif
global_data_all_processes All
Definition: main.cc:40
#define MAXLEN_PATH_EXTRA
Definition: constants.h:301
float MyFloat
Definition: dtypes.h:86
uint32_t MyIntPosType
Definition: dtypes.h:35
expr pow(half base, half exp)
Definition: half.hpp:2803
expr sqrt(half arg)
Definition: half.hpp:2777