rm compilation warnings,separate communicators

CMakeLists.txt

@@ -32,7 +32,10 @@ elseif (${CMAKE_SYSTEM_PROCESSOR} STREQUAL "x86_64") ## Xeon
      set(CMAKE_CXX_COMPILER "icpc")  
    endif()
    #set(CMAKE_CXX_FLAGS "-O3 -openmp -g -xHost -fno-exceptions -vec-report")
-   set(CMAKE_CXX_FLAGS  "-O3 -h noomp -DNO_HDF5  ") # -DPRINTPCL for printing particle subcycles  -DNO_HDF5 to disable output
+   # -DPRINTPCL for printing particle subcycles 
+   #-DNO_HDF5 to disable output 
+   #-h noomp if running on bluewaters
+   set(CMAKE_CXX_FLAGS  "-O3 -DNO_HDF5 ")
 else()
    set(CMAKE_CXX_FLAGS "-O3")
 endif()

communication/Com3DNonblk.cpp

@@ -23,15 +23,17 @@
 //isCenterFlag: 1 communicateCenter; 0 communicateNode
 void NBDerivedHaloComm(int nx, int ny, int nz, double ***vector,const VirtualTopology3D * vct, EMfields3D *EMf,bool isCenterFlag, bool isFaceOnlyFlag, bool needInterp, bool isParticle)
 {
-	MPI_Errhandler_set(MPI_COMM_WORLD,MPI_ERRORS_RETURN);
+    const MPI_Comm comm       = isParticle ?vct->getParticleComm()      :vct->getFieldComm();
+#ifdef DEBUG
+	MPI_Errhandler_set(comm,MPI_ERRORS_RETURN);
+#endif
     MPI_Status  stat[12];
     MPI_Request reqList[12];				  //at most 6 requests x 2 (send recv)
     int communicationCnt[6] = {0,0,0,0,0,0};  //1 if there communication on that dir
     int recvcnt = 0,sendcnt = 0;              //request counter
 
     const int tag_XL=1,tag_YL=2,tag_ZL=3,tag_XR=4,tag_YR=5,tag_ZR=6;//To address same rank as left and right neighbour in periodic case
     const int myrank    	  = vct->getCartesian_rank();
-    const MPI_Comm comm       = isParticle ?vct->getParticleComm()      :vct->getFieldComm();
     const int right_neighborX = isParticle ?vct->getXright_neighbor_P() :vct->getXright_neighbor();
     const int left_neighborX  = isParticle ?vct->getXleft_neighbor_P()  :vct->getXleft_neighbor();
     const int right_neighborY = isParticle ?vct->getYright_neighbor_P() :vct->getYright_neighbor();
@@ -134,12 +136,14 @@ void NBDerivedHaloComm(int nx, int ny, int nz, double ***vector,const VirtualTop
             int error_code = stat[si].MPI_ERROR;
             if (error_code != MPI_SUCCESS) {
                 stopFlag = true;
+#ifdef DEBUG
                 char error_string[100];
                 int length_of_error_string, error_class;
 
                 MPI_Error_class(error_code, &error_class);
                 MPI_Error_string(error_class, error_string, &length_of_error_string);
                 dprintf("MPI_Waitall error at %d  %s\n",si, error_string);
+#endif
             }
         }
         if(stopFlag) exit (EXIT_FAILURE);
@@ -360,12 +364,14 @@ void NBDerivedHaloComm(int nx, int ny, int nz, double ***vector,const VirtualTop
 				int error_code = stat[si].MPI_ERROR;
 				if (error_code != MPI_SUCCESS) {
 					stopFlag = true;
+#ifdef DEBUG
 					char error_string[100];
 					int length_of_error_string, error_class;
 
 					MPI_Error_class(error_code, &error_class);
 					MPI_Error_string(error_class, error_string, &length_of_error_string);
 					dprintf("MPI_Waitall error at %d  %s\n",si, error_string);
+#endif
 				}
 			}
 			if(stopFlag) exit (EXIT_FAILURE);;
@@ -467,12 +473,14 @@ void NBDerivedHaloComm(int nx, int ny, int nz, double ***vector,const VirtualTop
 				int error_code = stat[si].MPI_ERROR;
 				if (error_code != MPI_SUCCESS) {
 					stopFlag = true;
+#ifdef DEBUG
 					char error_string[100];
 					int length_of_error_string, error_class;
 
 					MPI_Error_class(error_code, &error_class);
 					MPI_Error_string(error_class, error_string, &length_of_error_string);
 					dprintf("MPI_Waitall error at %d  %s\n",si, error_string);
+#endif
 				}
 			}
 			if(stopFlag) exit (EXIT_FAILURE);

fields/EMfields3D.cpp

@@ -29,13 +29,13 @@
 #include "CG.h"
 #include "GMRES.h"
 #include "Particles3Dcomm.h"
-#include "TimeTasks.h"
 #include "Moments.h"
 #include "Parameters.h"
 #include "ompdefs.h"
 #include "debug.h"
 #include "string.h"
 #include "mic_particles.h"
+#include "TimeTasks.h"
 #include "ipicmath.h" // for roundup_to_multiple
 #include "Alloc.h"
 #include "asserts.h"

include/Basic.h

@@ -50,7 +50,7 @@ developers: Stefano Markidis, Giovanni Lapenta
 
 
 /** method to calculate the parallel dot product with vect1, vect2 having the ghost cells*/
-double dotP(const double *vect1, const double *vect2, int n);
+double dotP(const double *vect1, const double *vect2, int n,MPI_Comm* comm);
 /** method to calculate dot product */
 double dot(const double *vect1, const double *vect2, int n);
 /** method to calculate the square norm of a vector */
@@ -60,11 +60,11 @@ double norm2(const arr3_double vect, int nx, int ny);
 /** method to calculate the square norm of a vector */
 double norm2(const double *vect, int nx);
 /** method to calculate the parallel dot product */
-double norm2P(const arr3_double vect, int nx, int ny, int nz);
+//double norm2P(const arr3_double vect, int nx, int ny, int nz);
 /** method to calculate the parallel norm of a vector on different processors with the ghost cell */
-double norm2P(const double *vect, int n);
+//double norm2P(const double *vect, int n);
 /** method to calculate the parallel norm of a vector on different processors with the gost cell*/
-double normP(const double *vect, int n);
+double normP(const double *vect, int n,MPI_Comm* comm);
 /** method to calculate the difference of two vectors*/
 void sub(double *res, const double *vect1, const double *vect2, int n);
 /** method to calculate the sum of two vectors vector1 = vector1 + vector2*/

include/BlockCommunicator.h

@@ -159,7 +159,7 @@ struct Block
     //  listID, commID,
     //  dest.rank(), dest.tag_name(),
     //  nop);
-    MPI_Isend(&block[0], NUMBERS_PER_ELEMENT*block.size(), MPI_DOUBLE,
+    MPI_Isend((double*)&block[0], NUMBERS_PER_ELEMENT*block.size(), MPI_DOUBLE,
       dest.rank(), dest.tag(), dest.comm(), &request);
     //dprintf("finished sending block number %d", listID);
   }
@@ -188,7 +188,7 @@ struct Block
     // make sure that space exists to receive
     int newsize = signal_hack() ? capacity+1 : capacity;
     block.resize(newsize);
-    MPI_Irecv(&block[0], NUMBERS_PER_ELEMENT*block.size(), MPI_DOUBLE,
+    MPI_Irecv((double*)&block[0], NUMBERS_PER_ELEMENT*block.size(), MPI_DOUBLE,
       source.rank(), source.tag(), source.comm(), &request);
   }
   // processing received data

include/Com3DNonblk.h

@@ -35,11 +35,9 @@ developers           : Stefano Markidis, Ivy Bo Peng
 #include "mpi.h"
 #include "BcFields3D.h"
 #include "VCtopology3D.h"
-#include "TimeTasks.h"
 #include "ipicdefs.h"
 #include "Alloc.h"
 #include "debug.h"
-//#include "parallel.h"
 #include "EMfields3D.h"
 
 /** communicate ghost cells (FOR NODES) */

include/Connection.h

@@ -41,7 +41,7 @@
 // one way is to compare the positions of the centroids, first in
 // the x direction, then in the y direction, and then in the z
 // direction.  One might use tag=1 for downward communication and
-// tag=2 for upward communication, and use MPI_COMM_WORLD for the
+// tag=2 for upward communication, and all ways use particle communicator for the
 // group.
 //
 namespace Direction
@@ -83,14 +83,14 @@ class Connection
   Connection():
     _rank(0),
     _tag(0),
-    _comm(MPI_COMM_WORLD)
+    _comm(MPI_COMM_NULL)
   {}
-  Connection(int rank_, int tag_, MPI_Comm comm_ = MPI_COMM_WORLD)
+  Connection(int rank_, int tag_, MPI_Comm comm_ = MPI_COMM_NULL)
   {
     init(rank_,tag_,comm_);
   }
  private:
-  void init(int rank_, int tag_, MPI_Comm comm_ = MPI_COMM_WORLD)
+  void init(int rank_, int tag_, MPI_Comm comm_ = MPI_COMM_NULL)
   {
     _rank = rank_;
     _tag = tag_;

include/MPIdata.h

@@ -29,7 +29,9 @@ email                : [email protected], [email protected]
 
 #ifndef MPIDATA_H
 #define MPIDATA_H
-
+#ifndef NO_MPI
+#include <mpi.h>
+#endif
 /**
  * MPI Data Structure. This class contains:
  *
@@ -69,7 +71,10 @@ class MPIdata {
 public:
   static int get_rank(){return instance().rank;}
   static int get_nprocs(){return instance().nprocs;}
+  static int get_PicGlobalComm(){return instance().PIC_COMM;}
 private:
+    /*iPIC3D Global Communicator*/
+  static MPI_Comm PIC_COMM;
   /** rank of the process */
   static int rank;
   /** number of processes */

inputoutput/ParallelIO.cpp

@@ -301,7 +301,7 @@ void ReadFieldsH5hut(int nspec, EMfields3D *EMf, Collective *col, VCtopology3D *
   infile.CloseFieldsFile();
 
   // initialize B on centers
-  MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(MPIdata::get_PicGlobalComm());
 
   // Comm ghost nodes for B-field
   communicateNodeBC(grid->getNXN(), grid->getNYN(), grid->getNZN(), EMf->getBx(), col->bcBx[0],col->bcBx[1],col->bcBx[2],col->bcBx[3],col->bcBx[4],col->bcBx[5], vct);

main/iPIC3Dlib.cpp

@@ -111,7 +111,7 @@ int c_Solver::Init(int argc, char **argv) {
     }
   }
   // We create a new communicator with a 3D virtual Cartesian topology
-  vct->setup_vctopology(MPI_COMM_WORLD);
+    vct->setup_vctopology(MPIdata::get_PicGlobalComm());
   {
     stringstream num_proc_ss;
     num_proc_ss << vct->getCartesian_rank();

particles/Particles3D.cpp

@@ -1920,7 +1920,7 @@ void Particles3D::openbc_particles_inflow()
     for (int j = 1; j < grid->getNYC() - 1; j++)
       for (int k = 1; k < grid->getNZC() - 1; k++){
 
-    	/*
+    	
     	//The below is uniformly distributed in location
 		for (int ii=0; ii < npcelx; ii++)
 		  for (int jj=0; jj < npcely; jj++)
@@ -2281,7 +2281,7 @@ double Particles3D::delta_f(double u, double v, double w, double x, double y, do
     a3[l + lmax] += df0_dvpar(vpar, vperp) * bessel_Jn_array[l];
   }
 
-  deltaf = (0.0, 0.0);
+  //deltaf = (0.0, 0.0);
   for (register int l = -lmax; l < lmax + 1; l++) {
     deltaf += (a3[l + lmax] * Ex_mod * exp(I * Ex_phase) + a1[l + lmax] * Ey_mod * exp(I * Ey_phase) + a2[l + lmax] * Ez_mod * exp(I * Ez_phase)) / (kpar * vpar + l * om_c - omega) * exp(-I * phi * (double) l);
   }

particles/Particles3Dcomm.cpp

@@ -33,7 +33,6 @@ developers: Stefano Markidis, Giovanni Lapenta.
 #include "VCtopology3D.h"
 #include "Collective.h"
 #include "Alloc.h"
-#include "Basic.h"
 #include "Grid3DCU.h"
 #include "Field.h"
 #include "MPIdata.h"

performances/Timing.cpp

@@ -22,7 +22,7 @@
 #include "stdio.h"
 #include "Timing.h"
 #include "ipicdefs.h"
-
+#include "MPIdata.h"
 /**
  * 
  * series of methods for timing and profiling
@@ -56,7 +56,7 @@ Timing::Timing(int my_rank) {
   // MPE_Describe_state(event2a,event2b,"Field","blue"); // the mover is blue in the visualizer
   // MPE_Describe_state(event3a,event3b,"Interp P->G","yellow"); // the interpolation particle->Grid is yellow in the visualizer
   // }
-  former_MPI_Barrier(MPI_COMM_WORLD);
+    former_MPI_Barrier(MPIdata::get_PicGlobalComm());
   // start the log
   // MPE_Start_log();
 
@@ -65,12 +65,12 @@ Timing::Timing(int my_rank) {
 /** start the timer */
 void Timing::startTiming() {
   ttick = MPI_Wtick();
-  former_MPI_Barrier(MPI_COMM_WORLD);
+  former_MPI_Barrier(MPIdata::get_PicGlobalComm());
   tstart = MPI_Wtime();
 }
 /** stop the timer */
 void Timing::stopTiming() {
-  former_MPI_Barrier(MPI_COMM_WORLD);
+  former_MPI_Barrier(MPIdata::get_PicGlobalComm());
   tend = MPI_Wtime();
   texecution = tend - tstart;
   if (rank_id == 0) {

solvers/CG.cpp

@@ -23,6 +23,9 @@
 #include "CG.h"
 #include "Basic.h"
 #include "parallel.h"
+//#include "ipicdefs.h"
+#include "EMfields3D.h"
+#include "VCtopology3D.h"
 
 /**
  * 
@@ -44,6 +47,8 @@ bool CG(double *xkrylov, int xkrylovlen, double *b, int maxit, double tol, FIELD
   eqValue(0.0, v, xkrylovlen);
   eqValue(0.0, z, xkrylovlen);
   eqValue(0.0, im, xkrylovlen);
+
+  MPI_Comm fieldcomm = (field->get_vct()).getFieldComm();
 
   int i = 0;
   bool CONVERGED = false;
@@ -55,7 +60,7 @@ bool CG(double *xkrylov, int xkrylovlen, double *b, int maxit, double tol, FIELD
   sub(r, b, im, xkrylovlen);
   // v = r
   eq(v, r, xkrylovlen);
-  c = dotP(r, r, xkrylovlen);
+  c = dotP(r, r, xkrylovlen,&fieldcomm);
   initial_error = sqrt(c);
   if (is_output_thread())
     printf("CG Initial error: %g\n", initial_error);
@@ -64,12 +69,12 @@ bool CG(double *xkrylov, int xkrylovlen, double *b, int maxit, double tol, FIELD
     return (true);
   while (i < maxit) {
     (field->*FunctionImage) (z, v);
-    t = c / dotP(v, z, xkrylovlen);
+    t = c / dotP(v, z, xkrylovlen,&fieldcomm);
     // x(i+1) = x + t*v
     addscale(t, xkrylov, v, xkrylovlen);
     // r(i+1) = r - t*z
     addscale(-t, r, z, xkrylovlen);
-    d = dotP(r, r, xkrylovlen);
+    d = dotP(r, r, xkrylovlen,&fieldcomm);
     if (CGVERBOSE && is_output_thread())
     {
       printf("Iteration # %d - norm of residual relative to initial error %g\n",

solvers/GMRES.cpp

@@ -26,7 +26,9 @@
 #include "errors.h"
 #include "Alloc.h"
 #include "TimeTasks.h"
-#include "ipicdefs.h"
+//#include "ipicdefs.h"
+#include "EMfields3D.h"
+#include "VCtopology3D.h"
 
 void GMRES(FIELD_IMAGE FunctionImage, double *xkrylov, int xkrylovlen,
   const double *b, int m, int max_iter, double tol, Field * field)
@@ -73,7 +75,9 @@ void GMRES(FIELD_IMAGE FunctionImage, double *xkrylov, int xkrylovlen,
             "------------------------------------\n\n");
   }
 
-  double normb = normP(b, xkrylovlen);
+  MPI_Comm fieldcomm = (field->get_vct()).getFieldComm();
+
+  double normb = normP(b, xkrylovlen,&fieldcomm);
   if (normb == 0.0)
     normb = 1.0;
 
@@ -83,7 +87,7 @@ void GMRES(FIELD_IMAGE FunctionImage, double *xkrylov, int xkrylovlen,
     // r = b - A*x
     (field->*FunctionImage) (im, xkrylov);
     sub(r, b, im, xkrylovlen);
-    initial_error = normP(r, xkrylovlen);
+    initial_error = normP(r, xkrylovlen,&fieldcomm);
 
     if (itr == 0) {
       if (is_output_thread())
@@ -126,15 +130,15 @@ void GMRES(FIELD_IMAGE FunctionImage, double *xkrylov, int xkrylovlen,
       {
 
         MPI_Allreduce(MPI_IN_PLACE, y, (k+2),
-          MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+          MPI_DOUBLE, MPI_SUM, fieldcomm);
       }
       for (int j = 0; j <= k; j++) {
         H[j][k] = y[j];
         addscale(-H[j][k], V[k+1], V[j], xkrylovlen);
       }
       // Is there a numerically stable way to
       // eliminate this second all-reduce all?
-      H[k+1][k] = normP(V[k+1], xkrylovlen);
+      H[k+1][k] = normP(V[k+1], xkrylovlen,&fieldcomm);
       //
       // check that vectors are orthogonal
       //
@@ -150,11 +154,11 @@ void GMRES(FIELD_IMAGE FunctionImage, double *xkrylov, int xkrylovlen,
       if (av + delta * H[k + 1][k] == av)
       {
         for (int j = 0; j <= k; j++) {
-          const double htmp = dotP(w, V[j], xkrylovlen);
+          const double htmp = dotP(w, V[j], xkrylovlen,&fieldcomm);
           H[j][k] = H[j][k] + htmp;
           addscale(-htmp, w, V[j], xkrylovlen);
         }
-        H[k + 1][k] = normP(w, xkrylovlen);
+        H[k + 1][k] = normP(w, xkrylovlen,&fieldcomm);
       }
       // normalize the new vector
       scale(w, (1.0 / H[k + 1][k]), xkrylovlen);