From 7c6c1877dc096629aae9156f56ed05057f95e776 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Wed, 9 Dec 2015 17:56:12 +0100
Subject: [PATCH] partial fix for rFFTW part+interp

---
 bfps/cpp/interpolator.cpp       |  6 ++++++
 bfps/cpp/rFFTW_interpolator.cpp | 10 +++++++++-
 bfps/cpp/rFFTW_particles.cpp    | 12 +++---------
 3 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/bfps/cpp/interpolator.cpp b/bfps/cpp/interpolator.cpp
index 28312c31..cf04af66 100644
--- a/bfps/cpp/interpolator.cpp
+++ b/bfps/cpp/interpolator.cpp
@@ -163,9 +163,11 @@ void interpolator<rnumber, interp_neighbours>::operator()(
     }
     std::fill_n(dest, 3, 0);
     ptrdiff_t bigiz, bigiy, bigix;
+    double tval[3];
     for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
     {
         bigiz = ptrdiff_t(xg[2]+iz);
+        std::fill_n(tval, 3, 0);
         for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
         {
             bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
@@ -180,9 +182,13 @@ void interpolator<rnumber, interp_neighbours>::operator()(
                     dest[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
                                                                               by[iy+interp_neighbours]*
                                                                               bx[ix+interp_neighbours]);
+                    tval[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
+                                                                              by[iy+interp_neighbours]*
+                                                                              bx[ix+interp_neighbours]);
                 }
             }
         }
+        DEBUG_MSG("%ld %d %d %g %g %g\n", bigiz, xg[1], xg[0], tval[0], tval[1], tval[2]);
     }
 }
 
diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index 5e5f00ee..e4bd10a8 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -103,10 +103,13 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
     }
     std::fill_n(dest, 3, 0);
     ptrdiff_t bigiz, bigiy, bigix;
+    double tval[3];
     for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
     {
         bigiz = ptrdiff_t(((xg[2]+iz) + this->descriptor->sizes[0]) % this->descriptor->sizes[0]);
         if (this->descriptor->myrank == this->descriptor->rank[bigiz])
+        {
+            std::fill_n(tval, 3, 0);
             for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
             {
                 bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
@@ -116,14 +119,19 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
                     for (int c=0; c<3; c++)
                     {
                         ptrdiff_t tindex = ((bigiz *this->descriptor->sizes[1] +
-                                             bigiy)*this->descriptor->sizes[2] +
+                                             bigiy)*(this->descriptor->sizes[2]+2) +
                                              bigix)*3+c;
                         dest[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
                                                                                   by[iy+interp_neighbours]*
                                                                                   bx[ix+interp_neighbours]);
+                        tval[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
+                                                                                  by[iy+interp_neighbours]*
+                                                                                  bx[ix+interp_neighbours]);
                     }
                 }
             }
+            DEBUG_MSG("%ld %d %d %g %g %g\n", bigiz, xg[1], xg[0], tval[0], tval[1], tval[2]);
+        }
     }
 }
 
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
index 401aa72c..af225f95 100644
--- a/bfps/cpp/rFFTW_particles.cpp
+++ b/bfps/cpp/rFFTW_particles.cpp
@@ -141,13 +141,7 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::sample_vec_fiel
     this->get_grid_coordinates(x, xg, xx);
     /* perform interpolation */
     for (int p=0; p<this->nparticles; p++)
-    {
-        if (this->myrank == this->get_rank(x[3*p+2]))
-        {
-            DEBUG_MSG("particle %d\n", p);
-            (*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
-        }
-    }
+        (*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
     MPI_Allreduce(
             yy,
             y,
@@ -292,10 +286,10 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_grid_coordi
             xg[p*3+c] = MOD(int(tval), this->fs->rd->sizes[2-c]);
             xx[p*3+c] = (x[p*this->ncomponents+c] - tval*grid_size[c]) / grid_size[c];
         }
-        xg[p*3+2] -= this->fs->rd->starts[0];
+        /*xg[p*3+2] -= this->fs->rd->starts[0];
         if (this->myrank == this->fs->rd->rank[0] &&
             xg[p*3+2] > this->fs->rd->subsizes[0])
-            xg[p*3+2] -= this->fs->rd->sizes[0];
+            xg[p*3+2] -= this->fs->rd->sizes[0];*/
     }
 }
 
-- 
GitLab