kspace.cpp 39.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
/**********************************************************************
*                                                                     *
*  Copyright 2015 Max Planck Institute                                *
*                 for Dynamics and Self-Organization                  *
*                                                                     *
*  This file is part of bfps.                                         *
*                                                                     *
*  bfps is free software: you can redistribute it and/or modify       *
*  it under the terms of the GNU General Public License as published  *
*  by the Free Software Foundation, either version 3 of the License,  *
*  or (at your option) any later version.                             *
*                                                                     *
*  bfps is distributed in the hope that it will be useful,            *
*  but WITHOUT ANY WARRANTY; without even the implied warranty of     *
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
*  GNU General Public License for more details.                       *
*                                                                     *
*  You should have received a copy of the GNU General Public License  *
*  along with bfps.  If not, see <http://www.gnu.org/licenses/>       *
*                                                                     *
* Contact: Cristian.Lalescu@ds.mpg.de                                 *
*                                                                     *
**********************************************************************/


26

27
28
29
30
31
32
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <cassert>
#include "kspace.hpp"
#include "scope_timer.hpp"
33
#include "shared_array.hpp"
34
35
36
37
38
39
40
41
42
43

template <field_backend be,
          kspace_dealias_type dt>
template <field_components fc>
kspace<be, dt>::kspace(
        const field_layout<fc> *source_layout,
        const double DKX,
        const double DKY,
        const double DKZ)
{
44
    TIMEZONE("kspace::kspace");
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
    /* get layout */
    this->layout = new field_layout<ONE>(
            source_layout->sizes,
            source_layout->subsizes,
            source_layout->starts,
            source_layout->comm);

    /* store dk values */
    this->dkx = DKX;
    this->dky = DKY;
    this->dkz = DKZ;

    /* compute kx, ky, kz and compute kM values */
    switch(be)
    {
        case FFTW:
            this->kx.resize(this->layout->sizes[2]);
            this->ky.resize(this->layout->subsizes[0]);
            this->kz.resize(this->layout->sizes[1]);
            int i, ii;
            for (i = 0; i<int(this->layout->sizes[2]); i++)
                this->kx[i] = i*this->dkx;
            for (i = 0; i<int(this->layout->subsizes[0]); i++)
            {
                ii = i + this->layout->starts[0];
Cristian Lalescu's avatar
Cristian Lalescu committed
70
                if (ii <= int(this->layout->sizes[0]/2))
71
72
                    this->ky[i] = this->dky*ii;
                else
73
                    this->ky[i] = this->dky*(ii - int(this->layout->sizes[0]));
74
75
76
            }
            for (i = 0; i<int(this->layout->sizes[1]); i++)
            {
Cristian Lalescu's avatar
Cristian Lalescu committed
77
                if (i <= int(this->layout->sizes[1]/2))
78
79
                    this->kz[i] = this->dkz*i;
                else
Cristian Lalescu's avatar
Cristian Lalescu committed
80
                    this->kz[i] = this->dkz*(i - int(this->layout->sizes[1]));
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
            }
            switch(dt)
            {
                case TWO_THIRDS:
                    this->kMx = this->dkx*(int(2*(int(this->layout->sizes[2])-1)/3)-1);
                    this->kMy = this->dky*(int(this->layout->sizes[0] / 3)-1);
                    this->kMz = this->dkz*(int(this->layout->sizes[1] / 3)-1);
                    break;
                case SMOOTH:
                    this->kMx = this->dkx*(int(this->layout->sizes[2])-2);
                    this->kMy = this->dky*(int(this->layout->sizes[0] / 2)-1);
                    this->kMz = this->dkz*(int(this->layout->sizes[1] / 2)-1);
                    break;
            }
            break;
    }

    /* get global kM and dk */
    this->kM = this->kMx;
    if (this->kM < this->kMy) this->kM = this->kMy;
    if (this->kM < this->kMz) this->kM = this->kMz;
    this->kM2 = this->kM * this->kM;
    this->dk = this->dkx;
    if (this->dk > this->dky) this->dk = this->dky;
    if (this->dk > this->dkz) this->dk = this->dkz;
    this->dk2 = this->dk*this->dk;

    /* spectra stuff */
    this->nshells = int(this->kM / this->dk) + 2;
    this->kshell.resize(this->nshells, 0);
    this->nshell.resize(this->nshells, 0);
112
113
114
115
116
117
118

    shared_array<double> kshell_local_thread(this->nshells,[&](double* kshell_local){
        std::fill_n(kshell_local, this->nshells, 0);
    });
    shared_array<int64_t> nshell_local_thread(this->nshells,[&](int64_t* nshell_local){
        std::fill_n(nshell_local, this->nshells, 0);
    });
119

120
121
122
123
124
125
126
127
128
129
    this->CLOOP_K2_NXMODES(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2,
                int nxmodes){
            if (k2 < this->kM2)
            {
                double knorm = sqrt(k2);
130
131
                kshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes*knorm;
                nshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes;
132
            }
133
134
135
    });

    // Merge results
136
137
138
139

    kshell_local_thread.mergeParallel();
    nshell_local_thread.mergeParallel();

140
    MPI_Allreduce(
141
            nshell_local_thread.getMasterData(),
142
143
144
145
            &this->nshell.front(),
            this->nshells,
            MPI_INT64_T, MPI_SUM, this->layout->comm);
    MPI_Allreduce(
146
            kshell_local_thread.getMasterData(),
147
148
149
            &this->kshell.front(),
            this->nshells,
            MPI_DOUBLE, MPI_SUM, this->layout->comm);
Berenger Bramas's avatar
Berenger Bramas committed
150
151
152
153
154
    for (int n=0; n<this->nshells; n++){
		if(this->nshell[n] != 0){
	        this->kshell[n] /= this->nshell[n];
		}
    }
155
156
157
158
159
160
161
162
163
}

template <field_backend be,
          kspace_dealias_type dt>
kspace<be, dt>::~kspace()
{
    delete this->layout;
}

164
165
166
167
168
169
170
171
172
173
174
175
176
template <field_backend be,
          kspace_dealias_type dt>
int kspace<be, dt>::store(hid_t stat_file)
{
    TIMEZONE("kspace::store");
    assert(this->layout->myrank == 0);
    hsize_t dims[4];
    hid_t space, dset;
    // store kspace information
    dset = H5Dopen(stat_file, "/kspace/kshell", H5P_DEFAULT);
    space = H5Dget_space(dset);
    H5Sget_simple_extent_dims(space, dims, NULL);
    H5Sclose(space);
Cristian Lalescu's avatar
Cristian Lalescu committed
177
    if (this->nshells != int(dims[0]))
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
    {
        DEBUG_MSG(
                "ERROR: computed nshells %d not equal to data file nshells %d\n",
                this->nshells, dims[0]);
    }
    H5Dwrite(
            dset,
            H5T_NATIVE_DOUBLE,
            H5S_ALL,
            H5S_ALL,
            H5P_DEFAULT,
            &this->kshell.front());
    H5Dclose(dset);
    dset = H5Dopen(
            stat_file,
            "/kspace/nshell",
            H5P_DEFAULT);
    H5Dwrite(
            dset,
            H5T_NATIVE_INT64,
            H5S_ALL,
            H5S_ALL,
            H5P_DEFAULT,
            &this->nshell.front());
    H5Dclose(dset);
    dset = H5Dopen(stat_file, "/kspace/kM", H5P_DEFAULT);
    H5Dwrite(
            dset,
            H5T_NATIVE_DOUBLE,
            H5S_ALL,
            H5S_ALL,
            H5P_DEFAULT,
            &this->kM);
    H5Dclose(dset);
    dset = H5Dopen(stat_file, "/kspace/dk", H5P_DEFAULT);
    H5Dwrite(dset,
            H5T_NATIVE_DOUBLE,
            H5S_ALL,
            H5S_ALL,
            H5P_DEFAULT,
            &this->dk);
    H5Dclose(dset);
    return EXIT_SUCCESS;
}

223
224
225
226
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
Cristian Lalescu's avatar
Cristian Lalescu committed
227
228
229
void kspace<be, dt>::low_pass(
        typename fftw_interface<rnumber>::complex *__restrict__ a,
        const double kmax)
230
231
232
233
234
235
236
237
238
{
    const double km2 = kmax*kmax;
    this->CLOOP_K2(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2){
            if (k2 >= km2)
239
                std::fill_n((rnumber*)(a + ncomp(fc)*cindex), 2*ncomp(fc), 0);
240
241
242
                });
}

243
244
/** \brief Filter a field using a ball shaped top hat filter.
 *
Cristian Lalescu's avatar
Cristian Lalescu committed
245
246
247
248
249
250
 *  Filter's mathematical expression in real space is as follows:
 *  \f[
 *       \phi^b_\ell(r) =
 *           \frac{1}{\ell^3}\frac{6}{\pi} H(\ell/2 - r)
 *  \f]
 *  with the corresponding Fourier space expression:
251
 *  \f[
Cristian Lalescu's avatar
Cristian Lalescu committed
252
253
254
 *       \hat{\phi^b_\ell}(k) =
 *       \frac{3}{2(k\ell/2)^3}
 *       \left(2\sin (k \ell/2) - k \ell \cos (k \ell/2)\right)
255
256
257
258
259
260
261
262
 *  \f]
 */
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
void kspace<be, dt>::ball_filter(
        typename fftw_interface<rnumber>::complex *__restrict__ a,
Cristian Lalescu's avatar
Cristian Lalescu committed
263
        const double ell)
264
{
265
    const double prefactor0 = double(3) / pow(ell/2, 3);
266
267
268
269
270
271
272
273
    this->CLOOP_K2(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2){
                if (k2 > 0)
                {
274
275
                    double argument = sqrt(k2)*ell / 2;
                    double prefactor = prefactor0 / pow(k2, 1.5);
276
                    for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
Cristian Lalescu's avatar
Cristian Lalescu committed
277
278
279
                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
                            prefactor *
                            (sin(argument) - argument * cos(argument)));
280
281
282
283
                }
                });
}

Debarghya Banerjee's avatar
Debarghya Banerjee committed
284
285
286
287
288
/** \brief Filter a field using a M-filter to reproduce dissipation range.
 *
 *  Filter's  Fourier space expression:
 *  \f[
 *       \hat{\phi^M_\ell}(k) =
289
290
 *       \exp(-\frac{(3.54 k \ell)^{122 \ell^{0.0836}}}{2})
 *       \left( 1 + \frac{(k \eta/0.0636)^{3.44}}{1 + (k \eta/ 0.0621)^{3.44}} \right)^{1/2}
Debarghya Banerjee's avatar
Debarghya Banerjee committed
291
292
293
294
295
296
297
298
 *  \f]
 */
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
void kspace<be, dt>::general_M_filter(
        typename fftw_interface<rnumber>::complex *__restrict__ a,
299
        const double ell)
Debarghya Banerjee's avatar
Debarghya Banerjee committed
300
301
302
303
304
305
306
307
308
309
310
311
312
313
{
    const double prefactor0 = 1.0;
    this->CLOOP_K2(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2){
                if (k2 > 0)
                {
                    double argument = sqrt(k2)*ell;
                    double prefactor = prefactor0;
                    for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
314
                            prefactor * (exp(-0.5*pow((2.9*argument),(68.0*(pow(ell,0.74))))) * sqrt(1.0 + (pow((argument/0.06),3.8))/(1.0 + (pow((argument/0.057),3.8))))));
Debarghya Banerjee's avatar
Debarghya Banerjee committed
315
316
317
318
319
                }
                });
}


320
321
322
323
324
325
326
/** \brief Filter a field using a Gaussian kernel.
 *
 *  Filter's mathematical expression in Fourier space is as follows:
 *  \f[
 *      \hat{g}_\ell(\mathbf{k}) = \exp(-k^2 \sigma^2 / 2)
 *  \f]
 */
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
void kspace<be, dt>::Gauss_filter(
        typename fftw_interface<rnumber>::complex *__restrict__ a,
        const double sigma)
{
    const double prefactor = - sigma*sigma/2;
    this->CLOOP_K2(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2){
                {
                    for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= exp(prefactor*k2);
                }
                });
}

349
350
/** \brief Filter a field.
 *
351
352
 *  This is a wrapper that can choose between a sharp Fourier spherical filter,
 *  a Gaussian filter and a sharp real space spherical filter.
353
 *
Cristian Lalescu's avatar
Cristian Lalescu committed
354
 *  Filter expressions in real space are as follows:
355
356
357
358
359
360
 *  \f[
 *       \phi^b_\ell(r) =
 *          \frac{1}{\ell^3}\frac{6}{\pi} H(\ell/2 - r)
 *  \f]
 *  \f[
 *       \phi^g_\ell(r) =
Cristian Lalescu's avatar
Cristian Lalescu committed
361
 *           \frac{1}{\sigma_\ell^3}\frac{1}{(2\pi)^{3/2}}
362
363
364
365
 *           \exp\left(-\frac{1}{2}\left(\frac{r}{\sigma_\ell}\right)^2\right)
 *  \f]
 *  \f[
 *       \phi^s_\ell(r) =
Cristian Lalescu's avatar
Cristian Lalescu committed
366
367
 *           \frac{1}{2 \pi^2 r^3}
 *           \left(\sin k_\ell r - k_\ell r \cos k_\ell r\right)
368
 *  \f]
Cristian Lalescu's avatar
Cristian Lalescu committed
369
 *  and the corresponding expressions in Fourier space are:
370
371
 *  \f[
 *       \hat{\phi^b_\ell}(k) =
Cristian Lalescu's avatar
Cristian Lalescu committed
372
 *       \frac{3}{2(k\ell/2)^3}
373
374
375
376
377
378
379
380
381
 *       \left(2\sin (k \ell/2) - k \ell \cos (k \ell/2)\right)
 *  \f]
 *  \f[
 *       \hat{\phi^g_\ell}(k) =
 *       \exp\left(-\frac{1}{2}k^2 \sigma_\ell^2\right)
 *  \f]
 *  \f[
 *       \hat{\phi^s_\ell}(k) = H(k_\ell - k)
 *  \f]
Cristian Lalescu's avatar
Cristian Lalescu committed
382
383
 *
 *  \f$ k_\ell \f$ is given as a parameter, and then we use
384
 *  \f[
Cristian Lalescu's avatar
Cristian Lalescu committed
385
386
 *      \ell = \pi / k_\ell,
 *      \sigma_\ell = \pi / k_\ell
387
 *  \f]
Cristian Lalescu's avatar
Cristian Lalescu committed
388
389
390
391
392
 *
 *  For the Gaussian filter this is the same convention used in
 *  \cite Buzzicotti2017 .
 *
 *  See also `filter_calibrated_ell`.
393
 */
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
int kspace<be, dt>::filter(
        typename fftw_interface<rnumber>::complex *__restrict__ a,
        const double wavenumber,
        std::string filter_type)
{
    if (filter_type == std::string("sharp_Fourier_sphere"))
    {
        this->template low_pass<rnumber, fc>(
                a,
                wavenumber);
    }
    else if (filter_type == std::string("Gauss"))
    {
        this->template Gauss_filter<rnumber, fc>(
                a,
413
                2*acos(0.)/wavenumber);
414
    }
415
416
417
418
    else if (filter_type == std::string("ball"))
    {
        this->template ball_filter<rnumber, fc>(
                a,
Cristian Lalescu's avatar
Cristian Lalescu committed
419
420
                2*acos(0.)/wavenumber);
    }
Debarghya Banerjee's avatar
Debarghya Banerjee committed
421
422
423
424
    else if (filter_type == std::string("general_M"))
    {
        this->template general_M_filter<rnumber, fc>(
                a,
425
                2*acos(0.)/wavenumber);
Debarghya Banerjee's avatar
Debarghya Banerjee committed
426
    }
Cristian Lalescu's avatar
Cristian Lalescu committed
427
428
429
430
431
432
433
434
435
    return EXIT_SUCCESS;
}

/** \brief Filter a field.
 *
 *  This is a wrapper that can choose between a sharp Fourier spherical filter,
 *  a Gaussian filter and a sharp real space spherical filter.
 *
 *  Filter expressions in real space are as follows:
Cristian Lalescu's avatar
Cristian Lalescu committed
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
 *  \rst
 *  .. math::
 *      :nowrap:
 *
 *      \begin{eqnarray*}
 *          \phi^b_\ell(r) &=&
 *              \frac{1}{\ell^3}\frac{6}{\pi} H(\ell/2 - r) \\
 *          \phi^g_\ell(r) &=&
 *              \frac{1}{\sigma_\ell^3}\frac{1}{(2\pi)^{3/2}}
 *              \exp\left(-\frac{1}{2}\left(\frac{r}{\sigma_\ell}\right)^2\right) \\
 *          \phi^s_\ell(r) &=&
 *              \frac{1}{2 \pi^2 r^3}
 *              \left(\sin k_\ell r - k_\ell r \cos k_\ell r\right)
 *      \end{eqnarray*}
 *
 *  \endrst
 *
Cristian Lalescu's avatar
Cristian Lalescu committed
453
 *  and the corresponding expressions in Fourier space are:
Cristian Lalescu's avatar
Cristian Lalescu committed
454
455
456
457
458
459
460
461
462
463
464
465
466
467
 *  \rst
 *  .. math::
 *      :nowrap:
 *
 *      \begin{eqnarray*}
 *          \hat{\phi^b_\ell}(k) &=&
 *          \frac{3}{2(k\ell/2)^3}
 *          \left(2\sin (k \ell/2) - k \ell \cos (k \ell/2)\right) \\
 *          \hat{\phi^g_\ell}(k) &=&
 *          \exp\left(-\frac{1}{2}k^2 \sigma_\ell^2\right) \\
 *          \hat{\phi^s_\ell}(k) &=& H(k_\ell - k)
 *      \end{eqnarray*}
 *
 *  \endrst
Cristian Lalescu's avatar
Cristian Lalescu committed
468
469
 *
 *  \f$\sigma_\ell\f$ and \f$k_\ell\f$ are calibrated such that the energy of
Cristian Lalescu's avatar
Cristian Lalescu committed
470
 *  the large scales is approximately the same (within the inertial range)
Cristian Lalescu's avatar
Cristian Lalescu committed
471
472
473
474
475
476
 *  independently of the shape of the filter.
 *
 *  This was done by hand, see [INSERT CITATION HERE] for details, with the
 *  results:
 *
 *  \f[
Cristian Lalescu's avatar
Cristian Lalescu committed
477
478
 *      \sigma_\ell = 0.23 \ell,
 *      k_\ell = 2.8 / \ell
Cristian Lalescu's avatar
Cristian Lalescu committed
479
480
481
482
483
484
485
486
487
488
489
490
 *  \f]
 *
 */
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
int kspace<be, dt>::filter_calibrated_ell(
        typename fftw_interface<rnumber>::complex *__restrict__ a,
        const double ell,
        std::string filter_type)
{
491
    TIMEZONE("kspace::filter_calibrated_ell");
Cristian Lalescu's avatar
Cristian Lalescu committed
492
493
494
495
    if (filter_type == std::string("sharp_Fourier_sphere"))
    {
        this->template low_pass<rnumber, fc>(
                a,
Cristian Lalescu's avatar
Cristian Lalescu committed
496
                2.8 / ell);
Cristian Lalescu's avatar
Cristian Lalescu committed
497
498
499
500
501
    }
    else if (filter_type == std::string("Gauss"))
    {
        this->template Gauss_filter<rnumber, fc>(
                a,
Cristian Lalescu's avatar
Cristian Lalescu committed
502
                0.23*ell);
Cristian Lalescu's avatar
Cristian Lalescu committed
503
504
505
506
507
508
    }
    else if (filter_type == std::string("ball"))
    {
        this->template ball_filter<rnumber, fc>(
                a,
                ell);
509
    }
Debarghya Banerjee's avatar
Debarghya Banerjee committed
510
511
512
513
    else if (filter_type == std::string("general_M"))
    {
        this->template general_M_filter<rnumber, fc>(
                a,
514
                ell);
Debarghya Banerjee's avatar
Debarghya Banerjee committed
515
    }
516
517
518
    return EXIT_SUCCESS;
}

519
520
521
522
523
524
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
void kspace<be, dt>::dealias(typename fftw_interface<rnumber>::complex *__restrict__ a)
{
525
    TIMEZONE("kspace::dealias");
526
527
528
    switch(dt)
    {
        case TWO_THIRDS:
529
            this->low_pass<rnumber, fc>(a, this->kM);
530
531
            break;
        case SMOOTH:
532
            this->CLOOP(
533
534
535
                [&](ptrdiff_t cindex,
                    ptrdiff_t xindex,
                    ptrdiff_t yindex,
536
                    ptrdiff_t zindex){
Cristian Lalescu's avatar
Cristian Lalescu committed
537
538
539
                    double kk2 = (pow(this->kx[xindex]/this->kMx, 2) +
                                  pow(this->ky[yindex]/this->kMy, 2) +
                                  pow(this->kz[zindex]/this->kMz, 2));
540
                    double tval = exp(-36.0 * (pow(kk2, 18)));
541
                    for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
542
                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= tval;
543
                });
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
            break;
    }
}

template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber>
void kspace<be, dt>::force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a)
{
    TIMEZONE("kspace::force_divfree");
    this->CLOOP_K2(
                [&](ptrdiff_t cindex,
                    ptrdiff_t xindex,
                    ptrdiff_t yindex,
                    ptrdiff_t zindex,
                    double k2){
                if (k2 > 0)
        {
562
563
564
565
566
567
568
569
570
571
572
573
574
575
                    typename fftw_interface<rnumber>::complex tval;
                    tval[0] = (this->kx[xindex]*((*(a + cindex*3  ))[0]) +
                               this->ky[yindex]*((*(a + cindex*3+1))[0]) +
                               this->kz[zindex]*((*(a + cindex*3+2))[0]) ) / k2;
                    tval[1] = (this->kx[xindex]*((*(a + cindex*3  ))[1]) +
                               this->ky[yindex]*((*(a + cindex*3+1))[1]) +
                               this->kz[zindex]*((*(a + cindex*3+2))[1]) ) / k2;
                    for (int imag_part=0; imag_part<2; imag_part++)
                    {
                        a[cindex*3  ][imag_part] -= tval[imag_part]*this->kx[xindex];
                        a[cindex*3+1][imag_part] -= tval[imag_part]*this->ky[yindex];
                        a[cindex*3+2][imag_part] -= tval[imag_part]*this->kz[zindex];
                    }
           }
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
        }
    );
    if (this->layout->myrank == this->layout->rank[0][0])
        std::fill_n((rnumber*)(a), 6, 0.0);
}

template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
void kspace<be, dt>::cospectrum(
        const rnumber(* __restrict a)[2],
        const rnumber(* __restrict b)[2],
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset)
{
593
    TIMEZONE("field::cospectrum2");
Berenger Bramas's avatar
Berenger Bramas committed
594
    shared_array<double> spec_local_thread(this->nshells*ncomp(fc)*ncomp(fc),[&](double* spec_local){
Berenger Bramas's avatar
Berenger Bramas committed
595
        std::fill_n(spec_local, this->nshells*ncomp(fc)*ncomp(fc), 0);
596
597
    });

598
599
600
601
602
603
604
605
606
    this->CLOOP_K2_NXMODES(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2,
                int nxmodes){
            if (k2 <= this->kM2)
            {
607
                double* spec_local = spec_local_thread.getMine();
608
609
                int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
                for (hsize_t i=0; i<ncomp(fc); i++)
610
                for (hsize_t j=0; j<ncomp(fc); j++){
611
                    spec_local[tmp_int + i*ncomp(fc)+j] += nxmodes * (
612
613
614
                        (a[ncomp(fc)*cindex + i][0] * b[ncomp(fc)*cindex + j][0]) +
                        (a[ncomp(fc)*cindex + i][1] * b[ncomp(fc)*cindex + j][1]));
                }
615
616
            }
            });
617
618
619

    spec_local_thread.mergeParallel();

620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
    std::vector<double> spec;
    spec.resize(this->nshells*ncomp(fc)*ncomp(fc), 0);
    MPI_Allreduce(
            spec_local_thread.getMasterData(),
            &spec.front(),
            spec.size(),
            MPI_DOUBLE, MPI_SUM, this->layout->comm);
    if (this->layout->myrank == 0)
    {
        hid_t dset, wspace, mspace;
        hsize_t count[(ndim(fc)-2)*2], offset[(ndim(fc)-2)*2], dims[(ndim(fc)-2)*2];
        dset = H5Dopen(group, ("spectra/" + dset_name).c_str(), H5P_DEFAULT);
        wspace = H5Dget_space(dset);
        H5Sget_simple_extent_dims(wspace, dims, NULL);
        switch (fc)
        {
            case THREExTHREE:
                offset[4] = 0;
                offset[5] = 0;
                count[4] = 3;
                count[5] = 3;
            case THREE:
                offset[2] = 0;
                offset[3] = 0;
                count[2] = 3;
                count[3] = 3;
            default:
                offset[0] = toffset;
                offset[1] = 0;
                count[0] = 1;
                count[1] = this->nshells;
        }
        mspace = H5Screate_simple((ndim(fc)-2)*2, count, NULL);
        H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
        H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, &spec.front());
        H5Sclose(wspace);
        H5Sclose(mspace);
        H5Dclose(dset);
    }
}

template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
void kspace<be, dt>::cospectrum(
        const rnumber(* __restrict a)[2],
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset)
{
    TIMEZONE("field::cospectrum1");
    shared_array<double> spec_local_thread(this->nshells*ncomp(fc)*ncomp(fc),[&](double* spec_local){
        std::fill_n(spec_local, this->nshells*ncomp(fc)*ncomp(fc), 0);
    });

    this->CLOOP_K2_NXMODES(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2,
                int nxmodes){
            if (k2 <= this->kM2)
            {
                double* spec_local = spec_local_thread.getMine();
                int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
                for (hsize_t i=0; i<ncomp(fc); i++)
                for (hsize_t j=0; j<ncomp(fc); j++){
                    spec_local[tmp_int + i*ncomp(fc)+j] += nxmodes * (
                        (a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + j][0]) +
                        (a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + j][1]));
                }
            }
            });

    spec_local_thread.mergeParallel();

698
699
    std::vector<double> spec;
    spec.resize(this->nshells*ncomp(fc)*ncomp(fc), 0);
700
    MPI_Allreduce(
701
            spec_local_thread.getMasterData(),
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
            &spec.front(),
            spec.size(),
            MPI_DOUBLE, MPI_SUM, this->layout->comm);
    if (this->layout->myrank == 0)
    {
        hid_t dset, wspace, mspace;
        hsize_t count[(ndim(fc)-2)*2], offset[(ndim(fc)-2)*2], dims[(ndim(fc)-2)*2];
        dset = H5Dopen(group, ("spectra/" + dset_name).c_str(), H5P_DEFAULT);
        wspace = H5Dget_space(dset);
        H5Sget_simple_extent_dims(wspace, dims, NULL);
        switch (fc)
        {
            case THREExTHREE:
                offset[4] = 0;
                offset[5] = 0;
717
718
                count[4] = 3;
                count[5] = 3;
719
720
721
            case THREE:
                offset[2] = 0;
                offset[3] = 0;
722
723
                count[2] = 3;
                count[3] = 3;
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
            default:
                offset[0] = toffset;
                offset[1] = 0;
                count[0] = 1;
                count[1] = this->nshells;
        }
        mspace = H5Screate_simple((ndim(fc)-2)*2, count, NULL);
        H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
        H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, &spec.front());
        H5Sclose(wspace);
        H5Sclose(mspace);
        H5Dclose(dset);
    }
}

Cristian Lalescu's avatar
Cristian Lalescu committed
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
double kspace<be, dt>::L2norm(
        const rnumber(* __restrict a)[2])
{
    TIMEZONE("field::L2norm");
    shared_array<double> L2_local_thread(1,[&](double* spec_local){
        std::fill_n(spec_local, 1, 0);
    });

    this->CLOOP_K2_NXMODES(
            [&](ptrdiff_t cindex,
                ptrdiff_t xindex,
                ptrdiff_t yindex,
                ptrdiff_t zindex,
                double k2,
                int nxmodes){
            {
                double* L2_local = L2_local_thread.getMine();
760
                for (hsize_t i=0; i<ncomp(fc); i++){
Cristian Lalescu's avatar
Cristian Lalescu committed
761
                    L2_local[0] += nxmodes * (
762
763
                        (a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + i][0]) +
                        (a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + i][1]));
Cristian Lalescu's avatar
Cristian Lalescu committed
764
765
766
767
768
769
770
771
772
773
774
775
                }
            }
            });

    L2_local_thread.mergeParallel();

    double L2;
    MPI_Allreduce(
            L2_local_thread.getMasterData(),
            &L2,
            1,
            MPI_DOUBLE, MPI_SUM, this->layout->comm);
776
    return sqrt(L2 * this->dkx * this->dky * this->dkz);
Cristian Lalescu's avatar
Cristian Lalescu committed
777
778
}

779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803

template class kspace<FFTW, TWO_THIRDS>;
template class kspace<FFTW, SMOOTH>;

template kspace<FFTW, TWO_THIRDS>::kspace<>(
        const field_layout<ONE> *,
        const double, const double, const double);
template kspace<FFTW, TWO_THIRDS>::kspace<>(
        const field_layout<THREE> *,
        const double, const double, const double);
template kspace<FFTW, TWO_THIRDS>::kspace<>(
        const field_layout<THREExTHREE> *,
        const double, const double, const double);

template kspace<FFTW, SMOOTH>::kspace<>(
        const field_layout<ONE> *,
        const double, const double, const double);
template kspace<FFTW, SMOOTH>::kspace<>(
        const field_layout<THREE> *,
        const double, const double, const double);
template kspace<FFTW, SMOOTH>::kspace<>(
        const field_layout<THREExTHREE> *,
        const double, const double, const double);

template void kspace<FFTW, SMOOTH>::low_pass<float, ONE>(
804
        typename fftw_interface<float>::complex *__restrict__ a,
805
806
        const double kmax);
template void kspace<FFTW, SMOOTH>::low_pass<float, THREE>(
807
        typename fftw_interface<float>::complex *__restrict__ a,
808
809
        const double kmax);
template void kspace<FFTW, SMOOTH>::low_pass<float, THREExTHREE>(
810
        typename fftw_interface<float>::complex *__restrict__ a,
811
812
813
        const double kmax);

template void kspace<FFTW, SMOOTH>::low_pass<double, ONE>(
814
        typename fftw_interface<double>::complex *__restrict__ a,
815
816
        const double kmax);
template void kspace<FFTW, SMOOTH>::low_pass<double, THREE>(
817
        typename fftw_interface<double>::complex *__restrict__ a,
818
819
        const double kmax);
template void kspace<FFTW, SMOOTH>::low_pass<double, THREExTHREE>(
820
        typename fftw_interface<double>::complex *__restrict__ a,
821
822
        const double kmax);

Cristian Lalescu's avatar
Cristian Lalescu committed
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
template void kspace<FFTW, SMOOTH>::Gauss_filter<float, ONE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax);
template void kspace<FFTW, SMOOTH>::Gauss_filter<float, THREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax);
template void kspace<FFTW, SMOOTH>::Gauss_filter<float, THREExTHREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax);

template void kspace<FFTW, SMOOTH>::Gauss_filter<double, ONE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax);
template void kspace<FFTW, SMOOTH>::Gauss_filter<double, THREE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax);
template void kspace<FFTW, SMOOTH>::Gauss_filter<double, THREExTHREE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax);

843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
template int kspace<FFTW, SMOOTH>::filter<float, ONE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter<float, THREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter<float, THREExTHREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);

template int kspace<FFTW, SMOOTH>::filter<double, ONE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter<double, THREE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter<double, THREExTHREE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);

Cristian Lalescu's avatar
Cristian Lalescu committed
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<float, ONE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<float, THREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<float, THREExTHREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);

template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<double, ONE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<double, THREE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<double, THREExTHREE>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const double kmax,
        std::string filter_type);

895
896
897
898
899
900
901
902
903
904
905
906
907
908
template void kspace<FFTW, SMOOTH>::dealias<float, ONE>(
        typename fftw_interface<float>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::dealias<float, THREE>(
        typename fftw_interface<float>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::dealias<float, THREExTHREE>(
        typename fftw_interface<float>::complex *__restrict__ a);

template void kspace<FFTW, SMOOTH>::dealias<double, ONE>(
        typename fftw_interface<double>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::dealias<double, THREE>(
        typename fftw_interface<double>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::dealias<double, THREExTHREE>(
        typename fftw_interface<double>::complex *__restrict__ a);

909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, ONE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const typename fftw_interface<float>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const typename fftw_interface<float>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREExTHREE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const typename fftw_interface<float>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, ONE>(
        const typename fftw_interface<double>::complex *__restrict__ a,
        const typename fftw_interface<double>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREE>(
        const typename fftw_interface<double>::complex *__restrict__ a,
        const typename fftw_interface<double>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREExTHREE>(
        const typename fftw_interface<double>::complex *__restrict__ a,
        const typename fftw_interface<double>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);

template void kspace<FFTW, SMOOTH>::cospectrum<float, ONE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const typename fftw_interface<float>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const typename fftw_interface<float>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREExTHREE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const typename fftw_interface<float>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, SMOOTH>::cospectrum<double, ONE>(
        const typename fftw_interface<double>::complex *__restrict__ a,
        const typename fftw_interface<double>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREE>(
        const typename fftw_interface<double>::complex *__restrict__ a,
        const typename fftw_interface<double>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREExTHREE>(
        const typename fftw_interface<double>::complex *__restrict__ a,
        const typename fftw_interface<double>::complex *__restrict__ b,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);

983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, ONE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREExTHREE>(
        const typename fftw_interface<float>::complex *__restrict__ a,
        const hid_t group,
        const std::string dset_name,
        const hsize_t toffset);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, ONE>(
        const typename fftw_interface<double>::complex *__restrict__ a,
        const hid_t group,
For faster browsing, not all history is shown. View entire blame