Merge branch 'feature/new-solver' into develop

Add a new solver, "NSVorticityEquation".
Same algorithm as "NavierStokes", but it uses the field and kspace
classes throughout, and it can't handle particles.
The purpose is to have a bare-bone fluid solver that can take advantage
of optimizations and new backends created for the field class.

bfps/NavierStokes.py

@@ -196,6 +196,9 @@ class NavierStokes(_fluid_particle_base):
         self.stat_src += """
                 //begincpp
                 *tmp_vec_field = fs->cvorticity;
+                DEBUG_MSG("%g %g\\n",
+                    tmp_vec_field->cval(tmp_vec_field->get_cindex(2, 1, 1), 0, 0),
+                    tmp_vec_field->cval(tmp_vec_field->get_cindex(2, 1, 1), 0, 1));
                 switch(fs->dealias_type)
                 {
                     case 0:

bfps/__init__.py

@@ -49,4 +49,5 @@ from host_information import host_info
 from .FluidConvert import FluidConvert
 from .FluidResize import FluidResize
 from .NavierStokes import NavierStokes
+from .NSVorticityEquation import NSVorticityEquation
 

bfps/__main__.py

@@ -29,6 +29,7 @@ import argparse
 
 import bfps
 from .NavierStokes import NavierStokes
+from .NSVorticityEquation import NSVorticityEquation
 from .FluidResize import FluidResize
 from .FluidConvert import FluidConvert
 from .NSManyParticles import NSManyParticles
@@ -45,6 +46,12 @@ def main():
                  'NS',
                  'NS-single',
                  'NS-double']
+    NSVEoptions = ['NSVorticityEquation',
+                 'NSVorticityEquation-single',
+                 'NSVorticityEquation-double',
+                 'NSVE',
+                 'NSVE-single',
+                 'NSVE-double']
     FRoptions = ['FluidResize',
                  'FluidResize-single',
                  'FluidResize-double',
@@ -57,7 +64,7 @@ def main():
                'NSManyParticles-double']
     parser.add_argument(
             'base_class',
-            choices = NSoptions + FRoptions + FCoptions + NSMPopt,
+            choices = NSoptions + NSVEoptions + FRoptions + FCoptions + NSMPopt,
             type = str)
     # first option is the choice of base class or -h or -v
     # all other options are passed on to the base_class instance
@@ -70,6 +77,8 @@ def main():
         precision = 'single'
     if opt.base_class in NSoptions:
         base_class = NavierStokes
+    if opt.base_class in NSVEoptions:
+        base_class = NSVorticityEquation
     elif opt.base_class in FRoptions:
         base_class = FluidResize
     elif opt.base_class in FCoptions:

bfps/_base.py

@@ -83,8 +83,7 @@ class _base(object):
                    'hid_t dset, memtype, space;\n' +
                    'char fname[256];\n' +
                    'hsize_t dims[1];\n' +
-                   'char *string_data;\n' +
-                   #'char string_data[512];\n' +
+                   'char string_data[512];\n' +
                    'sprintf(fname, "%s.h5", simname);\n' +
                    'parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);\n')
         for i in range(len(key)):
@@ -95,14 +94,8 @@ class _base(object):
             elif type(parameters[key[i]]) == str:
                 src_txt += ('space = H5Dget_space(dset);\n' +
                             'memtype = H5Dget_type(dset);\n' +
-                            #'H5Sget_simple_extent_dims(space, dims, NULL);\n' +
-                            #'DEBUG_MSG("dims[0] = %ld\\n", dims[0]);\n' +
-                            'string_data = (char*)malloc(512*sizeof(char));\n' +
                             'H5Dread(dset, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &string_data);\n' +
                             'sprintf({0}, "%s", string_data);\n'.format(key[i]) +
-                            #'DEBUG_MSG("{0} = %s\\n", {0});\n'.format(key[i]) +
-                            #'delete[] string_data;\n' +
-                            'free(string_data);\n'
                             'H5Sclose(space);\n' +
                             'H5Tclose(memtype);\n')
             elif type(parameters[key[i]]) == np.ndarray:

bfps/cpp/distributed_particles.cpp

@@ -45,17 +45,17 @@ template <particle_types particle_type, class rnumber, int interp_neighbours>
 distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_particles(
         const char *NAME,
         const hid_t data_file_id,
-        interpolator<rnumber, interp_neighbours> *FIELD,
+        interpolator<rnumber, interp_neighbours> *VEL,
         const int TRAJ_SKIP,
         const int INTEGRATION_STEPS) : particles_io_base<particle_type>(
             NAME,
             TRAJ_SKIP,
             data_file_id,
-            FIELD->descriptor->comm)
+            VEL->descriptor->comm)
 {
     assert((INTEGRATION_STEPS <= 6) &&
            (INTEGRATION_STEPS >= 1));
-    this->vel = FIELD;
+    this->vel = VEL;
     this->rhs.resize(INTEGRATION_STEPS);
     this->integration_steps = INTEGRATION_STEPS;
     this->state.reserve(2*this->nparticles / this->nprocs);

bfps/cpp/field.hpp

@@ -24,110 +24,16 @@
 
 
 
-#include <mpi.h>
 #include <hdf5.h>
-#include <fftw3-mpi.h>
 #include <unordered_map>
 #include <vector>
 #include <string>
-#include "base.hpp"
+#include "kspace.hpp"
 
-#ifndef FIELD
+#ifndef FIELD_HPP
 
-#define FIELD
+#define FIELD_HPP
 
-enum field_backend {FFTW};
-enum field_components {ONE, THREE, THREExTHREE};
-enum kspace_dealias_type {TWO_THIRDS, SMOOTH};
-
-constexpr unsigned int ncomp(
-        field_components fc)
-    /* return actual number of field components for each enum value */
-{
-    return ((fc == THREE) ? 3 : (
-            (fc == THREExTHREE) ? 9 : 1));
-}
-
-constexpr unsigned int ndim(
-        field_components fc)
-    /* return actual number of field dimensions for each enum value */
-{
-    return ((fc == THREE) ? 4 : (
-            (fc == THREExTHREE) ? 5 : 3));
-}
-
-template <field_components fc>
-class field_layout
-{
-    public:
-        /* description */
-        hsize_t sizes[ndim(fc)];
-        hsize_t subsizes[ndim(fc)];
-        hsize_t starts[ndim(fc)];
-        hsize_t local_size, full_size;
-
-        int myrank, nprocs;
-        MPI_Comm comm;
-
-        std::vector<std::vector<int>> rank;
-        std::vector<std::vector<int>> all_start;
-        std::vector<std::vector<int>> all_size;
-
-        /* methods */
-        field_layout(
-                const hsize_t *SIZES,
-                const hsize_t *SUBSIZES,
-                const hsize_t *STARTS,
-                const MPI_Comm COMM_TO_USE);
-        ~field_layout(){}
-};
-
-template <field_backend be,
-          kspace_dealias_type dt>
-class kspace
-{
-    public:
-        /* relevant field layout */
-        field_layout<ONE> *layout;
-
-        /* physical parameters */
-        double dkx, dky, dkz, dk, dk2;
-
-        /* mode and dealiasing information */
-        double kMx, kMy, kMz, kM, kM2;
-        double kMspec, kMspec2;
-        std::vector<double> kx, ky, kz;
-        std::unordered_map<int, double> dealias_filter;
-        std::vector<double> kshell;
-        std::vector<int64_t> nshell;
-        int nshells;
-
-        /* methods */
-        template <field_components fc>
-        kspace(
-                const field_layout<fc> *source_layout,
-                const double DKX = 1.0,
-                const double DKY = 1.0,
-                const double DKZ = 1.0);
-        ~kspace();
-
-        template <typename rnumber,
-                  field_components fc>
-        void low_pass(rnumber *__restrict__ a, const double kmax);
-
-        template <typename rnumber,
-                  field_components fc>
-        void dealias(rnumber *__restrict__ a);
-
-        template <typename rnumber,
-                  field_components fc>
-        void cospectrum(
-                const rnumber(* __restrict__ a)[2],
-                const rnumber(* __restrict__ b)[2],
-                const hid_t group,
-                const std::string dset_name,
-                const hsize_t toffset);
-};
 
 template <typename rnumber,
           field_backend be,
@@ -136,10 +42,9 @@ class field
 {
     private:
         /* data arrays */
-        rnumber *data;
-        typedef rnumber cnumber[2];
-        hsize_t npoints;
+        rnumber *__restrict__ data;
     public:
+        hsize_t npoints;
         bool real_space_representation;
         /* basic MPI information */
         int myrank, nprocs;
@@ -153,8 +58,8 @@ class field
         field_layout<fc> *clayout, *rlayout, *rmemlayout;
 
         /* FFT plans */
-        void *c2r_plan;
-        void *r2c_plan;
+        typename fftw_interface<rnumber>::plan c2r_plan;
+        typename fftw_interface<rnumber>::plan r2c_plan;
         unsigned fftw_plan_rigor;
 
         /* HDF5 data types for arrays */
@@ -171,14 +76,22 @@ class field
 
         int io(
                 const std::string fname,
-                const std::string dset_name,
+                const std::string field_name,
+                const int iteration,
+                const bool read = true);
+        int io_database(
+                const std::string fname,
+                const std::string field_name,
                 const int toffset,
                 const bool read = true);
 
+        /* essential FFT stuff */
         void dft();
         void ift();
         void normalize();
+        void symmetrize();
 
+        /* stats */
         void compute_rspace_xincrement_stats(
                 const int xcells,
                 const hid_t group,
@@ -191,16 +104,57 @@ class field
                 const std::string dset_name,
                 const hsize_t toffset,
                 const std::vector<double> max_estimate);
-        inline rnumber *get_rdata()
+
+        /* acess data */
+        inline rnumber *__restrict__ get_rdata()
         {
             return this->data;
         }
-        inline cnumber *get_cdata()
+
+        inline typename fftw_interface<rnumber>::complex *__restrict__ get_cdata()
+        {
+            return (typename fftw_interface<rnumber>::complex*__restrict__)this->data;
+        }
+
+        inline rnumber &rval(ptrdiff_t rindex, int component = 0)
+        {
+            assert(fc == ONE || fc == THREE);
+            assert(component >= 0 && component <ncomp(fc));
+            return *(this->data + rindex*ncomp(fc) + component);
+        }
+
+        inline rnumber &rval(ptrdiff_t rindex, int comp1, int comp0)
+        {
+            assert(fc == THREExTHREE);
+            assert(comp1 >= 0 && comp1 < 3);
+            assert(comp0 >= 0 && comp0 < 3);
+            return *(this->data + ((rindex*3 + comp1)*3 + comp0));
+        }
+
+        inline rnumber &cval(ptrdiff_t cindex, int imag)
+        {
+            assert(fc == ONE);
+            assert(imag == 0 || imag == 1);
+            return *(this->data + cindex*2 + imag);
+        }
+
+        inline rnumber &cval(ptrdiff_t cindex, int component, int imag)
+        {
+            assert(fc == THREE);
+            assert(imag == 0 || imag == 1);
+            return *(this->data + (cindex*ncomp(fc) + component)*2 + imag);
+        }
+
+        inline rnumber &cval(ptrdiff_t cindex, int comp1, int comp0, int imag)
         {
-            return (cnumber*)this->data;
+            assert(fc == THREExTHREE);
+            assert(comp1 >= 0 && comp1 < 3);
+            assert(comp0 >= 0 && comp0 < 3);
+            assert(imag == 0 || imag == 1);
+            return *(this->data + ((cindex*3 + comp1)*3+comp0)*2 + imag);
         }
 
-        inline field<rnumber, be, fc>& operator=(const cnumber *__restrict__ source)
+        inline field<rnumber, be, fc>& operator=(const typename fftw_interface<rnumber>::complex *__restrict__ source)
         {
             std::copy((rnumber*)source,
                       (rnumber*)(source + this->clayout->local_size),
@@ -217,6 +171,15 @@ class field
             this->real_space_representation = true;
             return *this;
         }
+
+        inline field<rnumber, be, fc>& operator=(const rnumber value)
+        {
+            std::fill_n(this->data,
+                        this->rmemlayout->local_size,
+                        value);
+            return *this;
+        }
+
         template <kspace_dealias_type dt>
         void compute_stats(
                 kspace<be, dt> *kk,
@@ -224,6 +187,44 @@ class field
                 const std::string dset_name,
                 const hsize_t toffset,
                 const double max_estimate);
+        inline void impose_zero_mode()
+        {
+            if (this->clayout->myrank == this->clayout->rank[0][0] &&
+                this->real_space_representation == false)
+            {
+                std::fill_n(this->data, 2*ncomp(fc), 0.0);
+            }
+        }
+        template <class func_type>
+        void RLOOP(func_type expression)
+        {
+            switch(be)
+            {
+                case FFTW:
+                    for (hsize_t zindex = 0; zindex < this->rlayout->subsizes[0]; zindex++)
+                    for (hsize_t yindex = 0; yindex < this->rlayout->subsizes[1]; yindex++)
+                    {
+                        ptrdiff_t rindex = (
+                                zindex * this->rlayout->subsizes[1] + yindex)*(
+                                    this->rmemlayout->subsizes[2]);
+                        for (hsize_t xindex = 0; xindex < this->rlayout->subsizes[2]; xindex++)
+                        {
+                            expression(rindex, xindex, yindex, zindex);
+                            rindex++;
+                        }
+                    }
+                    break;
+            }
+        }
+        ptrdiff_t get_cindex(
+                ptrdiff_t xindex,
+                ptrdiff_t yindex,
+                ptrdiff_t zindex)
+        {
+            return ((yindex*this->clayout->subsizes[1] +
+                     zindex)*this->clayout->subsizes[2] +
+                    xindex);
+        }
 };
 
 template <typename rnumber,
@@ -236,73 +237,5 @@ void compute_gradient(
         field<rnumber, be, fc1> *source,
         field<rnumber, be, fc2> *destination);
 
-
-/* real space loop */
-#define FIELD_RLOOP(obj, expression) \
- \
-{ \
-    switch (be) \
-    { \
-        case FFTW: \
-            for (hsize_t zindex = 0; zindex < obj->rlayout->subsizes[0]; zindex++) \
-            for (hsize_t yindex = 0; yindex < obj->rlayout->subsizes[1]; yindex++) \
-            { \
-                ptrdiff_t rindex = ( \
-                        zindex * obj->rlayout->subsizes[1] + yindex)*( \
-                            obj->rmemlayout->subsizes[2]); \
-            for (hsize_t xindex = 0; xindex < obj->rlayout->subsizes[2]; xindex++) \
-                { \
-                    expression; \
-                    rindex++; \
-                } \
-            } \
-            break; \
-    } \
-}
-
-#define KSPACE_CLOOP_K2(obj, expression) \
- \
-{ \
-    double k2; \
-    ptrdiff_t cindex = 0; \
-    for (hsize_t yindex = 0; yindex < obj->layout->subsizes[0]; yindex++) \
-    for (hsize_t zindex = 0; zindex < obj->layout->subsizes[1]; zindex++) \
-    for (hsize_t xindex = 0; xindex < obj->layout->subsizes[2]; xindex++) \
-        { \
-            k2 = (obj->kx[xindex]*obj->kx[xindex] + \
-                  obj->ky[yindex]*obj->ky[yindex] + \
-                  obj->kz[zindex]*obj->kz[zindex]); \
-            expression; \
-            cindex++; \
-        } \
-}
-
-#define KSPACE_CLOOP_K2_NXMODES(obj, expression) \
- \
-{ \
-    double k2; \
-    ptrdiff_t cindex = 0; \
-    for (hsize_t yindex = 0; yindex < obj->layout->subsizes[0]; yindex++) \
-    for (hsize_t zindex = 0; zindex < obj->layout->subsizes[1]; zindex++) \
-    { \
-        int nxmodes = 1; \
-        hsize_t xindex = 0; \
-        k2 = (obj->kx[xindex]*obj->kx[xindex] + \
-              obj->ky[yindex]*obj->ky[yindex] + \
-              obj->kz[zindex]*obj->kz[zindex]); \
-        expression; \
-        cindex++; \
-        nxmodes = 2; \
-    for (xindex = 1; xindex < obj->layout->subsizes[2]; xindex++) \
-        { \
-            k2 = (obj->kx[xindex]*obj->kx[xindex] + \
-                  obj->ky[yindex]*obj->ky[yindex] + \
-                  obj->kz[zindex]*obj->kz[zindex]); \
-            expression; \
-            cindex++; \
-        } \
-    } \
-}
-
-#endif//FIELD
+#endif//FIELD_HPP
 

bfps/cpp/field_layout.cpp

+/**********************************************************************
+*                                                                     *
+*  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                                 *
+*                                                                     *
+**********************************************************************/
+
+
+#include <cassert>
+#include "field_layout.hpp"
+#include "scope_timer.hpp"
+
+template <field_components fc>
+field_layout<fc>::field_layout(
+        const hsize_t *SIZES,
+        const hsize_t *SUBSIZES,
+        const hsize_t *STARTS,
+        const MPI_Comm COMM_TO_USE)
+{
+    TIMEZONE("field_layout::field_layout");
+    this->comm = COMM_TO_USE;
+    MPI_Comm_rank(this->comm, &this->myrank);
+    MPI_Comm_size(this->comm, &this->nprocs);
+
+    std::copy(SIZES, SIZES + 3, this->sizes);
+    std::copy(SUBSIZES, SUBSIZES + 3, this->subsizes);
+    std::copy(STARTS, STARTS + 3, this->starts);
+    if (fc == THREE || fc == THREExTHREE)
+    {
+        this->sizes[3] = 3;
+        this->subsizes[3] = 3;
+        this->starts[3] = 0;
+    }
+    if (fc == THREExTHREE)
+    {
+        this->sizes[4] = 3;
+        this->subsizes[4] = 3;
+        this->starts[4] = 0;
+    }
+    this->local_size = 1;
+    this->full_size = 1;
+    for (unsigned int i=0; i<ndim(fc); i++)
+    {
+        this->local_size *= this->subsizes[i];
+        this->full_size *= this->sizes[i];
+    }
+
+    /*field will at most be distributed in 2D*/
+    this->rank.resize(2);
+    this->all_start.resize(2);
+    this->all_size.resize(2);
+    for (int i=0; i<2; i++)
+    {
+        this->rank[i].resize(this->sizes[i]);
+        std::vector<int> local_rank;
+        local_rank.resize(this->sizes[i], 0);
+        for (unsigned int ii=this->starts[i]; ii<this->starts[i]+this->subsizes[i]; ii++)
+            local_rank[ii] = this->myrank;
+        MPI_Allreduce(
+                &local_rank.front(),
+                &this->rank[i].front(),
+                this->sizes[i],
+                MPI_INT,
+                MPI_SUM,
+                this->comm);
+        this->all_start[i].resize(this->nprocs);
+        std::vector<int> local_start;
+        local_start.resize(this->nprocs, 0);
+        local_start[this->myrank] = this->starts[i];
+        MPI_Allreduce(
+                &local_start.front(),
+                &this->all_start[i].front(),
+                this->nprocs,
+                MPI_INT,
+                MPI_SUM,
+                this->comm);
+        this->all_size[i].resize(this->nprocs);
+        std::vector<int> local_subsize;
+        local_subsize.resize(this->nprocs, 0);
+        local_subsize[this->myrank] = this->subsizes[i];
+        MPI_Allreduce(
+                &local_subsize.front(),
+                &this->all_size[i].front(),
+                this->nprocs,
+                MPI_INT,
+                MPI_SUM,
+                this->comm);
+    }
+}
+
+template class field_layout<ONE>;
+template class field_layout<THREE>;
+template class field_layout<THREExTHREE>;
+

bfps/cpp/field_layout.hpp

+/**********************************************************************
+*                                                                     *
+*  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                                 *
+*                                                                     *
+**********************************************************************/
+
+
+
+#include <vector>
+#include "base.hpp"
+
+#ifndef FIELD_LAYOUT_HPP
+
+#define FIELD_LAYOUT_HPP
+
+enum field_components {ONE, THREE, THREExTHREE};
+
+constexpr unsigned int ncomp(
+        field_components fc)
+    /* return actual number of field components for each enum value */
+{
+    return ((fc == THREE) ? 3 : (
+            (fc == THREExTHREE) ? 9 : 1));
+}
+
+constexpr unsigned int ndim(