issue CmPA#52: use separate fieldForPcls array to push particles

communication/ComNodes3D.cpp

@@ -1,4 +1,5 @@
 
+#include "mpi.h"
 #include "ComNodes3D.h"
 #include "TimeTasks.h"
 #include "ipicdefs.h"

fields/EMfields3D.cpp

@@ -52,6 +52,7 @@ EMfields3D::EMfields3D(Collective * col, Grid * grid) :
   //
   // array allocation: nodes
   //
+  fieldForPcls  (nxn, nyn, nzn, 6),
   Ex   (nxn, nyn, nzn),
   Ey   (nxn, nyn, nzn),
   Ez   (nxn, nyn, nzn),
@@ -226,9 +227,8 @@ void EMfields3D::sumMoments(const Particles3Dcomm& pcls, Grid * grid, VirtualTop
   double const*const q = pcls.getQall();
   //
   const int is = pcls.get_ns();
-  bool bmoments10 = true;
 
-  // if b10moments
+  #ifdef TENMOMENTS
   double* rhons1d = &rhons[is][0][0][0];
   double* Jxs1d   = &Jxs  [is][0][0][0];
   double* Jys1d   = &Jys  [is][0][0][0];
@@ -239,6 +239,7 @@ void EMfields3D::sumMoments(const Particles3Dcomm& pcls, Grid * grid, VirtualTop
   double* pYYsn1d = &pYYsn[is][0][0][0];
   double* pYZsn1d = &pYZsn[is][0][0][0];
   double* pZZsn1d = &pZZsn[is][0][0][0];
+  #endif
   //
   const long long nop_ll = pcls.getNOP();
   const int nop = pcls.getNOP();
@@ -325,7 +326,7 @@ void EMfields3D::sumMoments(const Particles3Dcomm& pcls, Grid * grid, VirtualTop
       const double weight110 = qi * xi1 * eta1 * zeta0 * invVOL;
       const double weight111 = qi * xi1 * eta1 * zeta1 * invVOL;
 
-      if(bmoments10)
+      // add particle to moments
       {
         arr1_double_fetch moments000 = moments[ix  ][iy  ][iz  ];
         arr1_double_fetch moments001 = moments[ix  ][iy  ][iz-1];
@@ -1375,6 +1376,27 @@ void EMfields3D::ConstantChargePlanet(Grid * grid, VirtualTopology3D * vct, doub
 
 }
 
+/*! Populate the field data used to push particles */
+// 
+// One could add a background magnetic field B_ext at this point,
+// which was incompletely implemented in commit 05082fc8ad688
+//
+void EMfields3D::set_fieldForPcls()
+{
+  #pragma omp parallel for collapse(3)
+  for(int i=0;i<nxn;i++)
+  for(int j=0;j<nyn;j++)
+  for(int k=0;k<nzn;k++)
+  {
+    fieldForPcls[i][j][k][0] = (pfloat) Bxn[i][j][k];
+    fieldForPcls[i][j][k][1] = (pfloat) Byn[i][j][k];
+    fieldForPcls[i][j][k][2] = (pfloat) Bzn[i][j][k];
+    fieldForPcls[i][j][k][3] = (pfloat) Ex[i][j][k];
+    fieldForPcls[i][j][k][4] = (pfloat) Ey[i][j][k];
+    fieldForPcls[i][j][k][5] = (pfloat) Ez[i][j][k];
+  }
+}
+
 /*! Calculate Magnetic field with the implicit solver: calculate B defined on nodes With E(n+ theta) computed, the magnetic field is evaluated from Faraday's law */
 void EMfields3D::calculateB(Grid * grid, VirtualTopology3D * vct, Collective *col) {
   if (vct->getCartesian_rank() == 0)

grids/Grid3DCU.cpp

@@ -50,12 +50,18 @@ Grid3DCU::Grid3DCU(CollectiveIO * col, VirtualTopology3D * vct) {
   zEnd = zStart + (col->getLz() / (double) vct->getZLEN());
 
   // arrays allocation: nodes ---> the first node has index 1, the last has index nxn-2!
+  pfloat_node_xcoord = new pfloat[nxn];
+  pfloat_node_ycoord = new pfloat[nyn];
+  pfloat_node_zcoord = new pfloat[nzn];
   node_xcoord = new double[nxn];
   node_ycoord = new double[nyn];
   node_zcoord = new double[nzn];
   for (int i=0; i<nxn; i++) node_xcoord[i] = xStart + (i - 1) * dx;
   for (int j=0; j<nyn; j++) node_ycoord[j] = yStart + (j - 1) * dy;
   for (int k=0; k<nzn; k++) node_zcoord[k] = zStart + (k - 1) * dz;
+  for (int i=0; i<nxn; i++) pfloat_node_xcoord[i] = node_xcoord[i];
+  for (int j=0; j<nyn; j++) pfloat_node_ycoord[j] = node_ycoord[j];
+  for (int k=0; k<nzn; k++) pfloat_node_zcoord[k] = node_zcoord[k];
   // arrays allocation: cells ---> the first cell has index 1, the last has index ncn-2!
   center_xcoord = new double[nxc];
   center_ycoord = new double[nyc];

include/EMfields3D.h

@@ -118,6 +118,8 @@ class EMfields3D                // :public Field
     /*! smooth the electric field */
     void smoothE(double value, VirtualTopology3D * vct, Collective *col);
 
+    /*! copy the field data to the array used to move the particles */
+    void set_fieldForPcls();
     /*! communicate ghost for grid -> Particles interpolation */
     void communicateGhostP2G(int ns, int bcFaceXright, int bcFaceXleft, int bcFaceYright, int bcFaceYleft, VirtualTopology3D * vct);
     void sumMoments(const Particles3Dcomm& pcls, Grid * grid, VirtualTopology3D * vct);
@@ -187,6 +189,7 @@ class EMfields3D                // :public Field
     double getBy(int X, int Y, int Z) const { return Byn.get(X,Y,Z);}
     double getBz(int X, int Y, int Z) const { return Bzn.get(X,Y,Z);}
     //
+    const_arr4_pfloat get_fieldForPcls() { return fieldForPcls; }
     arr3_double getEx() { return Ex; }
     arr3_double getEy() { return Ey; }
     arr3_double getEz() { return Ez; }
@@ -339,6 +342,11 @@ class EMfields3D                // :public Field
     /*! PHI: electric potential (indexX, indexY, indexZ), defined on central points between nodes */
     array3_double PHI;
 
+    // Electric field component used to move particles
+    // organized for rapid access in mover_PC()
+    // [This is the information transferred from cluster to booster].
+    array4_pfloat fieldForPcls;
+
     // Electric field components defined on nodes
     //
     array3_double Ex;

include/Grid3DCU.h

@@ -7,8 +7,6 @@
 #ifndef GRID3DCU_H
 #define GRID3DCU_H
 
-#include <iostream>
-
 #include "Grid.h"
 #include "CollectiveIO.h"
 #include "ComInterpNodes3D.h"
@@ -142,6 +140,9 @@ class Grid3DCU                  // :public Grid
   /** invol = inverse of volume*/
   double invVOL;
   /** node coordinate */
+  pfloat *pfloat_node_xcoord;
+  pfloat *pfloat_node_ycoord;
+  pfloat *pfloat_node_zcoord;
   double *node_xcoord;
   double *node_ycoord;
   double *node_zcoord;
@@ -169,6 +170,9 @@ class Grid3DCU                  // :public Grid
   //const double &calcXN(int X) { return xStart+(X-1)*dx;}
   //const double &calcYN(int Y) { return yStart+(Y-1)*dy;}
   //const double &calcZN(int Z) { return zStart+(Z-1)*dz;}
+  const pfloat &get_pfloat_XN(int X) { return pfloat_node_xcoord[X];}
+  const pfloat &get_pfloat_YN(int Y) { return pfloat_node_ycoord[Y];}
+  const pfloat &get_pfloat_ZN(int Z) { return pfloat_node_zcoord[Z];}
   const double &getXN(int X) { return node_xcoord[X];}
   const double &getYN(int Y) { return node_ycoord[Y];}
   const double &getZN(int Z) { return node_zcoord[Z];}

include/Particles.h

@@ -21,10 +21,6 @@ developers: Stefano Markidis, Giovanni Lapenta
  *
  */
 
-// precision to use for particles
-//typedef float pfloat;
-typedef double pfloat;
-
 class Particles {
 public:
   /** allocate particles */

include/arraysfwd.h

@@ -1,6 +1,7 @@
 /* forward declaration for array classes */
 #ifndef arraysfwd_h
 #define arraysfwd_h
+#include "ipicdefs.h" // for pfloat
 
 namespace iPic3D
 {
@@ -41,6 +42,7 @@ namespace iPic3D
 //
 typedef iPic3D::const_array_ref3<double> const_arr3_double;
 typedef iPic3D::const_array_ref4<double> const_arr4_double;
+typedef iPic3D::const_array_ref4<pfloat> const_arr4_pfloat;
 typedef iPic3D::array_ref1<double> arr1_double;
 typedef iPic3D::array_ref2<double> arr2_double;
 typedef iPic3D::array_ref3<double> arr3_double;
@@ -49,13 +51,16 @@ typedef iPic3D::array1<double> array1_double;
 typedef iPic3D::array2<double> array2_double;
 typedef iPic3D::array3<double> array3_double;
 typedef iPic3D::array4<double> array4_double;
+typedef iPic3D::array4<pfloat> array4_pfloat;
 // This directive should be consistent with the directives in Alloc.h
 #if defined(FLAT_ARRAYS) || defined(CHECK_BOUNDS)
 typedef iPic3D::array_fetch1<double> arr1_double_fetch;
 typedef iPic3D::array_get1<double> arr1_double_get;
+typedef iPic3D::array_get1<pfloat> arr1_pfloat_get;
 #else
 typedef double* arr1_double_fetch;
 typedef double* arr1_double_get;
+typedef pfloat* arr1_pfloat_get;
 #endif
 
 #endif

include/ipicdefs.h

@@ -10,4 +10,23 @@
 // use precprocessor to remove MPI_Barrier() calls.
 #define MPI_Barrier(args...)
 
+//#define SINGLE_PRECISION_PCLS
+//
+// single precision does not seem to help on the MIC
+#ifdef SINGLE_PRECISION_PCLS
+  typedef float pfloat;
+  #ifdef __MIC__
+    #define VECTOR_WIDTH 16
+  #else
+    #define VECTOR_WIDTH 8
+  #endif
+#else
+  #ifdef __MIC__
+    #define VECTOR_WIDTH 8
+  #else
+    #define VECTOR_WIDTH 4
+  #endif
+  typedef double pfloat;
+#endif
+
 #endif

main/iPic3Dlib.cpp

@@ -218,9 +218,13 @@ bool c_Solver::ParticlesMover() {
   /*  -------------- */
 
   timeTasks.start(TimeTasks::PARTICLES);
+  // Should change this to add background guide field
+  EMf->set_fieldForPcls();
   for (int i = 0; i < ns; i++)  // move each species
   {
     // #pragma omp task inout(part[i]) in(grid) target_device(booster)
+    //
+    // should merely pass EMf->get_fieldForPcls() rather than EMf.
     mem_avail = part[i].mover_PC(grid, vct, EMf); // use the Predictor Corrector scheme 
   }
   timeTasks.end(TimeTasks::PARTICLES);