Skip to content
Snippets Groups Projects
fluid_solver.cpp 29.07 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                                 *
*                                                                     *
**********************************************************************/

// code is generally compiled via setuptools, therefore NDEBUG is present
//#ifdef NDEBUG
//#undef NDEBUG
//#endif//NDEBUG

#include <cassert>
#include <cmath>
#include <cstring>
#include "fluid_solver.hpp"
#include "fftw_tools.hpp"



template <class rnumber>
void fluid_solver<rnumber>::impose_zero_modes()
{
    if (this->cd->myrank == this->cd->rank[0])
    {
        std::fill_n((rnumber*)(this->cu), 6, 0.0);
        std::fill_n((rnumber*)(this->cv[0]), 6, 0.0);
        std::fill_n((rnumber*)(this->cv[1]), 6, 0.0);
        std::fill_n((rnumber*)(this->cv[2]), 6, 0.0);
    }
}
/*****************************************************************************/
/* macro for specializations to numeric types compatible with FFTW           */

#define FLUID_SOLVER_DEFINITIONS(FFTW, R, MPI_RNUM, MPI_CNUM) \
 \
template<> \
fluid_solver<R>::fluid_solver( \
        const char *NAME, \
        int nx, \
        int ny, \
        int nz, \
        double DKX, \
        double DKY, \
        double DKZ, \
        int DEALIAS_TYPE, \
        unsigned FFTW_PLAN_RIGOR) : fluid_solver_base<R>( \
                NAME, \
                nx , ny , nz, \
                DKX, DKY, DKZ, \
                DEALIAS_TYPE, \
                FFTW_PLAN_RIGOR) \
{ \
    this->cvorticity = FFTW(alloc_complex)(this->cd->local_size);\
    this->cvelocity  = FFTW(alloc_complex)(this->cd->local_size);\
    this->rvorticity = FFTW(alloc_real)(this->cd->local_size*2);\
    this->rvelocity  = (R*)(this->cvelocity);\
    /*this->rvelocity  = FFTW(alloc_real)(this->cd->local_size*2);*/\
 \
    this->ru = this->rvelocity;\
    this->cu = this->cvelocity;\
 \
    this->rv[0] = this->rvorticity;\
    this->rv[3] = this->rvorticity;\
    this->cv[0] = this->cvorticity;\
    this->cv[3] = this->cvorticity;\
 \
    this->cv[1] = FFTW(alloc_complex)(this->cd->local_size);\
    this->cv[2] = this->cv[1];\
    this->rv[1] = FFTW(alloc_real)(this->cd->local_size*2);\
    this->rv[2] = this->rv[1];\
 \
    this->c2r_vorticity = new FFTW(plan);\
    this->r2c_vorticity = new FFTW(plan);\
    this->c2r_velocity  = new FFTW(plan);\
    this->r2c_velocity  = new FFTW(plan);\
 \
    ptrdiff_t sizes[] = {nz, \
                         ny, \
                         nx};\
 \
    *(FFTW(plan)*)this->c2r_vorticity = FFTW(mpi_plan_many_dft_c2r)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->cvorticity, this->rvorticity, \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_IN); \
 \
    *(FFTW(plan)*)this->r2c_vorticity = FFTW(mpi_plan_many_dft_r2c)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->rvorticity, this->cvorticity, \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_OUT); \
 \
    *(FFTW(plan)*)this->c2r_velocity = FFTW(mpi_plan_many_dft_c2r)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->cvelocity, this->rvelocity, \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_IN); \
 \
    *(FFTW(plan)*)this->r2c_velocity = FFTW(mpi_plan_many_dft_r2c)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->rvelocity, this->cvelocity, \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_OUT); \
 \
    this->uc2r = this->c2r_velocity;\
    this->ur2c = this->r2c_velocity;\
    this->vc2r[0] = this->c2r_vorticity;\
    this->vr2c[0] = this->r2c_vorticity;\
 \
    this->vc2r[1] = new FFTW(plan);\
    this->vr2c[1] = new FFTW(plan);\
    this->vc2r[2] = new FFTW(plan);\
    this->vr2c[2] = new FFTW(plan);\
 \
    *(FFTW(plan)*)(this->vc2r[1]) = FFTW(mpi_plan_many_dft_c2r)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->cv[1], this->rv[1], \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_IN); \
 \
    *(FFTW(plan)*)this->vc2r[2] = FFTW(mpi_plan_many_dft_c2r)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->cv[2], this->rv[2], \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_IN); \
 \
    *(FFTW(plan)*)this->vr2c[1] = FFTW(mpi_plan_many_dft_r2c)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->rv[1], this->cv[1], \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_OUT); \
 \
    *(FFTW(plan)*)this->vr2c[2] = FFTW(mpi_plan_many_dft_r2c)( \
            3, sizes, 3, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, \
            this->rv[2], this->cv[2], \
            MPI_COMM_WORLD, this->fftw_plan_rigor | FFTW_MPI_TRANSPOSED_OUT); \
 \
    /* ``physical'' parameters etc, initialized here just in case */ \
 \
    this->nu = 0.1; \
    this->fmode = 1; \
    this->famplitude = 1.0; \
    this->fk0  = 0; \
    this->fk1 = 3.0; \
    /* initialization of fields must be done AFTER planning */ \
    std::fill_n((R*)this->cvorticity, this->cd->local_size*2, 0.0); \
    std::fill_n((R*)this->cvelocity, this->cd->local_size*2, 0.0); \
    std::fill_n(this->rvorticity, this->cd->local_size*2, 0.0); \
    std::fill_n((R*)this->cv[1], this->cd->local_size*2, 0.0); \
    std::fill_n(this->rv[1], this->cd->local_size*2, 0.0); \
    std::fill_n(this->rv[2], this->cd->local_size*2, 0.0); \
} \
 \
template<> \
fluid_solver<R>::~fluid_solver() \
{ \
    DEBUG_MSG("entered ~fluid_solver\n"); \
    FFTW(destroy_plan)(*(FFTW(plan)*)this->c2r_vorticity);\
    FFTW(destroy_plan)(*(FFTW(plan)*)this->r2c_vorticity);\
    FFTW(destroy_plan)(*(FFTW(plan)*)this->c2r_velocity );\
    FFTW(destroy_plan)(*(FFTW(plan)*)this->r2c_velocity );\
    FFTW(destroy_plan)(*(FFTW(plan)*)this->vc2r[1]);\
    FFTW(destroy_plan)(*(FFTW(plan)*)this->vr2c[1]);\
    FFTW(destroy_plan)(*(FFTW(plan)*)this->vc2r[2]);\
    FFTW(destroy_plan)(*(FFTW(plan)*)this->vr2c[2]);\
    DEBUG_MSG("destroyed plans\n"); \
 \
    delete (FFTW(plan)*)this->c2r_vorticity;\
    delete (FFTW(plan)*)this->r2c_vorticity;\
    delete (FFTW(plan)*)this->c2r_velocity ;\
    delete (FFTW(plan)*)this->r2c_velocity ;\
    delete (FFTW(plan)*)this->vc2r[1];\
    delete (FFTW(plan)*)this->vr2c[1];\
    delete (FFTW(plan)*)this->vc2r[2];\
    delete (FFTW(plan)*)this->vr2c[2];\
    DEBUG_MSG("deleted plans\n"); \
 \
    FFTW(free)(this->cv[1]);\
    FFTW(free)(this->rv[1]);\
    FFTW(free)(this->cvorticity);\
    FFTW(free)(this->rvorticity);\
    FFTW(free)(this->cvelocity);\
    DEBUG_MSG("freed arrays\n"); \
    DEBUG_MSG("exiting ~fluid_solver\n"); \
} \
 \
template<> \
void fluid_solver<R>::compute_vorticity() \
{ \
    ptrdiff_t tindex; \
    CLOOP_K2( \
            tindex = 3*cindex; \
            if (k2 <= this->kM2) \
            { \
                this->cvorticity[tindex+0][0] = -(this->ky[yindex]*this->cu[tindex+2][1] - this->kz[zindex]*this->cu[tindex+1][1]); \
                this->cvorticity[tindex+1][0] = -(this->kz[zindex]*this->cu[tindex+0][1] - this->kx[xindex]*this->cu[tindex+2][1]); \
                this->cvorticity[tindex+2][0] = -(this->kx[xindex]*this->cu[tindex+1][1] - this->ky[yindex]*this->cu[tindex+0][1]); \
                this->cvorticity[tindex+0][1] =  (this->ky[yindex]*this->cu[tindex+2][0] - this->kz[zindex]*this->cu[tindex+1][0]); \
                this->cvorticity[tindex+1][1] =  (this->kz[zindex]*this->cu[tindex+0][0] - this->kx[xindex]*this->cu[tindex+2][0]); \
                this->cvorticity[tindex+2][1] =  (this->kx[xindex]*this->cu[tindex+1][0] - this->ky[yindex]*this->cu[tindex+0][0]); \
            } \
            else \
                std::fill_n((R*)(this->cvorticity+tindex), 6, 0.0); \
            ); \
    this->symmetrize(this->cvorticity, 3); \
} \
 \
template<> \
void fluid_solver<R>::compute_velocity(FFTW(complex) *vorticity) \
{ \
    ptrdiff_t tindex; \
    CLOOP_K2( \
            tindex = 3*cindex; \
            if (k2 <= this->kM2 && k2 > 0) \
            { \
                this->cu[tindex+0][0] = -(this->ky[yindex]*vorticity[tindex+2][1] - this->kz[zindex]*vorticity[tindex+1][1]) / k2; \
                this->cu[tindex+1][0] = -(this->kz[zindex]*vorticity[tindex+0][1] - this->kx[xindex]*vorticity[tindex+2][1]) / k2; \
                this->cu[tindex+2][0] = -(this->kx[xindex]*vorticity[tindex+1][1] - this->ky[yindex]*vorticity[tindex+0][1]) / k2; \
                this->cu[tindex+0][1] =  (this->ky[yindex]*vorticity[tindex+2][0] - this->kz[zindex]*vorticity[tindex+1][0]) / k2; \
                this->cu[tindex+1][1] =  (this->kz[zindex]*vorticity[tindex+0][0] - this->kx[xindex]*vorticity[tindex+2][0]) / k2; \
                this->cu[tindex+2][1] =  (this->kx[xindex]*vorticity[tindex+1][0] - this->ky[yindex]*vorticity[tindex+0][0]) / k2; \
            } \
            else \
                std::fill_n((R*)(this->cu+tindex), 6, 0.0); \
            ); \
    /*this->symmetrize(this->cu, 3);*/ \
} \
 \
template<> \
void fluid_solver<R>::ift_velocity() \
{ \
    FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
} \
 \
template<> \
void fluid_solver<R>::ift_vorticity() \
{ \
    std::fill_n(this->rvorticity, this->cd->local_size*2, 0.0); \
    FFTW(execute)(*((FFTW(plan)*)this->c2r_vorticity )); \
} \
 \
template<> \
void fluid_solver<R>::dft_velocity() \
{ \
    FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
} \
 \
template<> \
void fluid_solver<R>::dft_vorticity() \
{ \
    std::fill_n((R*)this->cvorticity, this->cd->local_size*2, 0.0); \
    FFTW(execute)(*((FFTW(plan)*)this->r2c_vorticity )); \
} \
 \
template<> \
void fluid_solver<R>::add_forcing(\
        FFTW(complex) *acc_field, FFTW(complex) *vort_field, R factor) \
{ \
    if (strcmp(this->forcing_type, "none") == 0) \
        return; \
    if (strcmp(this->forcing_type, "Kolmogorov") == 0) \
    { \
        ptrdiff_t cindex; \
        if (this->cd->myrank == this->cd->rank[this->fmode]) \
        { \
            cindex = ((this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]*3; \
            acc_field[cindex+2][0] -= this->famplitude*factor/2; \
        } \
        if (this->cd->myrank == this->cd->rank[this->cd->sizes[0] - this->fmode]) \
        { \
            cindex = ((this->cd->sizes[0] - this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]*3; \
            acc_field[cindex+2][0] -= this->famplitude*factor/2; \
        } \
        return; \
    } \
    if (strcmp(this->forcing_type, "linear") == 0) \
    { \
        double knorm; \
        CLOOP( \
                knorm = sqrt(this->kx[xindex]*this->kx[xindex] + \
                             this->ky[yindex]*this->ky[yindex] + \
                             this->kz[zindex]*this->kz[zindex]); \
                if ((this->fk0 <= knorm) && \
                    (this->fk1 >= knorm)) \
                    for (int c=0; c<3; c++) \
                    for (int i=0; i<2; i++) \
                        acc_field[cindex*3+c][i] += this->famplitude*vort_field[cindex*3+c][i]*factor; \
             ); \
        return; \
    } \
} \
 \
template<> \
void fluid_solver<R>::omega_nonlin( \
        int src) \
{ \
    assert(src >= 0 && src < 3); \
    this->compute_velocity(this->cv[src]); \
    /* get fields from Fourier space to real space */ \
    FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity ));  \
    FFTW(execute)(*((FFTW(plan)*)this->vc2r[src]));      \
    /* compute cross product $u \times \omega$, and normalize */ \
    R tmp[3][2]; \
    ptrdiff_t tindex; \
    RLOOP ( \
            this, \
            tindex = 3*rindex; \
            for (int cc=0; cc<3; cc++) \
                tmp[cc][0] = (this->ru[tindex+(cc+1)%3]*this->rv[src][tindex+(cc+2)%3] - \
                              this->ru[tindex+(cc+2)%3]*this->rv[src][tindex+(cc+1)%3]); \
            for (int cc=0; cc<3; cc++) \
                this->ru[(3*rindex)+cc] = tmp[cc][0] / this->normalization_factor; \
            ); \
    /* go back to Fourier space */ \
    this->clean_up_real_space(this->ru, 3); \
    FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
    this->dealias(this->cu, 3); \
    /* $\imath k \times Fourier(u \times \omega)$ */ \
    CLOOP( \
            tindex = 3*cindex; \
            { \
                tmp[0][0] = -(this->ky[yindex]*this->cu[tindex+2][1] - this->kz[zindex]*this->cu[tindex+1][1]); \
                tmp[1][0] = -(this->kz[zindex]*this->cu[tindex+0][1] - this->kx[xindex]*this->cu[tindex+2][1]); \
                tmp[2][0] = -(this->kx[xindex]*this->cu[tindex+1][1] - this->ky[yindex]*this->cu[tindex+0][1]); \
                tmp[0][1] =  (this->ky[yindex]*this->cu[tindex+2][0] - this->kz[zindex]*this->cu[tindex+1][0]); \
                tmp[1][1] =  (this->kz[zindex]*this->cu[tindex+0][0] - this->kx[xindex]*this->cu[tindex+2][0]); \
                tmp[2][1] =  (this->kx[xindex]*this->cu[tindex+1][0] - this->ky[yindex]*this->cu[tindex+0][0]); \
            } \
            for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++) \
                this->cu[tindex+cc][i] = tmp[cc][i]; \
            ); \
    this->add_forcing(this->cu, this->cv[src], 1.0); \
    this->force_divfree(this->cu); \
} \
 \
template<> \
void fluid_solver<R>::step(double dt) \
{ \
    double factor0, factor1; \
    std::fill_n((R*)this->cv[1], this->cd->local_size*2, 0.0); \
    this->omega_nonlin(0); \
    CLOOP_K2( \
            if (k2 <= this->kM2) \
            { \
                factor0 = exp(-this->nu * k2 * dt); \
                for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++) \
                    this->cv[1][3*cindex+cc][i] = (this->cv[0][3*cindex+cc][i] + \
                                                   dt*this->cu[3*cindex+cc][i])*factor0; \
            } \
            ); \
 \
    this->omega_nonlin(1); \
    CLOOP_K2( \
            if (k2 <= this->kM2) \
            { \
                factor0 = exp(-this->nu * k2 * dt/2); \
                factor1 = exp( this->nu * k2 * dt/2); \
                for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++) \
                    this->cv[2][3*cindex+cc][i] = (3*this->cv[0][3*cindex+cc][i]*factor0 + \
                                                   (this->cv[1][3*cindex+cc][i] + \
                                                    dt*this->cu[3*cindex+cc][i])*factor1)*0.25; \
            } \
            ); \
 \
    this->omega_nonlin(2); \
    CLOOP_K2( \
            if (k2 <= this->kM2) \
            { \
                factor0 = exp(-this->nu * k2 * dt * 0.5); \
                for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++) \
                    this->cv[3][3*cindex+cc][i] = (this->cv[0][3*cindex+cc][i]*factor0 + \
                                                   2*(this->cv[2][3*cindex+cc][i] + \
                                                      dt*this->cu[3*cindex+cc][i]))*factor0/3; \
            } \
            ); \
 \
    this->force_divfree(this->cvorticity); \
    this->iteration++; \
} \
 \
template<> \
int fluid_solver<R>::read(char field, char representation) \
{ \
    char fname[512]; \
    int read_result; \
    if (field == 'v') \
    { \
        if (representation == 'c') \
        { \
            this->fill_up_filename("cvorticity", fname); \
            read_result = this->cd->read(fname, (void*)this->cvorticity); \
            if (read_result != EXIT_SUCCESS) \
                return read_result; \
        } \
        if (representation == 'r') \
        { \
            read_result = this->read_base("rvorticity", this->rvorticity); \
            if (read_result != EXIT_SUCCESS) \
                return read_result; \
            else \
                FFTW(execute)(*((FFTW(plan)*)this->r2c_vorticity )); \
        } \
        this->force_divfree(this->cvorticity); \
        this->symmetrize(this->cvorticity, 3); \
        return EXIT_SUCCESS; \
    } \
    if ((field == 'u') && (representation == 'c')) \
        return this->read_base("cvelocity", this->cvelocity); \
    if ((field == 'u') && (representation == 'r')) \
        return this->read_base("rvelocity", this->rvelocity); \
    return EXIT_FAILURE; \
} \
 \
template<> \
int fluid_solver<R>::write(char field, char representation) \
{ \
    char fname[512]; \
    if ((field == 'v') && (representation == 'c')) \
    { \
        this->fill_up_filename("cvorticity", fname); \
        return this->cd->write(fname, (void*)this->cvorticity); \
    } \
    if ((field == 'v') && (representation == 'r')) \
    { \
        FFTW(execute)(*((FFTW(plan)*)this->c2r_vorticity )); \
        clip_zero_padding<R>(this->rd, this->rvorticity, 3); \
        this->fill_up_filename("rvorticity", fname); \
        return this->rd->write(fname, this->rvorticity); \
    } \
    this->compute_velocity(this->cvorticity); \
    if ((field == 'u') && (representation == 'c')) \
    { \
        this->fill_up_filename("cvelocity", fname); \
        return this->cd->write(fname, this->cvelocity); \
    } \
    if ((field == 'u') && (representation == 'r')) \
    { \
        this->ift_velocity(); \
        clip_zero_padding<R>(this->rd, this->rvelocity, 3); \
        this->fill_up_filename("rvelocity", fname); \
        return this->rd->write(fname, this->rvelocity); \
    } \
    return EXIT_FAILURE; \
} \
 \
template<> \
void fluid_solver<R>::compute_Eulerian_acceleration(R *acceleration) \
{ \
    this->omega_nonlin(0); \
    CLOOP_K2( \
            for (int cc=0; cc<3; cc++) \
                for (int i=0; i<2; i++) \
                    this->cv[1][3*cindex+cc][i] = this->cu[3*cindex+cc][i] - this->nu*k2*this->cv[0][3*cindex+cc][i]; \
            ); \
    CLOOP_K2( \
            if (k2 <= this->kM2 && k2 > 0) \
            { \
                this->cv[2][3*cindex+0][0] = -(this->ky[yindex]*this->cv[1][3*cindex+2][1] - this->kz[zindex]*this->cv[1][3*cindex+1][1]) / k2; \
                this->cv[2][3*cindex+1][0] = -(this->kz[zindex]*this->cv[1][3*cindex+0][1] - this->kx[xindex]*this->cv[1][3*cindex+2][1]) / k2; \
                this->cv[2][3*cindex+2][0] = -(this->kx[xindex]*this->cv[1][3*cindex+1][1] - this->ky[yindex]*this->cv[1][3*cindex+0][1]) / k2; \
                this->cv[2][3*cindex+0][1] =  (this->ky[yindex]*this->cv[1][3*cindex+2][0] - this->kz[zindex]*this->cv[1][3*cindex+1][0]) / k2; \
                this->cv[2][3*cindex+1][1] =  (this->kz[zindex]*this->cv[1][3*cindex+0][0] - this->kx[xindex]*this->cv[1][3*cindex+2][0]) / k2; \
                this->cv[2][3*cindex+2][1] =  (this->kx[xindex]*this->cv[1][3*cindex+1][0] - this->ky[yindex]*this->cv[1][3*cindex+0][0]) / k2; \
            } \
            ); \
    if (this->cd->myrank == this->cd->rank[0]) \
        std::fill_n((R*)(this->cv[2]), 6, 0.0); \
    FFTW(execute)(*((FFTW(plan)*)this->vc2r[2])); \
    std::copy( \
            this->rv[2], \
            this->rv[2] + 2*this->cd->local_size, \
            acceleration); \
} \
 \
template<> \
void fluid_solver<R>::compute_pressure(FFTW(complex) *pressure) \
{ \
    /* assume velocity is already in real space representation */ \
    ptrdiff_t tindex; \
    \
    /* diagonal terms 11 22 33 */\
    RLOOP ( \
            this, \
            tindex = 3*rindex; \
            for (int cc=0; cc<3; cc++) \
                this->rv[1][tindex+cc] = this->ru[tindex+cc]*this->ru[tindex+cc]; \
            ); \
    this->clean_up_real_space(this->rv[1], 3); \
    FFTW(execute)(*((FFTW(plan)*)this->vr2c[1])); \
    this->dealias(this->cv[1], 3); \
    CLOOP_K2( \
            if (k2 <= this->kM2 && k2 > 0) \
            { \
                tindex = 3*cindex; \
                for (int i=0; i<2; i++) \
                { \
                    pressure[cindex][i] = -(this->kx[xindex]*this->kx[xindex]*this->cv[1][tindex+0][i] + \
                                            this->ky[yindex]*this->ky[yindex]*this->cv[1][tindex+1][i] + \
                                            this->kz[zindex]*this->kz[zindex]*this->cv[1][tindex+2][i]); \
                } \
            } \
            else \
                std::fill_n((R*)(pressure+cindex), 2, 0.0); \
            ); \
    /* off-diagonal terms 12 23 31 */\
    RLOOP ( \
            this, \
            tindex = 3*rindex; \
            for (int cc=0; cc<3; cc++) \
                this->rv[1][tindex+cc] = this->ru[tindex+cc]*this->ru[tindex+(cc+1)%3]; \
            ); \
    this->clean_up_real_space(this->rv[1], 3); \
    FFTW(execute)(*((FFTW(plan)*)this->vr2c[1])); \
    this->dealias(this->cv[1], 3); \
    CLOOP_K2( \
            if (k2 <= this->kM2 && k2 > 0) \
            { \
                tindex = 3*cindex; \
                for (int i=0; i<2; i++) \
                { \
                    pressure[cindex][i] -= 2*(this->kx[xindex]*this->ky[yindex]*this->cv[1][tindex+0][i] + \
                                              this->ky[yindex]*this->kz[zindex]*this->cv[1][tindex+1][i] + \
                                              this->kz[zindex]*this->kx[xindex]*this->cv[1][tindex+2][i]); \
                    pressure[cindex][i] /= this->normalization_factor*k2; \
                } \
            } \
            ); \
} \
 \
template<> \
void fluid_solver<R>::compute_gradient_statistics( \
        FFTW(complex) *vec, \
        double *moments, \
        ptrdiff_t *hist, \
        ptrdiff_t *QR2D_hist, \
        double max_estimates[], \
        int nbins, \
        int QR2D_nbins) \
{ \
    FFTW(complex) *ca; \
    R *ra; \
    ca = FFTW(alloc_complex)(this->cd->local_size*3); \
    ra = (R*)(ca); \
    this->compute_vector_gradient(ca, vec); \
    for (int cc=0; cc<3; cc++) \
    { \
        std::copy( \
                (R*)(ca + cc*this->cd->local_size), \
                (R*)(ca + (cc+1)*this->cd->local_size), \
                (R*)this->cv[1]); \
        FFTW(execute)(*((FFTW(plan)*)this->vc2r[1])); \
        std::copy( \
                this->rv[1], \
                this->rv[1] + this->cd->local_size*2, \
                ra + cc*this->cd->local_size*2); \
    } \
    /* velocity gradient is now stored, in real space, in ra */ \
    std::fill_n(this->rv[1], 2*this->cd->local_size, 0.0); \
    R *dx_u, *dy_u, *dz_u; \
    dx_u = ra; \
    dy_u = ra + 2*this->cd->local_size; \
    dz_u = ra + 4*this->cd->local_size; \
    double binsize[2]; \
    double tmp_max_estimate[4]; \
    tmp_max_estimate[0] = max_estimates[0]; \
    tmp_max_estimate[1] = max_estimates[1]; \
    tmp_max_estimate[2] = max_estimates[2]; \
    tmp_max_estimate[3] = 1.0; /* 4th component is gonna be disregarded anyway... */ \
    binsize[0] = 2*tmp_max_estimate[2] / QR2D_nbins; \
    binsize[1] = 2*tmp_max_estimate[1] / QR2D_nbins; \
    ptrdiff_t *local_hist = new ptrdiff_t[QR2D_nbins*QR2D_nbins]; \
    std::fill_n(local_hist, QR2D_nbins*QR2D_nbins, 0); \
    RLOOP( \
            this, \
            R AxxAxx; \
            R AyyAyy; \
            R AzzAzz; \
            R AxyAyx; \
            R AyzAzy; \
            R AzxAxz; \
            R Sxy; \
            R Syz; \
            R Szx; \
            ptrdiff_t tindex = 3*rindex; \
            AxxAxx = dx_u[tindex+0]*dx_u[tindex+0]; \
            AyyAyy = dy_u[tindex+1]*dy_u[tindex+1]; \
            AzzAzz = dz_u[tindex+2]*dz_u[tindex+2]; \
            AxyAyx = dx_u[tindex+1]*dy_u[tindex+0]; \
            AyzAzy = dy_u[tindex+2]*dz_u[tindex+1]; \
            AzxAxz = dz_u[tindex+0]*dx_u[tindex+2]; \
            this->rv[1][tindex+1] = - ((AxxAxx + AyyAyy + AzzAzz)/2 - AxyAyx - AyzAzy - AzxAxz); \
            this->rv[1][tindex+2] = - (dx_u[tindex+0]*(AxxAxx/3 + AxyAyx + AzxAxz) + \
                                       dy_u[tindex+1]*(AyyAyy/3 + AxyAyx + AyzAzy) + \
                                       dz_u[tindex+2]*(AzzAzz/3 + AzxAxz + AyzAzy) + \
                                       dx_u[tindex+1]*dy_u[tindex+2]*dz_u[tindex+0] + \
                                       dx_u[tindex+2]*dy_u[tindex+0]*dz_u[tindex+1]); \
            int bin0 = int((this->rv[1][tindex+2] + tmp_max_estimate[2]) / binsize[0]); \
            int bin1 = int((this->rv[1][tindex+1] + tmp_max_estimate[1]) / binsize[1]); \
            if ((bin0 >= 0 && bin0 < QR2D_nbins) && \
                (bin1 >= 0 && bin1 < QR2D_nbins)) \
                local_hist[bin1*QR2D_nbins + bin0]++; \
            Sxy = dx_u[tindex+1]+dy_u[tindex+0]; \
            Syz = dy_u[tindex+2]+dz_u[tindex+1]; \
            Szx = dz_u[tindex+0]+dx_u[tindex+2]; \
            this->rv[1][tindex] = (AxxAxx + AyyAyy + AzzAzz + \
                                   (Sxy*Sxy + Syz*Syz + Szx*Szx)/2); \
            ); \
    FFTW(free)(ca); \
    MPI_Allreduce( \
            local_hist, \
            QR2D_hist, \
            QR2D_nbins * QR2D_nbins, \
            MPI_INT64_T, MPI_SUM, this->cd->comm); \
    delete[] local_hist; \
    double *tmp_moments = new double[4*10]; \
    ptrdiff_t *tmp_hist = new ptrdiff_t[4*nbins]; \
    this->compute_rspace_stats( \
            this->rv[1], \
            tmp_moments, \
            tmp_hist, \
            tmp_max_estimate, \
            nbins); \
    for (int i=0; i<10; i++) \
        for (int j=0; j<3; j++) \
            moments[i*3+j] = tmp_moments[i*4+j]; \
    delete[] tmp_moments; \
    for (int i=0; i<nbins; i++) \
        for (int j=0; j<3; j++) \
            hist[i*3+j] = tmp_hist[i*4+j]; \
    delete[] tmp_hist; \
} \
 \
template<> \
void fluid_solver<R>::compute_Lagrangian_acceleration(R *acceleration) \
{ \
    ptrdiff_t tindex; \
    FFTW(complex) *pressure; \
    pressure = FFTW(alloc_complex)(this->cd->local_size/3); \
    this->compute_velocity(this->cvorticity); \
    this->ift_velocity(); \
    this->compute_pressure(pressure); \
    this->compute_velocity(this->cvorticity); \
    std::fill_n((R*)this->cv[1], 2*this->cd->local_size, 0.0); \
    CLOOP_K2( \
            if (k2 <= this->kM2) \
            { \
                tindex = 3*cindex; \
                for (int cc=0; cc<3; cc++) \
                    for (int i=0; i<2; i++) \
                        this->cv[1][tindex+cc][i] = - this->nu*k2*this->cu[tindex+cc][i]; \
                if (strcmp(this->forcing_type, "linear") == 0) \
                { \
                    double knorm = sqrt(k2); \
                    if ((this->fk0 <= knorm) && \
                        (this->fk1 >= knorm)) \
                        for (int c=0; c<3; c++) \
                            for (int i=0; i<2; i++) \
                                this->cv[1][tindex+c][i] += this->famplitude*this->cu[tindex+c][i]; \
                } \
                this->cv[1][tindex+0][0] += this->kx[xindex]*pressure[cindex][1]; \
                this->cv[1][tindex+1][0] += this->ky[yindex]*pressure[cindex][1]; \
                this->cv[1][tindex+2][0] += this->kz[zindex]*pressure[cindex][1]; \
                this->cv[1][tindex+0][1] -= this->kx[xindex]*pressure[cindex][0]; \
                this->cv[1][tindex+1][1] -= this->ky[yindex]*pressure[cindex][0]; \
                this->cv[1][tindex+2][1] -= this->kz[zindex]*pressure[cindex][0]; \
            } \
            ); \
    FFTW(execute)(*((FFTW(plan)*)this->vc2r[1])); \
    std::copy( \
            this->rv[1], \
            this->rv[1] + 2*this->cd->local_size, \
            acceleration); \
    FFTW(free)(pressure); \
} \

/*****************************************************************************/



/*****************************************************************************/
/* now actually use the macro defined above                                  */
FLUID_SOLVER_DEFINITIONS(
        FFTW_MANGLE_FLOAT,
        float,
        MPI_FLOAT,
        MPI_COMPLEX)
FLUID_SOLVER_DEFINITIONS(
        FFTW_MANGLE_DOUBLE,
        double,
        MPI_DOUBLE,
        BFPS_MPICXX_DOUBLE_COMPLEX)
/*****************************************************************************/



/*****************************************************************************/
/* finally, force generation of code for single precision                    */
template class fluid_solver<float>;
template class fluid_solver<double>;
/*****************************************************************************/