Skip to content
Snippets Groups Projects
kspace.cpp 32.11 KiB
/**********************************************************************
*                                                                     *
*  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                                 *
*                                                                     *
**********************************************************************/



#define NDEBUG

#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <cassert>
#include "kspace.hpp"
#include "scope_timer.hpp"
#include "shared_array.hpp"



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)
{
    TIMEZONE("kspace::kspace");
    /* 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];
                if (ii <= int(this->layout->sizes[0]/2))
                    this->ky[i] = this->dky*ii;
                else
                    this->ky[i] = this->dky*(ii - int(this->layout->sizes[0]));
            }
            for (i = 0; i<int(this->layout->sizes[1]); i++)
            {
                if (i <= int(this->layout->sizes[1]/2))
                    this->kz[i] = this->dkz*i;
                else
                    this->kz[i] = this->dkz*(i - int(this->layout->sizes[1]));
            }
            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);

    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);
    });

    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);
                kshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes*knorm;
                nshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes;
            }
    });

    // Merge results
    kshell_local_thread.mergeParallel();
    nshell_local_thread.mergeParallel();

    MPI_Allreduce(
            nshell_local_thread.getMasterData(),
            &this->nshell.front(),
            this->nshells,
            MPI_INT64_T, MPI_SUM, this->layout->comm);
    MPI_Allreduce(
            kshell_local_thread.getMasterData(),
            &this->kshell.front(),
            this->nshells,
            MPI_DOUBLE, MPI_SUM, this->layout->comm);
    for (int n=0; n<this->nshells; n++){
		if(this->nshell[n] != 0){
	        this->kshell[n] /= this->nshell[n];
		}
    }
}

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

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);
    if (this->nshells != int(dims[0]))
    {
        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;
}

template <field_backend be,
          kspace_dealias_type dt>
template <typename rnumber,
          field_components fc>
void kspace<be, dt>::low_pass(
        typename fftw_interface<rnumber>::complex *__restrict__ a,
        const double kmax)
{
    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)
                std::fill_n((rnumber*)(a + ncomp(fc)*cindex), 2*ncomp(fc), 0);
                });
}

/** \brief Filter a field using a ball shaped top hat filter.
 *
 *  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:
 *  \f[
 *       \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)
 *  \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,
        const double ell)
{
    const double prefactor0 = double(3) / pow(ell/2, 3);
    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 / 2;
                    double prefactor = prefactor0 / pow(k2, 1.5);
                    for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
                            prefactor *
                            (sin(argument) - argument * cos(argument)));
                }
                });
}

/** \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]
 */
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);
                }
                });
}

/** \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:
 *  \f{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)
 *  \f}
 *  and the corresponding expressions in Fourier space are:
 *  \f{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)
 *  \f}
 *
 *  \f$ k_\ell \f$ is given as a parameter, and then we use
 *  \f[
 *      \ell = \pi / k_\ell,
 *      \sigma_\ell = \pi / k_\ell
 *  \f]
 *
 *  For the Gaussian filter this is the same convention used in
 *  \cite Buzzicotti2017 .
 *
 *  See also `filter_calibrated_ell`.
 */
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,
                2*acos(0.)/wavenumber);
    }
    else if (filter_type == std::string("ball"))
    {
        this->template ball_filter<rnumber, fc>(
                a,
                2*acos(0.)/wavenumber);
    }
    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:
 *  \f{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)
 *  \f}
 *  and the corresponding expressions in Fourier space are:
 *  \f{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)
 *  \f}
 *
 *  \f$\sigma_\ell\f$ and \f$k_\ell\f$ are calibrated such that the energy of
 *  the large scales is approximately the same (within the inertial range)
 *  independently of the shape of the filter.
 *
 *  This was done by hand, see [INSERT CITATION HERE] for details, with the
 *  results:
 *
 *  \f[
 *      \sigma_\ell = 0.23 \ell,
 *      k_\ell = 2.8 / \ell
 *  \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)
{
    TIMEZONE("kspace::filter_calibrated_ell");
    if (filter_type == std::string("sharp_Fourier_sphere"))
    {
        this->template low_pass<rnumber, fc>(
                a,
                2.8 / ell);
    }
    else if (filter_type == std::string("Gauss"))
    {
        this->template Gauss_filter<rnumber, fc>(
                a,
                0.23*ell);
    }
    else if (filter_type == std::string("ball"))
    {
        this->template ball_filter<rnumber, fc>(
                a,
                ell);
    }
    return EXIT_SUCCESS;
}

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)
{
    TIMEZONE("kspace::dealias");
    switch(dt)
    {
        case TWO_THIRDS:
            this->low_pass<rnumber, fc>(a, this->kM);
            break;
        case SMOOTH:
            this->CLOOP(
                [&](ptrdiff_t cindex,
                    ptrdiff_t xindex,
                    ptrdiff_t yindex,
                    ptrdiff_t zindex){
                    double kk2 = (pow(this->kx[xindex]/this->kMx, 2) +
                                  pow(this->ky[yindex]/this->kMy, 2) +
                                  pow(this->kz[zindex]/this->kMz, 2));
                    double tval = exp(-36.0 * (pow(kk2, 18)));
                    for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= tval;
                });
            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)
        {
                    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];
                    }
           }
        }
    );
    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)
{
    TIMEZONE("field::cospectrum");
    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] * b[ncomp(fc)*cindex + j][0]) +
                        (a[ncomp(fc)*cindex + i][1] * b[ncomp(fc)*cindex + j][1]));
                }
            }
            });

    spec_local_thread.mergeParallel();

    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>
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();
                for (hsize_t i=0; i<ncomp(fc); i++){
                    L2_local[0] += nxmodes * (
                        (a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + i][0]) +
                        (a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + i][1]));
                }
            }
            });

    L2_local_thread.mergeParallel();

    double L2;
    MPI_Allreduce(
            L2_local_thread.getMasterData(),
            &L2,
            1,
            MPI_DOUBLE, MPI_SUM, this->layout->comm);
    return sqrt(L2 * this->dkx * this->dky * this->dkz);
}


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>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax);
template void kspace<FFTW, SMOOTH>::low_pass<float, THREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax);
template void kspace<FFTW, SMOOTH>::low_pass<float, THREExTHREE>(
        typename fftw_interface<float>::complex *__restrict__ a,
        const double kmax);

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

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);

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);

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);

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);

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);

template double kspace<FFTW, TWO_THIRDS>::L2norm<float, ONE>(
        const typename fftw_interface<float>::complex *__restrict__ a);
template double kspace<FFTW, TWO_THIRDS>::L2norm<float, THREE>(
        const typename fftw_interface<float>::complex *__restrict__ a);
template double kspace<FFTW, TWO_THIRDS>::L2norm<float, THREExTHREE>(
        const typename fftw_interface<float>::complex *__restrict__ a);
template double kspace<FFTW, TWO_THIRDS>::L2norm<double, ONE>(
        const typename fftw_interface<double>::complex *__restrict__ a);
template double kspace<FFTW, TWO_THIRDS>::L2norm<double, THREE>(
        const typename fftw_interface<double>::complex *__restrict__ a);
template double kspace<FFTW, TWO_THIRDS>::L2norm<double, THREExTHREE>(
        const typename fftw_interface<double>::complex *__restrict__ a);

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

template void kspace<FFTW, SMOOTH>::force_divfree<float>(
       typename fftw_interface<float>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::force_divfree<double>(
       typename fftw_interface<double>::complex *__restrict__ a);