add broken version of new class

the code runs through the first iteration, so now I just need to fix the
redistribution of particles. oh, and the I/O...

bfps/cpp/distributed_particles.cpp

@@ -56,6 +56,9 @@ distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_pa
     this->vel = FIELD;
     this->rhs.resize(INTEGRATION_STEPS);
     this->integration_steps = INTEGRATION_STEPS;
+    this->state.reserve(2*this->nparticles / this->nprocs);
+    for (unsigned int i=0; i<this->rhs.size(); i++)
+        this->rhs[i].reserve(2*this->nparticles / this->nprocs);
 }
 
 template <int particle_type, class rnumber, int interp_neighbours>
@@ -66,11 +69,12 @@ distributed_particles<particle_type, rnumber, interp_neighbours>::~distributed_p
 template <int particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
         interpolator<rnumber, interp_neighbours> *field,
+        const std::unordered_map<int, single_particle_state<particle_type>> &x,
         std::unordered_map<int, single_particle_state<POINT3D>> &y)
 {
     double *yy = new double[3];
     y.clear();
-    for (auto &pp: this->state)
+    for (auto &pp: x)
     {
         (*field)(pp.second.data, yy);
         y[pp.first] = yy;
@@ -83,14 +87,16 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::get_rhs(
         const std::unordered_map<int, single_particle_state<particle_type>> &x,
         std::unordered_map<int, single_particle_state<particle_type>> &y)
 {
-    double *yy = new double[this->ncomponents];
-    y.clear();
-    for (auto &pp: x)
+    std::unordered_map<int, single_particle_state<POINT3D>> yy;
+    switch(particle_type)
     {
-        (*this->vel)(pp.second.data, yy);
-        y[pp.first] = yy;
+        case VELOCITY_TRACER:
+            this->sample(this->vel, this->state, yy);
+            y.clear();
+            for (auto &pp: x)
+                y[pp.first] = yy[pp.first].data;
+            break;
     }
-    delete[] yy;
 }
 
 template <int particle_type, class rnumber, int interp_neighbours>
@@ -99,7 +105,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
         const char *dset_name)
 {
     std::unordered_map<int, single_particle_state<POINT3D>> y;
-    this->sample(field, y);
+    this->sample(field, this->state, y);
     this->write(dset_name, y);
 }
 

bfps/cpp/distributed_particles.hpp

@@ -74,6 +74,7 @@ class distributed_particles: public particles_io_base<particle_type>
                 const char *dset_name);
         void sample(
                 interpolator<rnumber, interp_neighbours> *field,
+                const std::unordered_map<int, single_particle_state<particle_type>> &x,
                 std::unordered_map<int, single_particle_state<POINT3D>> &y);
         void get_rhs(
                 const std::unordered_map<int, single_particle_state<particle_type>> &x,

bfps/cpp/interpolator.hpp

@@ -42,6 +42,7 @@ class interpolator:public interpolator_base<rnumber, interp_neighbours>
         rnumber *field;
 
     public:
+        using interpolator_base<rnumber, interp_neighbours>::operator();
         ptrdiff_t buffer_size;
 
         /* descriptor for buffered field */
@@ -67,17 +68,6 @@ class interpolator:public interpolator_base<rnumber, interp_neighbours>
                 const double *__restrict__ x,
                 double *__restrict__ y,
                 const int *deriv = NULL);
-        /* interpolate 1 point */
-        inline void operator()(
-                const double *__restrict__ x,
-                double *__restrict__ dest,
-                const int *deriv = NULL)
-        {
-            int xg[3];
-            double xx[3];
-            this->get_grid_coordinates(x, xg, xx);
-            (*this)(xg, xx, dest, deriv);
-        }
         void operator()(
                 const int *__restrict__ xg,
                 const double *__restrict__ xx,

bfps/cpp/interpolator_base.hpp

@@ -87,6 +87,18 @@ class interpolator_base
                 const double *__restrict__ xx,
                 double *__restrict__ dest,
                 const int *deriv = NULL) = 0;
+
+        /* interpolate 1 point */
+        inline void operator()(
+                const double *__restrict__ x,
+                double *__restrict__ dest,
+                const int *deriv = NULL)
+        {
+            int xg[3];
+            double xx[3];
+            this->get_grid_coordinates(x, xg, xx);
+            (*this)(xg, xx, dest, deriv);
+        }
 };
 
 #endif//INTERPOLATOR_BASE

bfps/cpp/rFFTW_distributed_particles.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 <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <unordered_map>
+#include <vector>
+#include <hdf5.h>
+#include "base.hpp"
+#include "particles_base.hpp"
+#include "fluid_solver_base.hpp"
+#include "rFFTW_interpolator.hpp"
+
+#ifndef RFFTW_DISTRIBUTED_PARTICLES
+
+#define RFFTW_DISTRIBUTED_PARTICLES
+
+template <int particle_type, class rnumber, int interp_neighbours>
+class rFFTW_distributed_particles: public particles_io_base<particle_type>
+{
+    private:
+        std::unordered_map<int, single_particle_state<particle_type> > state;
+        std::vector<std::unordered_map<int, single_particle_state<particle_type>>> rhs;
+        std::vector<MPI_Comm> interp_comm;
+        std::vector<int> interp_nprocs;
+
+    public:
+        int integration_steps;
+        // this class only works with rFFTW interpolator
+        rFFTW_interpolator<rnumber, interp_neighbours> *vel;
+
+        /* simulation parameters */
+        double dt;
+
+        /* methods */
+
+        /* constructor and destructor.
+         * allocate and deallocate:
+         *  this->state
+         *  this->rhs
+         * */
+        rFFTW_distributed_particles(
+                const char *NAME,
+                const hid_t data_file_id,
+                rFFTW_interpolator<rnumber, interp_neighbours> *FIELD,
+                const int TRAJ_SKIP,
+                const int INTEGRATION_STEPS = 2);
+        ~rFFTW_distributed_particles();
+
+        void sample(
+                rFFTW_interpolator<rnumber, interp_neighbours> *field,
+                const char *dset_name);
+        void sample(
+                rFFTW_interpolator<rnumber, interp_neighbours> *field,
+                const std::unordered_map<int, single_particle_state<particle_type>> &x,
+                std::unordered_map<int, single_particle_state<POINT3D>> &y);
+        void get_rhs(
+                const std::unordered_map<int, single_particle_state<particle_type>> &x,
+                std::unordered_map<int, single_particle_state<particle_type>> &y);
+
+        void redistribute(
+                std::unordered_map<int, single_particle_state<particle_type>> &x,
+                std::vector<std::unordered_map<int, single_particle_state<particle_type>>> &vals);
+
+
+        /* input/output */
+        void read();
+        void write(
+                const char *dset_name,
+                std::unordered_map<int, single_particle_state<POINT3D>> &y);
+        void write(
+                const char *dset_name,
+                std::unordered_map<int, single_particle_state<particle_type>> &y);
+        void write(const bool write_rhs = true);
+
+        /* solvers */
+        void step();
+        void roll_rhs();
+        void AdamsBashforth(const int nsteps);
+};
+
+#endif//RFFTW_DISTRIBUTED_PARTICLES
+

bfps/cpp/rFFTW_interpolator.cpp

@@ -54,6 +54,24 @@ rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
     delete[] this->compute;
 }
 
+template <class rnumber, int interp_neighbours>
+bool rFFTW_interpolator<rnumber, interp_neighbours>::get_rank_info(double z, int &maxz_rank, int &minz_rank)
+{
+    int zg = int(floor(z/this->dz));
+    minz_rank = this->descriptor->rank[MOD(
+             zg - interp_neighbours,
+            this->descriptor->sizes[0])];
+    maxz_rank = this->descriptor->rank[MOD(
+            zg + 1 + interp_neighbours,
+            this->descriptor->sizes[0])];
+    bool is_here = false;
+    for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
+        is_here = (is_here ||
+                   (this->descriptor->myrank ==
+                    this->descriptor->rank[MOD(zg+iz, this->descriptor->sizes[0])]));
+    return is_here;
+}
+
 template <class rnumber, int interp_neighbours>
 void rFFTW_interpolator<rnumber, interp_neighbours>::sample(
         const int nparticles,

bfps/cpp/rFFTW_interpolator.hpp

@@ -37,6 +37,7 @@ template <class rnumber, int interp_neighbours>
 class rFFTW_interpolator:public interpolator_base<rnumber, interp_neighbours>
 {
     public:
+        using interpolator_base<rnumber, interp_neighbours>::operator();
         /* size of field to interpolate */
         ptrdiff_t field_size;
 
@@ -62,6 +63,8 @@ class rFFTW_interpolator:public interpolator_base<rnumber, interp_neighbours>
             return EXIT_SUCCESS;
         }
 
+        bool get_rank_info(double z, int &maxz_rank, int &minz_rank);
+
         /* interpolate field at an array of locations */
         void sample(
                 const int nparticles,

bfps/tools.py

@@ -232,13 +232,16 @@ def particle_finite_diff_test(
     if interp_name not in pars.keys():
         # old format
         interp_name = 'tracers{0}_field'.format(species)
+    interp_name = pars[interp_name].value
+    if type(interp_name) == bytes:
+        interp_name = str(interp_name, 'ASCII')
     to_print = (
             'steps={0}, interp={1}, neighbours={2}, '.format(
                 pars['tracers{0}_integration_steps'.format(species)].value,
-                pars[str(pars[interp_name].value, 'ASCII') + '_type'].value,
-                pars[str(pars[interp_name].value, 'ASCII') + '_neighbours'].value))
-    if 'spline' in str(pars[interp_name].value, 'ASCII'):
-        to_print += 'smoothness = {0}, '.format(pars[str(pars[interp_name].value, 'ASCII') + '_smoothness'].value)
+                pars[interp_name + '_type'].value,
+                pars[interp_name + '_neighbours'].value))
+    if 'spline' in interp_name:
+        to_print += 'smoothness = {0}, '.format(pars[interp_name + '_smoothness'].value)
     to_print += (
             'SNR d1p-vel={0:.3f}'.format(np.mean(snr_vel1)))
     if acc_on:

setup.py

@@ -92,6 +92,7 @@ print('This is bfps version ' + VERSION)
 ### lists of files and MANIFEST.in
 src_file_list = ['field_descriptor',
                  'interpolator_base',
+                 'rFFTW_distributed_particles',
                  'distributed_particles',
                  'particles_base',
                  'interpolator',

tests/test_particle_classes.py

@@ -47,7 +47,7 @@ def scaling(opt):
     c1.compute_statistics()
     compare_stats(opt, c0, c1)
     opt.work_dir = wd + '/N{0:0>3x}_3'.format(opt.n)
-    c2 = launch(opt, dt = c0.parameters['dt'], particle_class = 'particles', interpolator_class = 'interpolator')
+    c2 = launch(opt, dt = c0.parameters['dt'], particle_class = 'rFFTW_distributed_particles')
     c2.compute_statistics()
     compare_stats(opt, c0, c2)
     compare_stats(opt, c1, c2)