Line data Source code
1 : #include "crpropa/magneticField/turbulentField/GridTurbulence.h"
2 : #include "crpropa/GridTools.h"
3 : #include "crpropa/Random.h"
4 :
5 : #ifdef CRPROPA_HAVE_FFTW3F
6 :
7 : namespace crpropa {
8 :
9 :
10 6 : GridTurbulence::GridTurbulence(const TurbulenceSpectrum &spectrum,
11 : const GridProperties &gridProp,
12 6 : unsigned int seed)
13 6 : : TurbulentField(spectrum), seed(seed) {
14 6 : initGrid(gridProp);
15 6 : checkGridRequirements(gridPtr, spectrum.getLmin(), spectrum.getLmax());
16 6 : initTurbulence();
17 6 : }
18 :
19 6 : void GridTurbulence::initGrid(const GridProperties &p) {
20 6 : gridPtr = new Grid3f(p);
21 6 : }
22 :
23 4 : Vector3d GridTurbulence::getField(const Vector3d &pos) const {
24 4 : return gridPtr->interpolate(pos);
25 : }
26 :
27 1 : const ref_ptr<Grid3f> &GridTurbulence::getGrid() const { return gridPtr; }
28 :
29 6 : void GridTurbulence::initTurbulence() {
30 :
31 : Vector3d spacing = gridPtr->getSpacing();
32 : size_t n = gridPtr->getNx(); // size of array
33 6 : size_t n2 = (size_t)floor(n / 2) +
34 6 : 1; // size array in z-direction in configuration space
35 :
36 : // arrays to hold the complex vector components of the B(k)-field
37 : fftwf_complex *Bkx, *Bky, *Bkz;
38 6 : Bkx = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
39 6 : Bky = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
40 6 : Bkz = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
41 :
42 6 : Random random;
43 6 : if (seed != 0)
44 4 : random.seed(seed); // use given seed
45 :
46 : // calculate the n possible discrete wave numbers
47 6 : double K[n];
48 390 : for (size_t i = 0; i < n; i++)
49 384 : K[i] = ((double)i / n - i / (n / 2));
50 :
51 : // double kMin = 2*M_PI / lMax; // * 2 * spacing.x; // spacing.x / lMax;
52 : // double kMax = 2*M_PI / lMin; // * 2 * spacing.x; // spacing.x / lMin;
53 6 : double kMin = spacing.x / spectrum.getLmax();
54 6 : double kMax = spacing.x / spectrum.getLmin();
55 6 : auto lambda = 1 / spacing.x * 2 * M_PI;
56 :
57 : Vector3f n0(1, 1, 1); // arbitrary vector to construct orthogonal base
58 :
59 390 : for (size_t ix = 0; ix < n; ix++) {
60 24960 : for (size_t iy = 0; iy < n; iy++) {
61 835584 : for (size_t iz = 0; iz < n2; iz++) {
62 :
63 : Vector3f ek, e1, e2; // orthogonal base
64 :
65 811008 : size_t i = ix * n * n2 + iy * n2 + iz;
66 811008 : ek.setXYZ(K[ix], K[iy], K[iz]);
67 811008 : double k = ek.getR();
68 :
69 : // wave outside of turbulent range -> B(k) = 0
70 811008 : if ((k < kMin) || (k > kMax)) {
71 395939 : Bkx[i][0] = 0;
72 395939 : Bkx[i][1] = 0;
73 395939 : Bky[i][0] = 0;
74 395939 : Bky[i][1] = 0;
75 395939 : Bkz[i][0] = 0;
76 395939 : Bkz[i][1] = 0;
77 : continue;
78 : }
79 :
80 : // construct an orthogonal base ek, e1, e2
81 415069 : if (ek.isParallelTo(n0, float(1e-3))) {
82 : // ek parallel to (1,1,1)
83 : e1.setXYZ(-1., 1., 0);
84 : e2.setXYZ(1., 1., -2.);
85 : } else {
86 : // ek not parallel to (1,1,1)
87 : e1 = n0.cross(ek);
88 : e2 = ek.cross(e1);
89 : }
90 : e1 /= e1.getR();
91 : e2 /= e2.getR();
92 :
93 : // random orientation perpendicular to k
94 415069 : double theta = 2 * M_PI * random.rand();
95 415069 : Vector3f b = e1 * std::cos(theta) + e2 * std::sin(theta); // real b-field vector
96 :
97 : // normal distributed amplitude with mean = 0
98 415069 : b *= std::sqrt(spectrum.energySpectrum(k*lambda));
99 :
100 : // uniform random phase
101 415069 : double phase = 2 * M_PI * random.rand();
102 415069 : double cosPhase = std::cos(phase); // real part
103 415069 : double sinPhase = std::sin(phase); // imaginary part
104 :
105 415069 : Bkx[i][0] = b.x * cosPhase;
106 415069 : Bkx[i][1] = b.x * sinPhase;
107 415069 : Bky[i][0] = b.y * cosPhase;
108 415069 : Bky[i][1] = b.y * sinPhase;
109 415069 : Bkz[i][0] = b.z * cosPhase;
110 415069 : Bkz[i][1] = b.z * sinPhase;
111 : } // for iz
112 : } // for iy
113 : } // for ix
114 :
115 6 : executeInverseFFTInplace(gridPtr, Bkx, Bky, Bkz);
116 :
117 6 : fftwf_free(Bkx);
118 6 : fftwf_free(Bky);
119 6 : fftwf_free(Bkz);
120 :
121 24 : scaleGrid(gridPtr, spectrum.getBrms() /
122 12 : rmsFieldStrength(gridPtr)); // normalize to Brms
123 12 : }
124 :
125 : // Check the grid properties before the FFT procedure
126 13 : void GridTurbulence::checkGridRequirements(ref_ptr<Grid3f> grid, double lMin,
127 : double lMax) {
128 : size_t Nx = grid->getNx();
129 : size_t Ny = grid->getNy();
130 : size_t Nz = grid->getNz();
131 : Vector3d spacing = grid->getSpacing();
132 :
133 13 : if ((Nx != Ny) or (Ny != Nz))
134 0 : throw std::runtime_error("turbulentField: only cubic grid supported");
135 13 : if ((spacing.x != spacing.y) or (spacing.y != spacing.z))
136 0 : throw std::runtime_error("turbulentField: only equal spacing suported");
137 13 : if (lMin < 2 * spacing.x)
138 1 : throw std::runtime_error("turbulentField: lMin < 2 * spacing");
139 12 : if (lMax > Nx * spacing.x) // before was (lMax > Nx * spacing.x / 2), why?
140 1 : throw std::runtime_error("turbulentField: lMax > size");
141 11 : if (lMax < lMin)
142 1 : throw std::runtime_error("lMax < lMin");
143 10 : }
144 :
145 : // Execute inverse discrete FFT in-place for a 3D grid, from complex to real
146 : // space
147 10 : void GridTurbulence::executeInverseFFTInplace(ref_ptr<Grid3f> grid,
148 : fftwf_complex *Bkx,
149 : fftwf_complex *Bky,
150 : fftwf_complex *Bkz) {
151 :
152 : size_t n = grid->getNx(); // size of array
153 10 : size_t n2 = (size_t)floor(n / 2) +
154 10 : 1; // size array in z-direction in configuration space
155 :
156 : // in-place, complex to real, inverse Fourier transformation on each
157 : // component note that the last elements of B(x) are unused now
158 : float *Bx = (float *)Bkx;
159 10 : fftwf_plan plan_x = fftwf_plan_dft_c2r_3d(n, n, n, Bkx, Bx, FFTW_ESTIMATE);
160 10 : fftwf_execute(plan_x);
161 10 : fftwf_destroy_plan(plan_x);
162 :
163 : float *By = (float *)Bky;
164 10 : fftwf_plan plan_y = fftwf_plan_dft_c2r_3d(n, n, n, Bky, By, FFTW_ESTIMATE);
165 10 : fftwf_execute(plan_y);
166 10 : fftwf_destroy_plan(plan_y);
167 :
168 : float *Bz = (float *)Bkz;
169 10 : fftwf_plan plan_z = fftwf_plan_dft_c2r_3d(n, n, n, Bkz, Bz, FFTW_ESTIMATE);
170 10 : fftwf_execute(plan_z);
171 10 : fftwf_destroy_plan(plan_z);
172 :
173 : // save to grid
174 650 : for (size_t ix = 0; ix < n; ix++) {
175 41600 : for (size_t iy = 0; iy < n; iy++) {
176 2662400 : for (size_t iz = 0; iz < n; iz++) {
177 2621440 : size_t i = ix * n * 2 * n2 + iy * 2 * n2 + iz;
178 : Vector3f &b = grid->get(ix, iy, iz);
179 2621440 : b.x = Bx[i];
180 2621440 : b.y = By[i];
181 2621440 : b.z = Bz[i];
182 : }
183 : }
184 : }
185 10 : }
186 :
187 0 : Vector3f GridTurbulence::getMeanFieldVector() const {
188 0 : return meanFieldVector(gridPtr);
189 : }
190 :
191 0 : double GridTurbulence::getMeanFieldStrength() const {
192 0 : return meanFieldStrength(gridPtr);
193 : }
194 :
195 0 : double GridTurbulence::getRmsFieldStrength() const {
196 0 : return rmsFieldStrength(gridPtr);
197 : }
198 :
199 0 : std::array<float, 3> GridTurbulence::getRmsFieldStrengthPerAxis() const {
200 0 : return rmsFieldStrengthPerAxis(gridPtr);
201 : }
202 :
203 0 : std::vector<std::pair<int, float>> GridTurbulence::getPowerSpectrum() const {
204 0 : return gridPowerSpectrum(gridPtr);
205 : }
206 :
207 0 : void GridTurbulence::dumpToFile(std::string filename) const {
208 0 : dumpGrid(gridPtr, filename);
209 0 : }
210 :
211 : } // namespace crpropa
212 :
213 : #endif // CRPROPA_HAVE_FFTW3F
|
