MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
maiafftw.h
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7// Copyright (C) 2024 The m-AIA AUTHORS
8//
9// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
10//
11// SPDX-License-Identifier: LGPL-3.0-only
12
13// Copyright (C) 2024 The m-AIA AUTHORS
14//
15// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
16//
17// SPDX-License-Identifier: LGPL-3.0-only
18
19
20// Copyright (C) 2024 The m-AIA AUTHORS
21//
22// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
23//
24// SPDX-License-Identifier: LGPL-3.0-only
25
26
27// Copyright (C) 2024 The m-AIA AUTHORS
28//
29// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
30//
31// SPDX-License-Identifier: LGPL-3.0-only
32
33
34// Copyright (C) 2019 The m-AIA AUTHORS
35//
36// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
37//
38// SPDX-License-Identifier: LGPL-3.0-only
39
40
41#ifndef MAIA_FFTW_H
42#define MAIA_FFTW_H
43
44#include <complex>
45#include <fftw3-mpi.h>
46#include <fftw3.h>
47#include <random>
48
49#include "../COMM/mpioverride.h"
50#include "../INCLUDE/maiaconstants.h"
51#include "../INCLUDE/maiatypes.h"
52#include "../IO/context.h"
53#include "../MEMORY/scratch.h"
54#include "debug.h"
55
56namespace maia {
57namespace math {
58
59inline MInt getGlobalPosFFTW(MInt i0, MInt i1, MInt i2, MInt ny, MInt nz) { return i2 + nz * (i1 + ny * i0); }
60
61inline MFloat triadImag(fftw_complex& a, fftw_complex& b, fftw_complex& c) {
62 return (a[0] * b[1] * c[0] + a[1] * b[0] * c[0] - a[0] * b[0] * c[1] + a[1] * b[1] * c[1]);
63}
64
65
83inline MUlong initFft(fftw_complex*& uPhysField, fftw_complex*& nabla2P, MInt nx, MInt ny, MInt nz,
84 const MFloat kpRatio, const MInt spectrumId, MIntScratchSpace& fftInfo, const MPI_Comm comm,
85 const MBool computeNabla2P = false) {
86 DEBUG("initFft entry", MAIA_DEBUG_TRACE_IN);
87
88 const MInt size = nx * ny * nz;
89 const ptrdiff_t rank = 3;
90 const MInt howmany = 3;
91 if(nx % 2 != 0 || ny % 2 != 0 || nz % 2 != 0) {
92 std::stringstream errorMessage;
93 errorMessage << " FFTInit: domainsize must NOT be an odd number! " << nx << " " << ny << " " << nz << std::endl;
94 mTerm(1, AT_, errorMessage.str());
95 }
96
97 m_log << " --- initializing FFTW --- " << std::endl;
98 m_log << " domain size = " << nx << "x" << ny << "x" << nz << std::endl;
99
100 MInt domainId;
101 MInt noDomains;
102 MPI_Comm_rank(comm, &domainId);
103 MPI_Comm_size(comm, &noDomains);
104 MPI_Comm MPI_COMM_FFTW;
105 MInt maxRank = 1;
106 maxRank = mMin(nz, noDomains);
107 while(nz % maxRank != 0) {
108 maxRank--;
109 }
110
111 MInt color = (domainId < maxRank) ? 0 : MPI_UNDEFINED;
112 MPI_Comm_split(comm, color, domainId, &MPI_COMM_FFTW, AT_, "MPI_COMM_FFTW");
113 m_log << "no ranks for fft: " << maxRank << std::endl;
114
115 MInt tmpHowMany = 2;
116 ptrdiff_t alloc_local, local_n0, local_0_start;
117 ptrdiff_t alloc_localP, local_n0P, local_0_startP;
118 ptrdiff_t alloc_localTmp, tmpLocal_n0, tmpLocal_0_start;
119
120 const ptrdiff_t n[3] = {nx, ny, nz};
121#ifndef MAIA_WINDOWS
122 alloc_local = (domainId < maxRank) ? fftw_mpi_local_size_many(rank, n, howmany, FFTW_MPI_DEFAULT_BLOCK, MPI_COMM_FFTW,
123 &local_n0, &local_0_start)
124 : 1;
125 alloc_localTmp = (domainId < maxRank) ? fftw_mpi_local_size_many(rank, n, tmpHowMany, FFTW_MPI_DEFAULT_BLOCK,
126 MPI_COMM_FFTW, &tmpLocal_n0, &tmpLocal_0_start)
127 : 1;
128 alloc_localP = (domainId < maxRank) ? fftw_mpi_local_size_many(rank, n, 1, FFTW_MPI_DEFAULT_BLOCK, MPI_COMM_FFTW,
129 &local_n0P, &local_0_startP)
130 : 1;
131#endif
132
133 fftInfo[0] = maxRank;
134 fftInfo[1] = (domainId < maxRank) ? local_n0 : 0;
135 fftInfo[2] = (domainId < maxRank) ? local_0_start : size;
136 fftInfo[3] = (domainId < maxRank) ? alloc_local : 3;
137
150 static const MInt referenceCubeSize =
151 ((Context::propertyExists("referenceCubeSize")) ? Context::getBasicProperty<MInt>("referenceCubeSize", AT_) : 1);
152
153 static const MInt kMinSpec = ((Context::propertyExists("kMin")) ? Context::getBasicProperty<MInt>("kMin", AT_) : 0);
154
155 static const MInt kMaxSpec = ((Context::propertyExists("kMax")) ? Context::getBasicProperty<MInt>("kMax", AT_) : 0);
156
157 static const int64_t seed0 = (int64_t)(
158 (Context::propertyExists("randomDeviceSeed")) ? Context::getBasicProperty<MFloat>("randomDeviceSeed", AT_) : -1);
159
160 static MUlong seed = (unsigned)((seed0 > -1) ? seed0 : std::random_device()());
161 MPI_Bcast(&seed, 1, MPI_UNSIGNED_LONG, 0, comm, AT_, "seed");
162 std::mt19937_64 gen(seed);
163 std::normal_distribution<> distr(F0, F1);
164
165 uPhysField = (fftw_complex*)fftw_malloc(alloc_local * sizeof(fftw_complex));
166
167 if(computeNabla2P) nabla2P = (fftw_complex*)fftw_malloc(alloc_localP * sizeof(fftw_complex));
168
169 if(domainId < maxRank) {
170 ptrdiff_t rSlice_0_start = nx - (MInt)local_0_start - (MInt)local_n0 + 1;
171 ScratchSpace<fftw_complex> uHatField(alloc_local, AT_, "uHatField");
172 ScratchSpace<fftw_complex> vHatField(alloc_local, AT_, "vHatField");
173 ScratchSpace<fftw_complex> wHatField(alloc_local, AT_, "wHatField");
174 ScratchSpace<fftw_complex> ruHatField(alloc_local, AT_, "ruHatField");
175 ScratchSpace<fftw_complex> rvHatField(alloc_local, AT_, "rvHatField");
176 ScratchSpace<fftw_complex> rwHatField(alloc_local, AT_, "rwHatField");
177 ScratchSpace<fftw_complex> tmpHatFieldMem(alloc_localTmp, AT_, "tmpHatFieldMem");
178
179 fftw_complex* hatField = uPhysField;
180 fftw_complex* tmpHat = &tmpHatFieldMem[0];
181
182#ifndef MAIA_WINDOWS
183 fftw_plan plan = fftw_mpi_plan_many_dft(rank, n, howmany, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, hatField,
184 hatField, MPI_COMM_FFTW, FFTW_BACKWARD, FFTW_ESTIMATE);
185 fftw_plan tmpPlan = fftw_mpi_plan_many_dft(rank, n, tmpHowMany, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
186 tmpHat, tmpHat, MPI_COMM_FFTW, FFTW_BACKWARD, FFTW_ESTIMATE);
187#else
188 fftw_plan plan;
189 fftw_plan tmpPlan;
190#endif
191
192 m_log << "random device seed is " << seed << std::endl;
193 const MInt refSize =
194 (nx <= 256) ? 256 : nx; // This won´t work if you want to compare LES to grids with more than 256^3 cells
195 for(MInt k0 = -refSize / 2 + 1; k0 <= refSize / 2; k0++) {
196 for(MInt k1 = -refSize / 2 + 1; k1 <= refSize / 2; k1++) {
197 for(MInt k2 = -refSize / 2 + 1; k2 <= refSize / 2; k2++) {
198 if(k0 < -nx / 2 + 1 || k0 > nx / 2 || k1 < -ny / 2 + 1 || k1 > ny / 2 || k2 < -nz / 2 + 1 || k2 > nz / 2) {
199 (void)distr(gen);
200 (void)distr(gen);
201 (void)distr(gen);
202 (void)distr(gen);
203 (void)distr(gen);
204 (void)distr(gen);
205 if(nx == refSize) {
206 mTerm(1, AT_, "wrong refSize.");
207 }
208 continue;
209 }
210
211 MInt i0 = (k0 >= 0) ? k0 : k0 + nx;
212 MInt i1 = (k1 >= 0) ? k1 : k1 + ny;
213 MInt i2 = (k2 >= 0) ? k2 : k2 + nz;
214 MInt pos = getGlobalPosFFTW(i0, i1, i2, ny, nz);
215 MInt localrPos = pos - (((MInt)rSlice_0_start) * ny * nz);
216 MInt localPos = pos - (((MInt)local_0_start) * ny * nz);
217
218 // Collect random numbers on the current FFTW slice
219 if(pos >= ((MInt)local_0_start) * nz * ny && pos < ((MInt)(local_0_start + local_n0) * ny * nz)) {
220 if(localPos >= ny * nz * local_n0 || localPos < 0) {
221 mTerm(1, AT_, "index exceeds array(1a)");
222 }
223 uHatField[localPos][0] = distr(gen);
224 uHatField[localPos][1] = distr(gen);
225 vHatField[localPos][0] = distr(gen);
226 vHatField[localPos][1] = distr(gen);
227 wHatField[localPos][0] = distr(gen);
228 wHatField[localPos][1] = distr(gen);
229 continue;
230 }
231 if(pos >= ((MInt)rSlice_0_start) * nz * ny && pos < ((MInt)(rSlice_0_start + local_n0) * ny * nz) && i0 != 0
232 && i0 != nx / 2) {
233 localrPos = pos - (((MInt)rSlice_0_start) * ny * nz);
234 if(localrPos >= ny * nz * local_n0 || localrPos < 0) {
235 mTerm(1, AT_, "index exceeds array(1c)");
236 }
237 ruHatField[localrPos][0] = distr(gen);
238 ruHatField[localrPos][1] = distr(gen);
239 rvHatField[localrPos][0] = distr(gen);
240 rvHatField[localrPos][1] = distr(gen);
241 rwHatField[localrPos][0] = distr(gen);
242 rwHatField[localrPos][1] = distr(gen);
243 continue;
244 }
245 if(!(pos >= ((MInt)local_0_start) * nz * ny && pos < ((MInt)(local_0_start + local_n0) * ny * nz))
246 && !(pos >= ((MInt)rSlice_0_start) * nz * ny && pos < ((MInt)(rSlice_0_start + local_n0) * ny * nz)
247 && i0 != 0 && i0 != nx / 2)) {
248 for(MInt k = 0; k < howmany; k++) {
249 (void)distr(gen);
250 (void)distr(gen);
251 }
252 }
253 }
254 }
255 }
256
257 const MFloat kMin = F2 * PI / referenceCubeSize;
258 const MFloat kp = kpRatio * kMin * referenceCubeSize;
259
260 MFloat eng = F0;
261 MFloat eps = F0;
262 MFloat length = F0;
263
264 for(MInt i0 = (MInt)local_0_start; i0 < (MInt)(local_0_start + local_n0); i0++) {
265 for(MInt i1 = 0; i1 < ny; i1++) {
266 for(MInt i2 = 0; i2 < nz; i2++) {
267 MFloat k[3];
268 MInt j0 = (i0 > nx / 2) ? i0 - nx : i0;
269 MInt j1 = (i1 > ny / 2) ? i1 - ny : i1;
270 MInt j2 = (i2 > nz / 2) ? i2 - nz : i2;
271
272 /*To limit the spectrum to a given number of wavenumbers.
273 kMinSpec and kMaxSpec are integers. Tested with spectrumId = 3*/
274 if(kMinSpec > 0 && kMaxSpec > 0) {
275 if(j0 > kMaxSpec || j0 < kMinSpec) j0 = 0;
276 if(j1 > kMaxSpec || j1 < kMinSpec) j1 = 0;
277 if(j2 > kMaxSpec || j2 < kMinSpec) j2 = 0;
278 }
279 k[0] = ((MFloat)j0) * kMin;
280 k[1] = ((MFloat)j1) * kMin;
281 k[2] = ((MFloat)j2) * kMin;
282
283 MInt pos = getGlobalPosFFTW(i0, i1, i2, ny, nz);
284 MInt rpos = ((nz - i2) % nz) + nz * (((ny - i1) % ny) + ny * ((nx - i0) % nx));
285 MInt localPos = pos - (((MInt)local_0_start) * ny * nz);
286 MInt localrPos = rpos - (((MInt)rSlice_0_start) * ny * nz);
287 ASSERT(pos > -1 && pos < size && rpos > -1 && rpos < size, "");
288
289 MFloat kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
290
291 MFloat energy = F0;
292
293 if(kAbs < 1e-15) {
294 hatField[howmany * localPos][0] = F0;
295 hatField[howmany * localPos][1] = F0;
296 hatField[howmany * localPos + 1][0] = F0;
297 hatField[howmany * localPos + 1][1] = F0;
298 hatField[howmany * localPos + 2][0] = F0;
299 hatField[howmany * localPos + 2][1] = F0;
300 continue;
301 }
302
303 // set spectral distribution
304 switch(spectrumId) {
305 case 0: {
306 energy = pow(kAbs / kp, 4.0) * exp(-2.0 * (kAbs / kp) * (kAbs / kp));
307 break;
308 }
309 case 1: {
310 energy = POW2(kAbs / kp) * exp(-F1B2 * POW2(kAbs / kp));
311 break;
312 }
313 case 2: {
314 // see Schumann and Patterson
315 energy = F3B2 * (kAbs / POW2(kp)) * exp(-kAbs / kp);
316 break;
317 }
318 default: {
319 mTerm(1, AT_, "Unknown spectrum.");
320 break;
321 }
322 }
323
324 if(j0 > 0 && j1 == 0 && j2 == 0) {
325 const MFloat dk = kMin;
326 eng += dk * energy;
327 eps += F2 * dk * POW2(k[0]) * energy;
328 length += F1B2 * PI * dk * energy / k[0];
329 }
330
331 // assemble uHat; see Orszag, Phys. Fluids (1969)
332 MFloat r[3], s[3];
333 const MFloat fac = PI * sqrt(energy) / (sqrt(referenceCubeSize * F2) * referenceCubeSize * kAbs);
334
335 if(i0 == 0 || i0 == nx / 2) {
336 localrPos = rpos - (((MInt)local_0_start) * ny * nz);
337 r[0] = fac * (uHatField[localPos][0] + uHatField[localrPos][0]);
338 s[0] = fac * (uHatField[localPos][1] - uHatField[localrPos][1]);
339 r[1] = fac * (vHatField[localPos][0] + vHatField[localrPos][0]);
340 s[1] = fac * (vHatField[localPos][1] - vHatField[localrPos][1]);
341 r[2] = fac * (wHatField[localPos][0] + wHatField[localrPos][0]);
342 s[2] = fac * (wHatField[localPos][1] - wHatField[localrPos][1]);
343 } else {
344 if(localPos >= ny * nz * local_n0 || localrPos >= ny * nz * local_n0 || localPos < 0 || localrPos < 0) {
345 mTerm(1, AT_, "index exceeds array(2)");
346 }
347 r[0] = fac * (uHatField[localPos][0] + ruHatField[localrPos][0]);
348 s[0] = fac * (uHatField[localPos][1] - ruHatField[localrPos][1]);
349 r[1] = fac * (vHatField[localPos][0] + rvHatField[localrPos][0]);
350 s[1] = fac * (vHatField[localPos][1] - rvHatField[localrPos][1]);
351 r[2] = fac * (wHatField[localPos][0] + rwHatField[localrPos][0]);
352 s[2] = fac * (wHatField[localPos][1] - rwHatField[localrPos][1]);
353 }
354
355
356 hatField[howmany * localPos][0] = (1.0 - k[0] * k[0] / (kAbs * kAbs)) * r[0]
357 - k[0] * k[1] / (kAbs * kAbs) * r[1] - k[0] * k[2] / (kAbs * kAbs) * r[2];
358 hatField[howmany * localPos][1] = (1.0 - k[0] * k[0] / (kAbs * kAbs)) * s[0]
359 - k[0] * k[1] / (kAbs * kAbs) * s[1] - k[0] * k[2] / (kAbs * kAbs) * s[2];
360
361 hatField[howmany * localPos + 1][0] = -k[1] * k[0] / (kAbs * kAbs) * r[0]
362 + (1.0 - k[1] * k[1] / (kAbs * kAbs)) * r[1]
363 - k[1] * k[2] / (kAbs * kAbs) * r[2];
364 hatField[howmany * localPos + 1][1] = -k[1] * k[0] / (kAbs * kAbs) * s[0]
365 + (1.0 - k[1] * k[1] / (kAbs * kAbs)) * s[1]
366 - k[1] * k[2] / (kAbs * kAbs) * s[2];
367
368 hatField[howmany * localPos + 2][0] = -k[2] * k[0] / (kAbs * kAbs) * r[0] - k[2] * k[1] / (kAbs * kAbs) * r[1]
369 + (1.0 - k[2] * k[2] / (kAbs * kAbs)) * r[2];
370 hatField[howmany * localPos + 2][1] = -k[2] * k[0] / (kAbs * kAbs) * s[0] - k[2] * k[1] / (kAbs * kAbs) * s[1]
371 + (1.0 - k[2] * k[2] / (kAbs * kAbs)) * s[2];
372 }
373 }
374 }
375 if(computeNabla2P) {
376 for(MInt i = 0; i < local_n0 * nz * ny; i++) {
377 if(i >= alloc_local / 3) mTerm(1, AT_, "index exceeds array(6)");
378 nabla2P[i][0] = F0;
379 nabla2P[i][1] = F0;
380 }
381
382
383 for(MInt i = 0; i < 3; i++) {
384 for(MInt j = i; j < 3; j++) {
385 for(MInt i0 = (MInt)tmpLocal_0_start; i0 < (MInt)(tmpLocal_0_start + tmpLocal_n0); i0++) {
386 for(MInt i1 = 0; i1 < ny; i1++) {
387 for(MInt i2 = 0; i2 < nz; i2++) {
388 MFloat k[3];
389 MInt j0 = (i0 > nx / 2) ? i0 - nx : i0;
390 MInt j1 = (i1 > ny / 2) ? i1 - ny : i1;
391 MInt j2 = (i2 > nz / 2) ? i2 - nz : i2;
392 k[0] = ((MFloat)j0) * kMin;
393 k[1] = ((MFloat)j1) * kMin;
394 k[2] = ((MFloat)j2) * kMin;
395 MInt pos = getGlobalPosFFTW(i0, i1, i2, ny, nz);
396 MInt localPos = pos - (((MInt)tmpLocal_0_start) * ny * nz);
397 if(localPos * tmpHowMany + tmpHowMany - 1 >= alloc_localTmp) mTerm(1, AT_, "index exceeds array(3)");
398 const MInt dimIndex0 = (i == 0) ? 0 : ((i == 1) ? 1 : 2);
399 const MInt dimIndex1 = (j == 0) ? 0 : ((j == 1) ? 1 : 2);
400 tmpHat[localPos * tmpHowMany][0] = -k[j] * hatField[howmany * localPos + dimIndex0][1];
401 tmpHat[localPos * tmpHowMany][1] = k[j] * hatField[howmany * localPos + dimIndex0][0];
402 tmpHat[localPos * tmpHowMany + 1][0] = -k[i] * hatField[howmany * localPos + dimIndex1][1];
403 tmpHat[localPos * tmpHowMany + 1][1] = k[i] * hatField[howmany * localPos + dimIndex1][0];
404 }
405 }
406 }
407
408 fftw_execute(tmpPlan);
409 MFloat fac = (i == j) ? F1 : F2;
410
411 for(MInt i0 = (MInt)tmpLocal_0_start; i0 < (MInt)(tmpLocal_0_start + tmpLocal_n0); i0++) {
412 for(MInt i1 = 0; i1 < ny; i1++) {
413 for(MInt i2 = 0; i2 < nz; i2++) {
414 MInt pos = getGlobalPosFFTW(i0, i1, i2, ny, nz);
415 MInt localPos = pos - (((MInt)tmpLocal_0_start) * ny * nz);
416 if(localPos * tmpHowMany + tmpHowMany - 1 >= alloc_localTmp) mTerm(1, AT_, "index exceeds array(4)");
417 if(localPos >= alloc_local / 3) mTerm(1, AT_, "index exceeds array(5)");
418 nabla2P[localPos][0] += fac * tmpHat[localPos * tmpHowMany][0] * tmpHat[localPos * tmpHowMany + 1][0];
419 tmpHat[localPos * tmpHowMany][0] = 0;
420 tmpHat[localPos * tmpHowMany][1] = 0;
421 tmpHat[localPos * tmpHowMany + 1][0] = 0;
422 tmpHat[localPos * tmpHowMany + 1][1] = 0;
423 }
424 }
425 }
426 }
427 }
428
429 fftw_destroy_plan(tmpPlan);
430
431 for(MInt i = 0; i < local_n0 * nz * ny; i++) {
432 if(i >= alloc_local / 3) mTerm(1, AT_, "index exceeds array(6)");
433 nabla2P[i][1] = F0;
434 }
435
436 fftw_plan planP = fftw_mpi_plan_many_dft(rank, n, 1, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, nabla2P,
437 nabla2P, MPI_COMM_FFTW, FFTW_FORWARD, FFTW_ESTIMATE);
438 fftw_execute(planP);
439
440 for(MInt k = 0; k < local_n0 * nz * ny; k++) {
441 if(k >= alloc_local / 3) mTerm(1, AT_, "index exceeds array(7)");
442 nabla2P[k][0] /= ((MFloat)size);
443 nabla2P[k][1] /= ((MFloat)size);
444 }
445
446 fftw_destroy_plan(planP);
447 for(MInt i0 = (MInt)local_0_start; i0 < (MInt)(local_0_start + local_n0); i0++) {
448 for(MInt i1 = 0; i1 < ny; i1++) {
449 for(MInt i2 = 0; i2 < nz; i2++) {
450 MFloat k[3];
451 MInt j0 = (i0 > nx / 2) ? i0 - nx : i0;
452 MInt j1 = (i1 > ny / 2) ? i1 - ny : i1;
453 MInt j2 = (i2 > nz / 2) ? i2 - nz : i2;
454 k[0] = ((MFloat)j0) * kMin;
455 k[1] = ((MFloat)j1) * kMin;
456 k[2] = ((MFloat)j2) * kMin;
457 MInt pos = getGlobalPosFFTW(i0, i1, i2, ny, nz);
458 MInt localPos = pos - (((MInt)local_0_start) * ny * nz);
459 MFloat k2 = k[0] * k[0] + k[1] * k[1] + k[2] * k[2];
460 if(localPos >= alloc_local / 3) mTerm(1, AT_, "index exceeds array(8)");
461 if(k2 < 1e-15) {
462 nabla2P[localPos][0] = F0;
463 nabla2P[localPos][1] = F0;
464 continue;
465 }
466 nabla2P[localPos][0] = nabla2P[localPos][0] / k2;
467 nabla2P[localPos][1] = nabla2P[localPos][1] / k2;
468 }
469 }
470 }
471
472#ifndef MAIA_WINDOWS
473 fftw_plan planP2 = fftw_mpi_plan_many_dft(rank, n, 1, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, nabla2P,
474 nabla2P, MPI_COMM_FFTW, FFTW_BACKWARD, FFTW_ESTIMATE);
475 fftw_execute(planP2);
476 fftw_destroy_plan(planP2);
477#endif
478 }
479 MPI_Allreduce(MPI_IN_PLACE, &eng, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE", "eng");
480 MPI_Allreduce(MPI_IN_PLACE, &eps, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE", "eps");
481 MPI_Allreduce(MPI_IN_PLACE, &length, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE", "length");
482
483 if(domainId == 0)
484 std::cerr << "SPECTRUM: kinetic energy: " << eng << " (theoretical: 1.5)"
485 << ", dissipation rate: " << eps << " (theoretical: " << 72.0 * POW2(kpRatio * PI) << ")"
486 << ", integral length/unit length: " << length << " (theoretical: " << F3B8 / kpRatio << ")"
487 << std::endl;
488
489 fftw_execute(plan);
490 }
491 m_log << "fft finished" << std::endl;
492
493 return seed;
494
495 DEBUG("initFft return", MAIA_DEBUG_TRACE_OUT);
496}
497
498/* brief returns a normal distributed random-
499 * number with mu=mean and sigma=standard deviation
500 *
501 */
502inline MFloat randnormal(MFloat mu, MFloat sigma) {
503 TRACE();
504
505 static MBool deviateAvailable = false; // flag
506 static float storedDeviate; // deviate from previous calculation
507 MFloat polar, rsquared, var1, var2;
508
509 // If no deviate has been stored, the polar Box-Muller transformation is
510 // performed, producing two independent normally-distributed random
511 // deviates. One is stored for the next round, and one is returned.
512 if(!deviateAvailable) {
513 // choose pairs of uniformly distributed deviates, discarding those
514 // that don't fall within the unit circle
515 do {
516 var1 = 2.0 * (MFloat(rand()) / MFloat(RAND_MAX)) - 1.0;
517 var2 = 2.0 * (MFloat(rand()) / MFloat(RAND_MAX)) - 1.0;
518 rsquared = var1 * var1 + var2 * var2;
519 } while(rsquared >= 1.0 || approx(rsquared, 0.0, MFloatEps));
520
521 // calculate polar tranformation for each deviate
522 polar = sqrt(-2.0 * log(rsquared) / rsquared);
523
524 // store first deviate and set flag
525 storedDeviate = var1 * polar;
526 deviateAvailable = true;
527
528 // return second deviate
529
530 return var2 * polar * sigma + mu;
531 }
532
533 // If a deviate is available from a previous call to this function, it is
534 // returned, and the flag is set to false.
535 else {
536 deviateAvailable = false;
537
538 return storedDeviate * sigma + mu;
539 }
540}
541
548inline std::complex<MFloat>* getFourierCoefficient(MFloat* k, MFloat k0, const MFloat Ma) {
549 TRACE();
550
551 MFloat r[6], s[6], kAbs, energy;
552 std::complex<MFloat> uHat, vHat, wHat;
553 std::complex<MFloat>* fourierCoefficient;
554
555 kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
556
557 // the zero-frequency component is always set to zero, so there is no offset
558 if(approx(kAbs, 0.0, MFloatEps)) {
559 fourierCoefficient = new std::complex<MFloat>[3];
560
561 fourierCoefficient[0] = std::complex<MFloat>(0, 0);
562 fourierCoefficient[1] = std::complex<MFloat>(0, 0);
563 fourierCoefficient[2] = std::complex<MFloat>(0, 0);
564
565
566 return fourierCoefficient;
567
568 } else {
569 // energy = (kAbs/k0)*(kAbs/k0)*(kAbs/k0)*(kAbs/k0) * exp(-2.0*(kAbs/k0)*(kAbs/k0));
570 // energy = pow(kAbs/k0,8.0) * exp(-4.0*(kAbs/k0)*(kAbs/k0)); // set spectral distribution
571 energy = pow(kAbs / k0, 4.0) * exp(-2.0 * (kAbs / k0) * (kAbs / k0)); // set spectral distribution
572 energy *= exp(2.0) * 0.499 * (Ma * LBCS)
573 * (Ma * LBCS); // set maximal fluctuation amplitude to 20% of the freestream velocity (for 128^3: 0.88)
574
575 // determine Fourier coefficients:
576 // r and s are Independant random vector fields with independant
577 // components (zero mean and rms according to energy spectrum).
578 // Each vector has three components for k and another three for -k.
579
580 for(MInt i = 0; i < 6; i++) {
581 r[i] = randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
582 // r[i] = randNumGen.randNorm(0.0, PI*sqrt(energy)/(SQRT2*kAbs));
583 s[i] = randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
584 // s[i] = randNumGen.randNorm(0.0, PI*sqrt(energy)/(SQRT2*kAbs));
585 }
586
587 uHat = (1.0 - k[0] * k[0] / (kAbs * kAbs)) * std::complex<MFloat>(r[0] + r[3], s[0] - s[3])
588 - k[0] * k[1] / (kAbs * kAbs) * std::complex<MFloat>(r[1] + r[4], s[1] - s[4])
589 - k[0] * k[2] / (kAbs * kAbs) * std::complex<MFloat>(r[2] + r[5], s[2] - s[5]);
590
591
592 vHat = -k[1] * k[0] / (kAbs * kAbs) * std::complex<MFloat>(r[0] + r[3], s[0] - s[3])
593 + (1.0 - k[1] * k[1] / (kAbs * kAbs)) * std::complex<MFloat>(r[1] + r[4], s[1] - s[4])
594 - k[1] * k[2] / (kAbs * kAbs) * std::complex<MFloat>(r[2] + r[5], s[2] - s[5]);
595
596
597 wHat = -k[2] * k[0] / (kAbs * kAbs) * std::complex<MFloat>(r[0] + r[3], s[0] - s[3])
598 - k[2] * k[1] / (kAbs * kAbs) * std::complex<MFloat>(r[1] + r[4], s[1] - s[4])
599 + (1.0 - k[2] * k[2] / (kAbs * kAbs)) * std::complex<MFloat>(r[2] + r[5], s[2] - s[5]);
600
601 fourierCoefficient = new std::complex<MFloat>[3];
602
603 // fourierCoefficient[0] = std::complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);// uHat;
604 // fourierCoefficient[1] = std::complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);//vHat;
605 // fourierCoefficient[2] = std::complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);//wHat;
606
607 fourierCoefficient[0] = uHat;
608 fourierCoefficient[1] = vHat;
609 fourierCoefficient[2] = wHat;
610 return fourierCoefficient;
611 }
612}
613
614
620inline std::complex<MFloat>* getFourierCoefficient2(MFloat* k, MFloat k0) {
621 static std::mt19937 randNumGen;
622 static std::uniform_real_distribution<> distrib{0.0, 1.0};
623
624 MFloat r[3], kAbs, energy;
625 std::complex<MFloat> uHat, vHat, wHat;
626 std::complex<MFloat>* fourierCoefficient;
627
628 kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
629
630 // the zero-frequency component is always set to zero, so there is no offset
631 if(approx(kAbs, 0.0, MFloatEps)) {
632 fourierCoefficient = new std::complex<MFloat>[3];
633
634 fourierCoefficient[0] = std::complex<MFloat>(0, 0);
635 fourierCoefficient[1] = std::complex<MFloat>(0, 0);
636 fourierCoefficient[2] = std::complex<MFloat>(0, 0);
637
638
639 return fourierCoefficient;
640
641 } else {
642 // energy = pow(kAbs/k0,4.0) * exp(-2.0*(kAbs/k0)*(kAbs/k0));
643 // energy = ( m_Ma*LBCS * m_Ma*LBCS ) * pow(kAbs/k0,4.0) * exp(-2.0*(kAbs/k0)*(kAbs/k0));
644 // energy = 0.135335;
645 energy = (kAbs / k0) * (kAbs / k0) * (kAbs / k0) * (kAbs / k0) * exp(-2.0 * (kAbs / k0) * (kAbs / k0));
646
647 // determine Fourier coefficients:
648 for(MInt i = 0; i < 3; i++) {
649 r[i] = 2.0 * PI * distrib(randNumGen);
650 }
651
652 uHat = std::complex<MFloat>(cos(r[0]) * energy, sin(r[0]) * energy);
653 vHat = std::complex<MFloat>(cos(r[1]) * energy, sin(r[1]) * energy);
654 wHat = std::complex<MFloat>(cos(r[2]) * energy, sin(r[2]) * energy);
655
656 fourierCoefficient = new std::complex<MFloat>[3];
657
658 fourierCoefficient[0] = uHat;
659 fourierCoefficient[1] = vHat;
660 fourierCoefficient[2] = wHat;
661
662
663 return fourierCoefficient;
664 }
665}
666
681inline void initFftFilter(MFloat* uPhysFieldCoarse,
682 MFloat* vPhysFieldCoarse,
683 MFloat* wPhysFieldCoarse,
684 MInt lx,
685 MInt ly,
686 MInt lz,
687 MInt lxCoarse,
688 MInt lyCoarse,
689 MInt lzCoarse,
690 MInt noPeakModes,
691 const MFloat Ma) {
692 TRACE();
693
694
695 MFloat waveVector[3], k0;
696
697 std::complex<MFloat>* fourierCoefficient;
698
699 fftw_complex *uHatField, *vHatField, *wHatField;
700 fftw_complex *uPhysField, *vPhysField, *wPhysField;
701
702 fftw_plan planU, planV, planW;
703
704 MInt xRatio = lx / lxCoarse;
705 MInt yRatio = ly / lyCoarse;
706 MInt zRatio = lz / lzCoarse;
707
708 if(lx % 2 != 0 || ly % 2 != 0 || lz % 2 != 0) {
709 mTerm(1, AT_, " FFTInit: domainsize must NOT be an odd number! ");
710 }
711
712 m_log << " --- initializing FFTW --- " << std::endl;
713 m_log << " domain size = " << lx << "x" << ly << "x" << lz << std::endl;
714
715 // allocate field of velocities (is deleted at the end of the method)
716 uPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
717 vPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
718 wPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
719
720 // allocate field of Fourier coefficients (is deleted after the transformation)
721 uHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
722 vHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
723 wHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
724
725 // create plans for FFTW
726 planU = fftw_plan_dft_3d(lx, ly, lz, uHatField, uPhysField, FFTW_BACKWARD, FFTW_MEASURE);
727 planV = fftw_plan_dft_3d(lx, ly, lz, vHatField, vPhysField, FFTW_BACKWARD, FFTW_MEASURE);
728 planW = fftw_plan_dft_3d(lx, ly, lz, wHatField, wPhysField, FFTW_BACKWARD, FFTW_MEASURE);
729
730 // reset coefficients and velocities
731 for(MInt p = 0; p < lx; p++) {
732 for(MInt q = 0; q < ly; q++) {
733 for(MInt r = 0; r < lz; r++) {
734 uHatField[r + lz * (q + ly * p)][0] = 0.0;
735 uHatField[r + lz * (q + ly * p)][1] = 0.0;
736 uPhysField[r + lz * (q + ly * p)][0] = 0.0;
737 uPhysField[r + lz * (q + ly * p)][1] = 0.0;
738
739 vHatField[r + lz * (q + ly * p)][0] = 0.0;
740 vHatField[r + lz * (q + ly * p)][1] = 0.0;
741 vPhysField[r + lz * (q + ly * p)][0] = 0.0;
742 vPhysField[r + lz * (q + ly * p)][1] = 0.0;
743
744 wHatField[r + lz * (q + ly * p)][0] = 0.0;
745 wHatField[r + lz * (q + ly * p)][1] = 0.0;
746 wPhysField[r + lz * (q + ly * p)][0] = 0.0;
747 wPhysField[r + lz * (q + ly * p)][1] = 0.0;
748 }
749 }
750 }
751
752 // Set coefficients:
753 // FFTW stores the coefficients for positive wavenumbers in the first half of the array,
754 // and those for negative wavenumbers in reverse order in the second half.
755 // [0, 1, ... , N/2-1, N/2, ... , N-1]
756 // - the entry for zero-wavenumber is at position 0
757 // - the k-th entry and the (N-k)th entry correspond to wavenumbers with opposite sign
758 // - the entry at position N/2 corresponds to the Nyquist wavenumber and appears only once
759
760 // peak wave number of energy spectrum
761 k0 = 2.0 * PI / (lx / noPeakModes);
762
763 for(MInt p = 0; p <= lx / 2; p++) {
764 for(MInt q = 0; q <= ly / 2; q++) {
765 for(MInt r = 0; r <= lz / 2; r++) {
766 // wave-vector: k(p,q,r) = (2 \pi p / lx, 2 \pi q / ly, 2 \pi r / lz)
767 waveVector[0] = (p)*2.0 * PI / lx;
768 waveVector[1] = (q)*2.0 * PI / ly;
769 waveVector[2] = (r)*2.0 * PI / lz;
770
771 fourierCoefficient = getFourierCoefficient(waveVector, k0, Ma);
772
773 // 1. Positive frequencies:
774 uHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[0]);
775 uHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[0]);
776
777 vHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[1]);
778 vHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[1]);
779
780 wHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[2]);
781 wHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[2]);
782
783 // 2. Negative frequencies:
784 if(p > 1 && q > 1 && r > 1) {
785 if(p < lx / 2 && q < ly / 2 && r < lz / 2) {
786 // since the physical velocity field is real, the coefficients for negative frequencies
787 // are the std::complex conjugate of those for positive frequencies
788 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = uHatField[r + lz * (q + ly * p)][0];
789 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -uHatField[r + lz * (q + ly * p)][1];
790
791 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = vHatField[r + lz * (q + ly * p)][0];
792 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -vHatField[r + lz * (q + ly * p)][1];
793
794 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = wHatField[r + lz * (q + ly * p)][0];
795 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -wHatField[r + lz * (q + ly * p)][1];
796 }
797 }
798 }
799 }
800 }
801
802 // Do Fourier transform (backward, see plan definition)
803 // Definition in one dimension:
804 // u(x) = \sum_{j=0}^{lx-1} \hat{u}_j exp(i 2 \pi j x / lx)
805 fftw_execute(planU);
806 fftw_execute(planV);
807 fftw_execute(planW);
808
809 // normalize (this preserves the norm of the basis functions)
810 for(MInt p = 0; p < lx; p++) {
811 for(MInt q = 0; q < ly; q++) {
812 for(MInt r = 0; r < lz; r++) {
813 uPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
814 vPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
815 wPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
816
817 uPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
818 vPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
819 wPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
820 }
821 }
822 }
823
824 // // test:
825 // // cos(x) in a cube with 32x32x32 cells
826 // uHatField[(0+lz*(0+ly*2)=2048][0]=90.5;
827 // uHatField[(0+lz*(0+ly*30)=30720][0]=90.5;
828
829 // reset coarse fields
830 for(MInt p = 0; p < lxCoarse; p++) {
831 for(MInt q = 0; q < lyCoarse; q++) {
832 for(MInt r = 0; r < lzCoarse; r++) {
833 uPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] = F0;
834 vPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] = F0;
835 wPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] = F0;
836 }
837 }
838 }
839
840 // transfer physField to physFieldCoarse via arithmetic averaging
841 for(MInt p = 0; p < lxCoarse; p++) {
842 for(MInt q = 0; q < lyCoarse; q++) {
843 for(MInt r = 0; r < lzCoarse; r++) {
844 for(MInt i = p * xRatio; i < p * zRatio + xRatio; i++) {
845 for(MInt j = q * yRatio; j < q * zRatio + yRatio; j++) {
846 for(MInt k = r * zRatio; k < r * zRatio + zRatio; k++) {
847 uPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] += uPhysField[k + lz * (j + ly * i)][0];
848 vPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] += vPhysField[k + lz * (j + ly * i)][0];
849 wPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] += wPhysField[k + lz * (j + ly * i)][0];
850 }
851 }
852 }
853 }
854 }
855 }
856 for(MInt p = 0; p < lxCoarse; p++) {
857 for(MInt q = 0; q < lyCoarse; q++) {
858 for(MInt r = 0; r < lzCoarse; r++) {
859 uPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] /= (xRatio * yRatio * zRatio);
860 vPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] /= (xRatio * yRatio * zRatio);
861 wPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] /= (xRatio * yRatio * zRatio);
862 }
863 }
864 }
865
866 fftw_destroy_plan(planU);
867 fftw_destroy_plan(planV);
868 fftw_destroy_plan(planW);
869
870 fftw_free(uHatField);
871 fftw_free(vHatField);
872 fftw_free(wHatField);
873
874 fftw_free(uPhysField);
875 fftw_free(vPhysField);
876 fftw_free(wPhysField);
877}
878
892inline void initFft(fftw_complex* uPhysField,
893 fftw_complex* vPhysField,
894 fftw_complex* wPhysField,
895 MInt lx,
896 MInt ly,
897 MInt lz,
898 MInt noPeakModes,
899 const MFloat Ma) {
900 TRACE();
901
902 MFloat waveVector[3], k0;
903
904 std::complex<MFloat>* fourierCoefficient;
905
906 fftw_complex *uHatField, *vHatField, *wHatField;
907
908 fftw_plan planU, planV, planW;
909
910 if(lx % 2 != 0 || ly % 2 != 0 || lz % 2 != 0) {
911 std::stringstream errorMessage;
912 errorMessage << " FFTInit: domainsize must NOT be an odd number! " << std::endl;
913 mTerm(1, AT_, errorMessage.str());
914 }
915
916 m_log << " --- initializing FFTW --- " << std::endl;
917 m_log << " domain size = " << lx << "x" << ly << "x" << lz << std::endl;
918
919 // allocate field of Fourier coefficients (is deleted after the transformation)
920 uHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
921 vHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
922 wHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
923
924 // create plans for FFTW
925 planU = fftw_plan_dft_3d(lx, ly, lz, uHatField, uPhysField, FFTW_BACKWARD, FFTW_MEASURE);
926 planV = fftw_plan_dft_3d(lx, ly, lz, vHatField, vPhysField, FFTW_BACKWARD, FFTW_MEASURE);
927 planW = fftw_plan_dft_3d(lx, ly, lz, wHatField, wPhysField, FFTW_BACKWARD, FFTW_MEASURE);
928
929 // reset coefficients and velocities
930 for(MInt p = 0; p < lx; p++) {
931 for(MInt q = 0; q < ly; q++) {
932 for(MInt r = 0; r < lz; r++) {
933 uHatField[r + lz * (q + ly * p)][0] = 0.0;
934 uHatField[r + lz * (q + ly * p)][1] = 0.0;
935 uPhysField[r + lz * (q + ly * p)][0] = 0.0;
936 uPhysField[r + lz * (q + ly * p)][1] = 0.0;
937
938 vHatField[r + lz * (q + ly * p)][0] = 0.0;
939 vHatField[r + lz * (q + ly * p)][1] = 0.0;
940 vPhysField[r + lz * (q + ly * p)][0] = 0.0;
941 vPhysField[r + lz * (q + ly * p)][1] = 0.0;
942
943 wHatField[r + lz * (q + ly * p)][0] = 0.0;
944 wHatField[r + lz * (q + ly * p)][1] = 0.0;
945 wPhysField[r + lz * (q + ly * p)][0] = 0.0;
946 wPhysField[r + lz * (q + ly * p)][1] = 0.0;
947 }
948 }
949 }
950
951 // Set coefficients:
952 // FFTW stores the coefficients for positive wavenumbers in the first half of the array,
953 // and those for negative wavenumbers in reverse order in the second half.
954 // [0, 1, ... , N/2-1, N/2, ... , N-1]
955 // - the entry for zero-wavenumber is at position 0
956 // - the k-th entry and the (N-k)th entry correspond to wavenumbers with opposite sign
957 // - the entry at position N/2 corresponds to the Nyquist wavenumber and appears only once
958
959 // peak wave number of energy spectrum
960 k0 = 2.0 * PI / (lx / noPeakModes);
961
962 for(MInt p = 0; p <= lx / 2; p++) {
963 for(MInt q = 0; q <= ly / 2; q++) {
964 for(MInt r = 0; r <= lz / 2; r++) {
965 // wave-vector: k(p,q,r) = (2 \pi p / lx, 2 \pi q / ly, 2 \pi r / lz)
966 waveVector[0] = (p)*2.0 * PI / lx;
967 waveVector[1] = (q)*2.0 * PI / ly;
968 waveVector[2] = (r)*2.0 * PI / lz;
969
970 fourierCoefficient = getFourierCoefficient(waveVector, k0, Ma);
971
972 // 1. Positive frequencies:
973 uHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[0]);
974 uHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[0]);
975
976 vHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[1]);
977 vHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[1]);
978
979 wHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[2]);
980 wHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[2]);
981
982 // 2. Negative frequencies:
983 if(p > 1 && q > 1 && r > 1) {
984 if(p < lx / 2 && q < ly / 2 && r < lz / 2) {
985 // since the physical velocity field is real, the coefficients for negative frequencies
986 // are the std::complex conjugate of those for positive frequencies
987 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = uHatField[r + lz * (q + ly * p)][0];
988 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -uHatField[r + lz * (q + ly * p)][1];
989
990 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = vHatField[r + lz * (q + ly * p)][0];
991 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -vHatField[r + lz * (q + ly * p)][1];
992
993 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = wHatField[r + lz * (q + ly * p)][0];
994 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -wHatField[r + lz * (q + ly * p)][1];
995 }
996 }
997 }
998 }
999 }
1000
1001 // // test:
1002 // // cos(x) in a cube with 32x32x32 cells
1003 // uHatField[(0+lz*(0+ly*2)=2048][0]=90.5;
1004 // uHatField[(0+lz*(0+ly*30)=30720][0]=90.5;
1005
1006 // Do Fourier transform (backward, see plan definition)
1007 // Definition in one dimension:
1008 // u(x) = \sum_{j=0}^{lx-1} \hat{u}_j exp(i 2 \pi j x / lx)
1009 fftw_execute(planU);
1010 fftw_execute(planV);
1011 fftw_execute(planW);
1012
1013 // normalize (this preserves the norm of the basis functions)
1014 for(MInt p = 0; p < lx; p++) {
1015 for(MInt q = 0; q < ly; q++) {
1016 for(MInt r = 0; r < lz; r++) {
1017 uPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
1018 vPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
1019 wPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
1020
1021 uPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
1022 vPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
1023 wPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
1024 }
1025 }
1026 }
1027
1028 fftw_destroy_plan(planU);
1029 fftw_destroy_plan(planV);
1030 fftw_destroy_plan(planW);
1031 fftw_free(uHatField);
1032 fftw_free(vHatField);
1033 fftw_free(wHatField);
1034}
1035
1036
1051inline void computeEnergySpectrum(MFloatScratchSpace& velocity, MIntScratchSpace& indices, const ptrdiff_t howmany,
1052 MInt locSize, MInt nx, MInt ny, MInt nz, MFloat viscosity, const MPI_Comm comm) {
1053 DEBUG("computeEnergySpectrum entry", MAIA_DEBUG_TRACE_IN);
1054 TRACE();
1055
1056 const MInt size = nx * ny * nz;
1057 const MFloat fsize = (MFloat)size;
1058 const ptrdiff_t rank = 3;
1059 if(nx % 2 != 0 || ny % 2 != 0 || nz % 2 != 0) {
1060 std::stringstream errorMessage;
1061 errorMessage << " FFTInit: domainsize must NOT be an odd number! " << nx << " " << ny << " " << nz << std::endl;
1062 mTerm(1, AT_, errorMessage.str());
1063 }
1064 if(howmany < 3) mTerm(1, AT_, "expecting at least 3D velocity field");
1065
1066 m_log << " --- initializing FFTW --- " << std::endl;
1067 m_log << " domain size = " << nx << "x" << ny << "x" << nz << " (" << howmany << ")" << std::endl;
1068
1069 // caution: nz must be evenly divisble by number of ranks!
1070 MInt domainId;
1071 MInt noDomains;
1072 MPI_Comm_rank(comm, &domainId);
1073 MPI_Comm_size(comm, &noDomains);
1074 MPI_Comm MPI_COMM_FFTW;
1075 MInt maxRank = 1;
1076 maxRank = mMin(nz, noDomains);
1077 while(nz % maxRank != 0) {
1078 maxRank--;
1079 }
1080 MInt color = (domainId < maxRank) ? 0 : MPI_UNDEFINED;
1081 MPI_Comm_split(comm, color, domainId, &MPI_COMM_FFTW, AT_, "MPI_COMM_FFTW");
1082 m_log << "no ranks for fft: " << maxRank << std::endl;
1083
1084 MFloat urms0 = F0;
1085 for(MInt i = 0; i < locSize; i++) {
1086 urms0 += POW2(velocity(i, 0)) + POW2(velocity(i, 1)) + POW2(velocity(i, 2));
1087 }
1088 MPI_Allreduce(MPI_IN_PLACE, &urms0, 1, MPI_DOUBLE, MPI_SUM, comm, AT_, "MPI_IN_PLACE", "urms0");
1089 urms0 = sqrt(F1B3 * urms0 / fsize);
1090
1091 MInt fftLocSize = 0;
1092 ptrdiff_t alloc_local = 0, local_n0 = 0, local_0_start = 0;
1093 fftw_complex* hatField = nullptr;
1094 fftw_plan plan;
1095
1096 const MFloat time0 = MPI_Wtime();
1097
1098 const ptrdiff_t n[3] = {nx, ny, nz};
1099#ifndef MAIA_WINDOWS
1100 alloc_local = (domainId < maxRank) ? fftw_mpi_local_size_many(rank, n, howmany, FFTW_MPI_DEFAULT_BLOCK, MPI_COMM_FFTW,
1101 &local_n0, &local_0_start)
1102 : 1;
1103#endif
1104
1105 ScratchSpace<fftw_complex> hatFieldMem(alloc_local, AT_, "hatFieldMem");
1106 if(domainId < maxRank) {
1107 hatField = &hatFieldMem[0];
1108 } else {
1109 alloc_local = 0;
1110 local_n0 = 0;
1111 local_0_start = nx;
1112 }
1113
1114 ScratchSpace<MLong> offsets(noDomains + 1, AT_, "offsets");
1115 ScratchSpace<MInt> noRecv(noDomains, AT_, "noRecv");
1116 ScratchSpace<MInt> noSend(noDomains, AT_, "noSend");
1117 noRecv.fill(0);
1118 noSend.fill(0);
1119 MLong locOffset = (domainId < maxRank) ? ((MInt)local_0_start) * ny * nz : size;
1120
1121 MPI_Allgather(&locOffset, 1, MPI_LONG, &offsets[0], 1, MPI_LONG, comm, AT_, "locOffset", "offsets[0]");
1122 offsets(noDomains) = size;
1123 if(domainId < maxRank) {
1124 if(offsets(domainId) != locOffset || offsets(domainId + 1) != ((MInt)(local_0_start + local_n0)) * ny * nz)
1125 mTerm(1, AT_, "wrong size 0");
1126 }
1127 for(MInt i = 0; i < locSize; i++) {
1128 MInt j = indices(i);
1129 MInt nghbrDomain = mMin(maxRank - 1, j / (size / maxRank));
1130 while(j < offsets(nghbrDomain) || j >= offsets(nghbrDomain + 1)) {
1131 if(j < offsets(nghbrDomain)) nghbrDomain--;
1132 if(j >= offsets(nghbrDomain + 1)) nghbrDomain++;
1133 }
1134 if(nghbrDomain >= maxRank) mTerm(1, AT_, "wrong domain");
1135 noSend(nghbrDomain)++;
1136 }
1137 MPI_Alltoall(&noSend[0], 1, MPI_INT, &noRecv[0], 1, MPI_INT, comm, AT_, "noSend[0]", "noRecv[0]");
1138 MInt sendSize = 0;
1139 MInt recvSize = 0;
1140 ScratchSpace<MInt> recvOffsets(noDomains + 1, AT_, "recvOffsets");
1141 ScratchSpace<MInt> sendOffsets(noDomains + 1, AT_, "sendOffsets");
1142 for(MInt i = 0; i < noDomains; i++) {
1143 recvOffsets(i) = recvSize;
1144 sendOffsets(i) = sendSize;
1145 sendSize += noSend(i);
1146 recvSize += noRecv(i);
1147 }
1148 recvOffsets(noDomains) = recvSize;
1149 sendOffsets(noDomains) = sendSize;
1150 ScratchSpace<MFloat> recvData(mMax(1, recvSize), howmany, AT_, "recvData");
1151 ScratchSpace<MFloat> sendData(mMax(1, sendSize), howmany, AT_, "sendData");
1152 ScratchSpace<MInt> recvIndices(mMax(1, recvSize), AT_, "recvIndices");
1153 ScratchSpace<MInt> sendIndices(mMax(1, sendSize), AT_, "sendIndices");
1154 noSend.fill(0);
1155 for(MInt i = 0; i < locSize; i++) {
1156 MInt j = indices(i);
1157 MInt nghbrDomain = mMin(maxRank - 1, j / (size / maxRank));
1158 while(j < offsets(nghbrDomain) || j >= offsets(nghbrDomain + 1)) {
1159 if(j < offsets(nghbrDomain)) nghbrDomain--;
1160 if(j >= offsets(nghbrDomain + 1)) nghbrDomain++;
1161 }
1162 if(nghbrDomain >= maxRank) mTerm(1, AT_, "wrong domain");
1163 if(sendOffsets(nghbrDomain) + noSend(nghbrDomain) >= sendOffsets(nghbrDomain + 1)) mTerm(1, AT_, "wrong size");
1164 for(MInt k = 0; k < howmany; k++) {
1165 sendData(sendOffsets(nghbrDomain) + noSend(nghbrDomain), k) = velocity(i, k);
1166 }
1167 sendIndices(sendOffsets(nghbrDomain) + noSend(nghbrDomain)) = j;
1168 if(j < offsets(nghbrDomain) || j >= offsets(nghbrDomain + 1)) mTerm(1, AT_, "wrong size 01");
1169 noSend(nghbrDomain)++;
1170 }
1171 for(MInt i = 0; i < noDomains; i++) {
1172 if(noSend(i) != sendOffsets(i + 1) - sendOffsets(i)) mTerm(1, AT_, "wrong size 1");
1173 if(noRecv(i) != recvOffsets(i + 1) - recvOffsets(i)) mTerm(1, AT_, "wrong size 1");
1174 }
1175
1176 ScratchSpace<MPI_Request> sendReq(noDomains, AT_, "sendReq");
1177 sendReq.fill(MPI_REQUEST_NULL);
1178 MInt cnt = 0;
1179 for(MInt i = 0; i < noDomains; i++) {
1180 if(noSend(i) == 0) continue;
1181 MPI_Issend(&sendData[howmany * sendOffsets(i)], howmany * noSend(i), MPI_DOUBLE, i, 123, comm, &sendReq[cnt], AT_,
1182 "sendData[howmany * sendOffsets(i)]");
1183 cnt++;
1184 }
1185 for(MInt i = 0; i < noDomains; i++) {
1186 if(noRecv(i) == 0) continue;
1187 MPI_Recv(&recvData[howmany * recvOffsets(i)], howmany * noRecv(i), MPI_DOUBLE, i, 123, comm, MPI_STATUS_IGNORE, AT_,
1188 "recvData[howmany * recvOffsets(i)]");
1189 }
1190 if(cnt > 0) MPI_Waitall(cnt, &sendReq[0], MPI_STATUSES_IGNORE, AT_);
1191 sendReq.fill(MPI_REQUEST_NULL);
1192 cnt = 0;
1193 for(MInt i = 0; i < noDomains; i++) {
1194 if(noSend(i) == 0) continue;
1195 MPI_Issend(&sendIndices[sendOffsets(i)], noSend(i), MPI_INT, i, 124, comm, &sendReq[cnt], AT_,
1196 "sendIndices[sendOffsets(i)]");
1197 cnt++;
1198 }
1199 for(MInt i = 0; i < noDomains; i++) {
1200 if(noRecv(i) == 0) continue;
1201 MPI_Recv(&recvIndices[recvOffsets(i)], noRecv(i), MPI_INT, i, 124, comm, MPI_STATUS_IGNORE, AT_,
1202 "recvIndices[recvOffsets(i)]");
1203 }
1204 if(cnt > 0) MPI_Waitall(cnt, &sendReq[0], MPI_STATUSES_IGNORE, AT_);
1205
1206 if(domainId < maxRank) {
1207 for(MInt i = 0; i < recvSize; i++) {
1208 MInt j = recvIndices(i);
1209 if(j < ((MInt)local_0_start) * ny * nz || j >= ((MInt)(local_0_start + local_n0)) * ny * nz)
1210 mTerm(1, AT_, "wrong size 2");
1211 MInt pos = j - (((MInt)local_0_start) * ny * nz);
1212 for(MInt k = 0; k < howmany; k++) {
1213 hatField[howmany * pos + k][0] = recvData(i, k); // Real part of fourier transformed velocity
1214 hatField[howmany * pos + k][1] = F0; // Imaginary part of fourier transformed velocity
1215 }
1216 }
1217 fftLocSize = recvSize;
1218
1219 MFloat couplcheck = F0;
1220 for(MInt i = 0; i < fftLocSize; i++) {
1221 for(MInt k = 0; k < 3; k++) {
1222 couplcheck -= hatField[howmany * i + k][0] * hatField[howmany * i + 3 + k][0];
1223 }
1224 }
1225 MPI_Allreduce(MPI_IN_PLACE, &couplcheck, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE", "couplcheck");
1226 if(domainId == 0)
1227 std::cerr << "couplecheck fft " << couplcheck << " " << couplcheck * POW3(0.25 / 64.0) << std::endl;
1228 }
1229
1230 // Perform Fourier transform
1231 if(domainId < maxRank) {
1232#ifndef MAIA_WINDOWS
1233 plan = fftw_mpi_plan_many_dft(rank, n, howmany, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, hatField, hatField,
1234 MPI_COMM_FFTW, FFTW_FORWARD, FFTW_ESTIMATE);
1235 fftw_execute(plan);
1236 fftw_destroy_plan(plan);
1237#endif
1238 for(MInt i = 0; i < fftLocSize; i++) {
1239 for(MInt j = 0; j < howmany; j++) {
1240 hatField[howmany * i + j][0] /= fsize;
1241 hatField[howmany * i + j][1] /= fsize;
1242 }
1243 }
1244 }
1245
1246 const MFloat time1 = MPI_Wtime();
1247
1248 m_log << "parallel fftw time " << time1 - time0 << std::endl;
1249
1250 // Set coefficients:
1251 // FFTW stores the coefficients for positive wavenumbers in the first half of the array,
1252 // and those for negative wavenumbers in reverse order in the second half.
1253 // [0, 1, ... , N/2-1, N/2, ... , N-1]
1254 // - the entry for zero-wavenumber is at position 0
1255 // - the k-th entry and the (N-k)th entry correspond to wavenumbers with opposite sign
1256 // - the entry at position N/2 corresponds to the Nyquist wavenumber and appears only once
1257
1258 static const MInt referenceCubeSize =
1259 ((Context::propertyExists("referenceCubeSize")) ? Context::getBasicProperty<MInt>("referenceCubeSize", AT_) : 1);
1260
1261 static const MInt kMinSpec = ((Context::propertyExists("kMin")) ? Context::getBasicProperty<MInt>("kMin", AT_) : 0);
1262
1263 static const MInt kMaxSpec = ((Context::propertyExists("kMax")) ? Context::getBasicProperty<MInt>("kMax", AT_) : 0);
1264
1265 const MFloat k0 = F2 * PI / referenceCubeSize;
1266 MInt noBins = 0;
1267 MFloat kMin = 0.0;
1268 MFloat kMax = 0.0;
1269 MFloat deltaK = 0.0;
1270 if(kMinSpec > 0 && kMaxSpec > 0) {
1271 noBins = kMaxSpec;
1272 kMin = static_cast<float>(kMinSpec) * k0;
1273 kMax = static_cast<float>(kMaxSpec) * k0;
1274 deltaK = (kMax - kMin) / ((MFloat)(kMaxSpec - 1));
1275 } else {
1276 noBins = mMax(nx, mMax(ny, nz)) / 2;
1277 kMin = F1 * k0;
1278 kMax = noBins * k0;
1279 deltaK = (kMax - kMin) / ((MFloat)(noBins - 1));
1280 }
1281
1282 // const MInt noBins = 200;
1283 // const MFloat kMax = 200.0*k0;
1284 ScratchSpace<MFloat> spectrumBnd(noBins + 1, AT_, "spectrumBnd");
1285 ScratchSpace<MFloat> spectrum(noBins, AT_, "spectrum");
1286 ScratchSpace<MFloat> coupling(noBins, AT_, "coupling");
1287 ScratchSpace<MFloat> transfer(noBins, AT_, "transfer");
1288 ScratchSpace<MFloat> transfer2(noBins, AT_, "transfer2");
1289 ScratchSpace<MFloat> rate(noBins, AT_, "rate");
1290 ScratchSpace<MFloat> flux(noBins, 2, AT_, "flux");
1291 ScratchSpace<MFloat> flux2(noBins, 2, AT_, "flux2");
1292 ScratchSpace<MFloat> spectrumSum(noBins, AT_, "spectrumSum");
1293 ScratchSpace<MFloat> spectrumSum2(noBins, AT_, "spectrumSum2");
1294 spectrumBnd[0] = kMin - F1B2 * deltaK;
1295 const MFloat spectrumBnd00 = kMin - F3B4 * deltaK;
1296 for(MInt i = 0; i < noBins; i++) {
1297 spectrumBnd[i + 1] = spectrumBnd[i] + deltaK;
1298 }
1299 spectrum.fill(F0);
1300 coupling.fill(F0);
1301 transfer.fill(F0);
1302 transfer2.fill(F0);
1303 rate.fill(F0);
1304 flux.fill(F0);
1305 flux2.fill(F0);
1306 spectrumSum.fill(F0);
1307 spectrumSum2.fill(F0);
1308
1309 // beware: very expensive routine
1310 static const MBool computeTransferRate =
1311 ((Context::propertyExists("computeTransferRate")) ? Context::getBasicProperty<MBool>("computeTransferRate", AT_)
1312 : false);
1313 if(computeTransferRate) {
1314 MIntScratchSpace noDatPerDomain(noDomains, AT_, "noDatPerDomain");
1315 MIntScratchSpace dataOffsets(noDomains + 1, AT_, "dataOffsets");
1316 MPI_Allgather(&fftLocSize, 1, MPI_INT, &noDatPerDomain[0], 1, MPI_INT, comm, AT_, "fftLocSize",
1317 "noDatPerDomain[0]");
1318 dataOffsets.fill(0);
1319 for(MInt d = 0; d < noDomains; d++) {
1320 dataOffsets(d + 1) = dataOffsets(d) + noDatPerDomain(d);
1321 }
1322
1323 MPI_Request sendReqLoc = MPI_REQUEST_NULL;
1324 fftw_complex* uHatGlobal = fftw_alloc_complex(3 * size);
1325 fftw_complex* uHatGlobalTmp = nullptr;
1326 for(MInt i = 0; i < fftLocSize; i++) {
1327 for(MInt j = 0; j < 3; j++) {
1328 uHatGlobal[3 * i + j][0] = hatField[howmany * i + j][0];
1329 uHatGlobal[3 * i + j][1] = hatField[howmany * i + j][1];
1330 }
1331 }
1332 if(fftLocSize > 0) {
1333 MPI_Issend(&uHatGlobal[0], 6 * fftLocSize, MPI_DOUBLE, 0, 18, comm, &sendReqLoc, AT_, "uHatGlobal[0]");
1334 }
1335 if(domainId == 0) {
1336 uHatGlobalTmp = fftw_alloc_complex(3 * size);
1337 for(MInt d = 0; d < noDomains; d++) {
1338 if(noDatPerDomain[d] > 0) {
1339 MPI_Recv(&(uHatGlobalTmp[3 * dataOffsets[d]][0]), 6 * noDatPerDomain[d], MPI_DOUBLE, d, 18, comm,
1340 MPI_STATUSES_IGNORE, AT_, "(uHatGlobalTmp[3 * dataOffsets[d]][0])");
1341 }
1342 }
1343 }
1344 if(fftLocSize > 0) {
1345 MPI_Wait(&sendReqLoc, MPI_STATUSES_IGNORE, AT_);
1346 }
1347
1348 if(domainId == 0) {
1349 // re-sort by ascending wave numbers
1350 for(MInt i0 = 0; i0 < nx; i0++) {
1351 for(MInt i1 = 0; i1 < ny; i1++) {
1352 for(MInt i2 = 0; i2 < nz; i2++) {
1353 MInt j0 = (i0 > nx / 2) ? i0 - nx : i0;
1354 MInt j1 = (i1 > ny / 2) ? i1 - ny : i1;
1355 MInt j2 = (i2 > nz / 2) ? i2 - nz : i2;
1356 j0 += nx / 2 - 1;
1357 j1 += nx / 2 - 1;
1358 j2 += nx / 2 - 1;
1359 MInt pos0 = getGlobalPosFFTW(i0, i1, i2, ny, nz);
1360 MInt pos1 = getGlobalPosFFTW(j0, j1, j2, ny, nz);
1361 for(MInt j = 0; j < 3; j++) {
1362 uHatGlobal[3 * pos1 + j][0] = uHatGlobalTmp[3 * pos0 + j][0];
1363 uHatGlobal[3 * pos1 + j][1] = uHatGlobalTmp[3 * pos0 + j][1];
1364 }
1365 }
1366 }
1367 }
1368 fftw_free(uHatGlobalTmp);
1369 }
1370
1371 MPI_Bcast(&uHatGlobal[0], 6 * size, MPI_DOUBLE, 0, comm, AT_, "uHatGlobal");
1372
1373
1374 // MInt domainSize = size/noDomains;
1375 // MInt posMin = domainSize*domainId;
1376 // MInt posMax = mMin(size,domainSize*(domainId+1));
1377 // cerr << domainId << ": pos " << posMin << " " << posMax << std::endl;
1378
1379 MInt noWaveNumbers = 0;
1380 MLong weightSum = 0;
1381 for(MInt pos = 0; pos < size; pos++) {
1382 MInt i2 = pos % nz;
1383 MInt i1 = (pos / nz) % ny;
1384 MInt i0 = pos / (ny * nz);
1385 if(pos != i2 + nz * (i1 + ny * (i0))) mTerm(1, AT_, "index mismatch");
1386 MFloat k[3];
1387 MInt j0 = i0 - nx / 2 + 1;
1388 MInt j1 = i1 - ny / 2 + 1;
1389 MInt j2 = i2 - nz / 2 + 1;
1390 k[0] = ((MFloat)j0) * k0;
1391 k[1] = ((MFloat)j1) * k0;
1392 k[2] = ((MFloat)j2) * k0;
1393
1394 MFloat kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
1395 if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
1396 if(kAbs < spectrumBnd[0] - 1e-12) continue;
1397
1398 MInt p0Min = std::max(0, i0 - nx / 2);
1399 MInt p0Max = std::min(nx, i0 + nx / 2);
1400 MInt p1Min = std::max(0, i1 - ny / 2);
1401 MInt p1Max = std::min(ny, i1 + ny / 2);
1402 MInt p2Min = std::max(0, i2 - nz / 2);
1403 MInt p2Max = std::min(nz, i2 + nz / 2);
1404 MLong weight = ((p0Max - p0Min) * (p1Max - p1Min) * (p2Max - p2Min));
1405
1406 weightSum += weight;
1407 noWaveNumbers++;
1408 }
1409 // MInt domainSize = noWaveNumbers/noDomains;
1410 MLong noDomainsLong = (MLong)noDomains;
1411 MLong domainIdLong = (MLong)domainId;
1412 MLong domainSize = weightSum / noDomainsLong;
1413 MInt posMin = 0;
1414 MInt posMax = size;
1415 noWaveNumbers = 0;
1416 weightSum = 0;
1417 for(MInt pos = 0; pos < size; pos++) {
1418 MInt i2 = pos % nz;
1419 MInt i1 = (pos / nz) % ny;
1420 MInt i0 = pos / (ny * nz);
1421 if(pos != i2 + nz * (i1 + ny * (i0))) mTerm(1, AT_, "index mismatch");
1422 MFloat k[3];
1423 MInt j0 = i0 - nx / 2 + 1;
1424 MInt j1 = i1 - ny / 2 + 1;
1425 MInt j2 = i2 - nz / 2 + 1;
1426 k[0] = ((MFloat)j0) * k0;
1427 k[1] = ((MFloat)j1) * k0;
1428 k[2] = ((MFloat)j2) * k0;
1429
1430 MFloat kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
1431 if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
1432 if(kAbs < spectrumBnd[0] - 1e-12) continue;
1433
1434 MInt p0Min = std::max(0, i0 - nx / 2);
1435 MInt p0Max = std::min(nx, i0 + nx / 2);
1436 MInt p1Min = std::max(0, i1 - ny / 2);
1437 MInt p1Max = std::min(ny, i1 + ny / 2);
1438 MInt p2Min = std::max(0, i2 - nz / 2);
1439 MInt p2Max = std::min(nz, i2 + nz / 2);
1440 MLong weight = ((p0Max - p0Min) * (p1Max - p1Min) * (p2Max - p2Min));
1441
1442 if(weightSum <= (domainSize * domainIdLong) && (weightSum + weight) > (domainSize * domainIdLong)) posMin = pos;
1443 if(weightSum <= (domainSize * (domainIdLong + 1)) && (weightSum + weight) > (domainSize * (domainIdLong + 1)))
1444 posMax = pos;
1445
1446 weightSum += weight;
1447 noWaveNumbers++;
1448 }
1449
1450 weightSum = 0;
1451 for(MInt pos = posMin; pos < posMax; pos++) {
1452 const MInt i2 = pos % nz;
1453 const MInt i1 = (pos / nz) % ny;
1454 const MInt i0 = pos / (ny * nz);
1455
1456 const MInt j0 = i0 - nx / 2 + 1;
1457 const MInt j1 = i1 - ny / 2 + 1;
1458 const MInt j2 = i2 - nz / 2 + 1;
1459 const MFloat k[3] = {((MFloat)j0) * k0, ((MFloat)j1) * k0, ((MFloat)j2) * k0};
1460 const MFloat kAbs2 = mMax(1e-12, k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
1461 const MFloat kAbs = mMax(1e-12, sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]));
1462
1463 if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
1464 const MInt bin = ((MInt)((kAbs - spectrumBnd[0] + 1e-12) / deltaK));
1465 const MInt bin2 = ((MInt)((kAbs - spectrumBnd00 + 1e-12) / deltaK));
1466
1467 const MFloat kjl[3][3] = {
1468 {1.0 - (k[0] * k[0] / kAbs2), 0.0 - (k[0] * k[1] / kAbs2), 0.0 - (k[0] * k[2] / kAbs2)},
1469 {0.0 - (k[1] * k[0] / kAbs2), 1.0 - (k[1] * k[1] / kAbs2), 0.0 - (k[1] * k[2] / kAbs2)},
1470 {0.0 - (k[2] * k[0] / kAbs2), 0.0 - (k[2] * k[1] / kAbs2), 1.0 - (k[2] * k[2] / kAbs2)}};
1471
1472 const MInt p0Min = std::max(0, i0 - nx / 2);
1473 const MInt p0Max = std::min(nx, i0 + nx / 2);
1474 const MInt p1Min = std::max(0, i1 - ny / 2);
1475 const MInt p1Max = std::min(ny, i1 + ny / 2);
1476 const MInt p2Min = std::max(0, i2 - nz / 2);
1477 const MInt p2Max = std::min(nz, i2 + nz / 2);
1478
1479 MFloat trans = F0;
1480
1481 for(MInt p0 = p0Min; p0 < p0Max; p0++) {
1482 for(MInt p1 = p1Min; p1 < p1Max; p1++) {
1483 for(MInt p2 = p2Min; p2 < p2Max; p2++) {
1484 const MInt pos1 = p2 + nz * (p1 + ny * (p0));
1485 const MInt s0 = i0 - p0 + nx / 2 - 1;
1486 const MInt s1 = i1 - p1 + ny / 2 - 1;
1487 const MInt s2 = i2 - p2 + nz / 2 - 1;
1488 const MInt pos2 = s2 + nz * (s1 + ny * (s0));
1489
1490 for(MInt i = 0; i < 3; i++) {
1491 for(MInt j = 0; j < 3; j++) {
1492 for(MInt l = 0; l < 3; l++) {
1493 trans += k[i] * kjl[j][l]
1494 * triadImag(uHatGlobal[3 * pos1 + l], uHatGlobal[3 * pos2 + i], uHatGlobal[3 * pos + j]);
1495 }
1496 }
1497 }
1498 }
1499 }
1500 }
1501
1502 if(bin > -1 && bin < noBins) transfer(bin) += trans;
1503 if(bin2 > -1 && bin2 < noBins) transfer2(bin2) += trans;
1504 }
1505
1506 MPI_Allreduce(MPI_IN_PLACE, &transfer[0], noBins, MPI_DOUBLE, MPI_SUM, comm, AT_, "MPI_IN_PLACE", "transfer[0]");
1507 MPI_Allreduce(MPI_IN_PLACE, &transfer2[0], noBins, MPI_DOUBLE, MPI_SUM, comm, AT_, "MPI_IN_PLACE", "transfer2[0]");
1508 fftw_free(uHatGlobal);
1509 }
1510
1511 const MFloat time2 = MPI_Wtime();
1512
1513 if(domainId < maxRank) {
1514 for(MInt i0 = (MInt)local_0_start; i0 < (MInt)(local_0_start + local_n0); i0++) {
1515 for(MInt i1 = 0; i1 < ny; i1++) {
1516 for(MInt i2 = 0; i2 < nz; i2++) {
1517 MFloat k[3];
1518 MInt j0 = (i0 > nx / 2) ? i0 - nx : i0;
1519 MInt j1 = (i1 > ny / 2) ? i1 - ny : i1;
1520 MInt j2 = (i2 > nz / 2) ? i2 - nz : i2;
1521 k[0] = ((MFloat)j0) * k0;
1522 k[1] = ((MFloat)j1) * k0;
1523 k[2] = ((MFloat)j2) * k0;
1524
1525 MInt pos = i2 + nz * (i1 + ny * (i0 - ((MInt)local_0_start)));
1526 if(pos < 0 || pos >= fftLocSize) mTerm(1, AT_, "wrong size 3");
1527
1528 MFloat kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
1529 MFloat eng = F0;
1530 MFloat coupl = F0;
1531 MFloat dedt = F0;
1532 for(MInt q = 0; q < 3; q++) {
1533 eng += F1B2 * (POW2(hatField[howmany * pos + q][0]) + POW2(hatField[howmany * pos + q][1]));
1534 }
1535
1536 // if not only velocity field is passed as argument, also coupling rate Psi(k) and change of Energy dE/dt is
1537 // computed
1538 if(howmany > 3) {
1539 for(MInt q = 0; q < 3; q++) {
1540 coupl += (hatField[howmany * pos + q][0] * hatField[howmany * pos + 3 + q][0]
1541 + hatField[howmany * pos + q][1] * hatField[howmany * pos + 3 + q][1]);
1542 dedt += (hatField[howmany * pos + 6 + q][0] * hatField[howmany * pos + q][0]
1543 + hatField[howmany * pos + 6 + q][1] * hatField[howmany * pos + q][1]);
1544 }
1545 }
1546 if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
1547
1548 MInt bin = ((MInt)((kAbs - spectrumBnd[0] + 1e-12) / deltaK));
1549 MInt bin2 = ((MInt)((kAbs - spectrumBnd00 + 1e-12) / deltaK));
1550 if(bin > -1 && bin < noBins) {
1551 spectrum(bin) += eng;
1552 coupling(bin) += coupl;
1553 rate(bin) += dedt;
1554 spectrumSum(bin) += F1;
1555 }
1556 if(bin2 > -1 && bin2 < noBins) {
1557 spectrumSum2(bin2) += F1;
1558 }
1559 }
1560 }
1561 }
1562
1563 MPI_Allreduce(MPI_IN_PLACE, &spectrum[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
1564 "spectrum[0]");
1565 MPI_Allreduce(MPI_IN_PLACE, &coupling[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
1566 "coupling[0]");
1567 MPI_Allreduce(MPI_IN_PLACE, &rate[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE", "rate[0]");
1568 MPI_Allreduce(MPI_IN_PLACE, &spectrumSum[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
1569 "spectrumSum[0]");
1570 MPI_Allreduce(MPI_IN_PLACE, &spectrumSum2[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
1571 "spectrumSum2[0]");
1572
1573 for(MInt i = 0; i < noBins; i++) {
1574 MFloat db = spectrumBnd00 - spectrumBnd[0];
1575 MFloat vol =
1576 F4B3 * PI * (POW3(spectrumBnd[i + 1]) - POW3(spectrumBnd[i])) / (spectrumBnd[i + 1] - spectrumBnd[i]);
1577 MFloat vol2 = F4B3 * PI * (POW3(spectrumBnd[i + 1] + db) - POW3(spectrumBnd[i] + db))
1578 / (spectrumBnd[i + 1] - spectrumBnd[i]);
1579 spectrum(i) *= vol / (spectrumSum(i) * POW3(k0));
1580 coupling(i) *= vol / (spectrumSum(i) * POW3(k0));
1581 transfer(i) *= vol / (spectrumSum(i) * POW3(k0));
1582 transfer2(i) *= vol2 / (spectrumSum2(i) * POW3(k0));
1583 rate(i) *= vol / (spectrumSum(i) * POW3(k0));
1584 }
1585
1586 for(MInt i = 0; i < noBins; i++) {
1587 flux(i, 0) = F0;
1588 flux(i, 1) = F0;
1589 for(MInt j = 0; j < i; j++) {
1590 flux(i, 0) += transfer(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
1591 }
1592 for(MInt j = i; j < noBins; j++) {
1593 flux(i, 1) += transfer(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
1594 }
1595 flux2(i, 0) = F0;
1596 flux2(i, 1) = F0;
1597 for(MInt j = 0; j < i; j++) {
1598 flux2(i, 0) += transfer2(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
1599 }
1600 for(MInt j = i; j < noBins; j++) {
1601 flux2(i, 1) += transfer2(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
1602 }
1603 }
1604
1605 if(domainId == 0) {
1606 std::ofstream ofl;
1607 std::ofstream ofl2;
1608 ofl.open("./out/energySpectrum_00" + std::to_string(globalTimeStep), std::ios_base::out | std::ios_base::trunc);
1609 ofl2.open("./out/energySpectrum_coarse_00" + std::to_string(globalTimeStep),
1610 std::ios_base::out | std::ios_base::trunc);
1611 if(ofl.is_open() && ofl.good() && ofl2.is_open() && ofl2.good()) {
1612 MFloat energy = F0;
1613 MFloat dissip = F0;
1614 MFloat length = F0;
1615 MFloat dedt = F0;
1616 MFloat force = F0;
1617 MFloat urms = F0;
1618 MFloat transf = F0;
1619 MFloat transf2 = F0;
1620 for(MInt i = 0; i < noBins; i++) {
1621 urms += spectrum(i) * (spectrumBnd[i + 1] - spectrumBnd[i]);
1622 }
1623 urms = sqrt(F2B3 * urms); // warum hier 2/3 und nicht 1/3
1624 for(MInt i = 0; i < noBins; i++) {
1625 MFloat k = F1B2 * (spectrumBnd[i + 1] + spectrumBnd[i]);
1626 MFloat Ek = spectrum(i);
1627 if(std::isnan(Ek) || std::isnan(coupling(i)) || std::isnan(rate(i))) continue;
1628 energy += (spectrumBnd[i + 1] - spectrumBnd[i]) * Ek;
1629 dissip += F2 * viscosity * POW2(k) * (spectrumBnd[i + 1] - spectrumBnd[i]) * Ek;
1630 length += (F1B2 * PI / POW2(urms)) * (spectrumBnd[i + 1] - spectrumBnd[i]) * Ek / (k);
1631 force += (spectrumBnd[i + 1] - spectrumBnd[i]) * coupling(i);
1632 dedt += (spectrumBnd[i + 1] - spectrumBnd[i]) * rate(i);
1633 transf += (spectrumBnd[i + 1] - spectrumBnd[i]) * transfer(i);
1634 transf2 += (spectrumBnd[i + 1] - spectrumBnd[i]) * transfer2(i);
1635 }
1636 ofl << "# Energy spectrum at time step " << globalTimeStep << ", total energy Ek=" << energy
1637 << ", dissipation rate=" << dissip << ", u_rms=" << urms << " (" << urms0
1638 << "), integral length scale=" << length << ", Kolmogorov time scale=" << sqrt(viscosity / dissip)
1639 << ", coupling=" << force << ", dedt=" << dedt << ", transfer=" << transf << ", transfer2=" << transf2
1640 << std::endl;
1641 ofl << "# 1: wave number k" << std::endl;
1642 ofl << "# 2: energy E(k)" << std::endl;
1643 ofl << "# 3: viscous dissipation rate D(k)" << std::endl;
1644 ofl << "# 4: transfer rate T(k)" << std::endl;
1645 ofl << "# 5: coupling rate Psi(k)" << std::endl;
1646 ofl << "# 6: flux -F(k)" << std::endl;
1647 ofl << "# 7: flux F(k)" << std::endl;
1648 ofl << "# 8: dE/dt(k)" << std::endl;
1649 ofl << "# 9: alternative transfer rate T(k)" << std::endl;
1650 ofl << "# 10: alternative flux -F(k)" << std::endl;
1651 ofl << "# 11: alternative flux F(k)" << std::endl;
1652 ofl << "# 12: number of samples in corresponding bin" << std::endl;
1653 ofl << "# 13: lower boundary of corresponding bin" << std::endl;
1654 ofl << "# 14: upper boundary of corresponding bin" << std::endl;
1655 ofl << "# 15: number of samples other spectrum: " << std::endl;
1656 ofl << "# 16: 4 * PI * k^2" << std::endl;
1657 ofl << "# 17: 4/3 * PI * (upperBndry^3 - lowerBndry^3)" << std::endl;
1658 for(MInt i = 0; i < noBins; i++) {
1659 MFloat k = F1B2 * (spectrumBnd[i + 1] + spectrumBnd[i]);
1660 MFloat Ek = spectrum(i);
1661 MFloat coupl = coupling(i);
1662 MFloat trans = transfer(i);
1663 MFloat trans2 = transfer2(i);
1664 ofl << k / k0 << " " << Ek << " " << F2 * POW2(k) * viscosity * Ek << " " << trans << " " << coupl << " "
1665 << flux(i, 0) << " " << flux(i, 1) << " " << rate(i) << " " << trans2 << " " << flux2(i, 0) << " "
1666 << flux2(i, 1) << " " << spectrumSum(i) << " " << spectrumBnd[i] << " " << spectrumBnd[i + 1] << " "
1667 << spectrumSum2(i) << " " << F4 * PI * POW2(k) << " "
1668 << F4B3 * PI * (POW3(spectrumBnd[i + 1]) - POW3(spectrumBnd[i])) << std::endl;
1669 }
1670 for(MInt i = 0; i < noBins / 2; i++) {
1671 const MInt i0 = 2 * i;
1672 const MInt i1 = 2 * i + 1;
1673 const MFloat k = F1B2 * (spectrumBnd[i1 + 1] + spectrumBnd[i0]);
1674 MFloat vol0 =
1675 F4B3 * PI * (POW3(spectrumBnd[i0 + 1]) - POW3(spectrumBnd[i0])) / (spectrumBnd[i0 + 1] - spectrumBnd[i0]);
1676 MFloat vol1 =
1677 F4B3 * PI * (POW3(spectrumBnd[i1 + 1]) - POW3(spectrumBnd[i1])) / (spectrumBnd[i1 + 1] - spectrumBnd[i1]);
1678 MFloat vol =
1679 F4B3 * PI * (POW3(spectrumBnd[i1 + 1]) - POW3(spectrumBnd[i0])) / (spectrumBnd[i1 + 1] - spectrumBnd[i0]);
1680 MFloat f0 = spectrumSum(i0) * vol / (vol0 * (spectrumSum(i0) + spectrumSum(i1)));
1681 MFloat f1 = spectrumSum(i1) * vol / (vol1 * (spectrumSum(i0) + spectrumSum(i1)));
1682
1683 ofl2 << k / k0 << " " << f1 * spectrum(i1) + f0 * spectrum(i0) << " "
1684 << F2 * POW2(k) * viscosity * (f1 * spectrum(i1) + f0 * spectrum(i0)) << " "
1685 << f1 * transfer(i1) + f0 * transfer(i0) << " " << f1 * coupling(i1) + f0 * coupling(i0) << " "
1686 << f1 * flux(i1, 0) + f0 * flux(i0, 0) << " " << f1 * flux(i1, 1) + f0 * flux(i0, 1) << " "
1687 << f1 * rate(i1) + f0 * rate(i0) << " " << f1 * transfer2(i1) + f0 * transfer2(i0) << " "
1688 << f1 * flux2(i1, 0) + f0 * flux2(i0, 0) << " " << f1 * flux2(i1, 1) + f0 * flux2(i0, 1) << " "
1689 << spectrumSum(i0) + spectrumSum(i1) << " " << spectrumBnd[i0] << " " << spectrumBnd[i1 + 1] << " "
1690 << spectrumSum2(i0) + spectrumSum2(i1) << " " << F4 * PI * POW2(k) << " "
1691 << F4B3 * PI * (POW3(spectrumBnd[i1 + 1]) - POW3(spectrumBnd[i0])) << std::endl;
1692 }
1693 ofl.close();
1694 ofl2.close();
1695 }
1696 }
1697
1698 // cleanup
1699 // fftw_cleanup(); problems occur when calling clean up function!!!
1700 }
1701
1702 const MFloat time3 = MPI_Wtime();
1703
1704 if(domainId == 0)
1705 std::cerr << "fft time " << time1 - time0 << " " << time2 - time1 << " " << time3 - time2 << std::endl;
1706
1707 m_log << "fft finished" << std::endl;
1708
1709 DEBUG("computeEnergySpectrum return", MAIA_DEBUG_TRACE_OUT);
1710}
1711
1712} // namespace math
1713} // namespace maia
1714#endif // MAIA_FFTW_H
Context::propertyExists
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
Definition: context.cpp:494
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
ScratchSpace::fill
void fill(T val)
fill the scratch with a given value
Definition: scratch.h:311
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
POW3
constexpr Real POW3(const Real x)
Definition: functions.h:123
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119
approx
MBool approx(const T &, const U &, const T)
Definition: functions.h:272
mMin
constexpr T mMin(const T &x, const T &y)
Definition: functions.h:90
mMax
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
globalTimeStep
MInt globalTimeStep
Definition: globalvariables.cpp:33
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
MInt
int32_t MInt
Definition: maiatypes.h:62
MFloat
double MFloat
Definition: maiatypes.h:52
MLong
int64_t MLong
Definition: maiatypes.h:64
MBool
bool MBool
Definition: maiatypes.h:58
MUlong
uint64_t MUlong
Definition: maiatypes.h:65
MPI_Recv
int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status *status, const MString &name, const MString &varname)
same as MPI_Recv
Definition: mpioverride.cpp:565
MPI_Issend
int MPI_Issend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Issend
Definition: mpioverride.cpp:523
MPI_Wait
int MPI_Wait(MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Wait
Definition: mpioverride.cpp:732
MPI_Waitall
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
Definition: mpioverride.cpp:766
MPI_Allreduce
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
Definition: mpioverride.cpp:952
MPI_Alltoall
int MPI_Alltoall(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Alltoall
Definition: mpioverride.cpp:1356
MPI_Comm_split
int MPI_Comm_split(MPI_Comm comm, int color, int key, MPI_Comm *newcomm, const MString &name, const MString &varname)
same as MPI_Comm_split, but updates the number of MPI communicators
Definition: mpioverride.cpp:307
MPI_Bcast
int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm, const MString &name, const MString &varname)
same as MPI_Bcast
Definition: mpioverride.cpp:1114
MPI_Allgather
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
Definition: mpioverride.cpp:1274
maia::math::randnormal
MFloat randnormal(MFloat mu, MFloat sigma)
Definition: maiafftw.h:502
maia::math::getFourierCoefficient2
std::complex< MFloat > * getFourierCoefficient2(MFloat *k, MFloat k0)
Definition: maiafftw.h:620
maia::math::triadImag
MFloat triadImag(fftw_complex &a, fftw_complex &b, fftw_complex &c)
Definition: maiafftw.h:61
maia::math::initFftFilter
void initFftFilter(MFloat *uPhysFieldCoarse, MFloat *vPhysFieldCoarse, MFloat *wPhysFieldCoarse, MInt lx, MInt ly, MInt lz, MInt lxCoarse, MInt lyCoarse, MInt lzCoarse, MInt noPeakModes, const MFloat Ma)
Generates a velocity field from Fourier-modes using FFTW The disturbances are generated and filtered ...
Definition: maiafftw.h:681
maia::math::initFft
void initFft(fftw_complex *uPhysField, fftw_complex *vPhysField, fftw_complex *wPhysField, MInt lx, MInt ly, MInt lz, MInt noPeakModes, const MFloat Ma)
Generates a velocity field from Fourier-modes using FFTW.
Definition: maiafftw.h:892
maia::math::getGlobalPosFFTW
MInt getGlobalPosFFTW(MInt i0, MInt i1, MInt i2, MInt ny, MInt nz)
Definition: maiafftw.h:59
maia::math::getFourierCoefficient
std::complex< MFloat > * getFourierCoefficient(MFloat *k, MFloat k0, const MFloat Ma)
Definition: maiafftw.h:548
maia::math::computeEnergySpectrum
void computeEnergySpectrum(MFloatScratchSpace &velocity, MIntScratchSpace &indices, const ptrdiff_t howmany, MInt locSize, MInt nx, MInt ny, MInt nz, MFloat viscosity, const MPI_Comm comm)
Compute energy spectrum on unity cube.
Definition: maiafftw.h:1051
maia
Namespace for auxiliary functions/classes.
Definition: acaobserverdatacollector.h:76
a
Definition: contexttypes.h:19