move linear interpolation to slab_field_particles

src/slab_field_particles.cpp

@@ -20,7 +20,7 @@
 
 
 
-#define NDEBUG
+//#define NDEBUG
 
 #include <cmath>
 #include <cassert>
@@ -276,6 +276,20 @@ void slab_field_particles<rnumber>::get_grid_coordinates(double *x, int *xg, dou
     }
 }
 
+template <class rnumber>
+void slab_field_particles<rnumber>::linear_interpolation(float *field, int *xg, double *xx, double *dest)
+{
+    for (int c=0; c<3; c++)
+        dest[c] = (field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(1-xx[1])*(1-xx[2])) +
+                   field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(1-xx[1])*(1-xx[2])) +
+                   field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(  xx[1])*(1-xx[2])) +
+                   field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(  xx[1])*(1-xx[2])) +
+                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(1-xx[1])*(  xx[2])) +
+                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(1-xx[1])*(  xx[2])) +
+                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(  xx[1])*(  xx[2])) +
+                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(  xx[1])*(  xx[2])));
+}
+
 template <class rnumber>
 void slab_field_particles<rnumber>::read()
 {

src/slab_field_particles.hpp

@@ -95,6 +95,7 @@ class slab_field_particles
         void rFFTW_to_buffered(rnumber *src, rnumber *dst);
         ptrdiff_t buffered_local_size();
         void get_grid_coordinates(double *x, int *xg, double *xx);
+        void linear_interpolation(float *field, int *xg, double *xx, double *dest);
         void read();
         void write();
 

src/tracers.cpp

@@ -27,20 +27,6 @@
 #include "fftw_tools.hpp"
 #include "tracers.hpp"
 
-template <class rnumber>
-void tracers<rnumber>::linear_interpolation(float *field, int *xg, double *xx, double *dest)
-{
-    for (int c=0; c<3; c++)
-        dest[c] = (field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(1-xx[1])*(1-xx[2])) +
-                   field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(1-xx[1])*(1-xx[2])) +
-                   field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(  xx[1])*(1-xx[2])) +
-                   field[((ptrdiff_t(xg[2]  )*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(  xx[1])*(1-xx[2])) +
-                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(1-xx[1])*(  xx[2])) +
-                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]  )*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(1-xx[1])*(  xx[2])) +
-                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]  )*3+c]*((1-xx[0])*(  xx[1])*(  xx[2])) +
-                   field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[2]+xg[0]+1)*3+c]*((  xx[0])*(  xx[1])*(  xx[2])));
-}
-
 template <class rnumber>
 void tracers<rnumber>::jump_estimate(double *jump)
 {

src/tracers.hpp

@@ -47,7 +47,6 @@ class tracers:public slab_field_particles<rnumber>
         void update_field(bool clip_on = true);
         virtual void get_rhs(double *x, double *rhs);
         virtual void jump_estimate(double *jump_length);
-        void linear_interpolation(float *field, int *xg, double *xx, double *dest);
 };
 
 

test.py

@@ -47,9 +47,9 @@ def main(opt):
     c.parameters['dt'] = 2e-3
     c.parameters['niter_todo'] = opt.nsteps
     c.parameters['famplitude'] = 0.0
-    c.parameters['nparticles'] = 32
+    c.parameters['nparticles'] = 1
     if opt.run:
-        c.execute(ncpu = opt.ncpu)
+        c.execute(ncpu = opt.ncpu, particle_rseed = 1)
     dtype = pickle.load(open(c.name + '_dtype.pickle'))
     stats1 = np.fromfile('test1_stats.bin', dtype = dtype)
     stats2 = np.fromfile('test2_stats.bin', dtype = dtype)