Call system() only if needed (#214)

CMakeLists.txt

@@ -174,6 +174,7 @@ include(config/blaslapack.cmake)
 set(SRC 
     src/lib/nekrs.cpp
     src/io/writeFld.cpp
+    src/io/utils.cpp
     src/core/utils/mysort.cpp
     src/core/utils/parallelSort.cpp
     src/core/utils/occaHelpers.cpp

config/nek5000.cmake

@@ -1,6 +1,6 @@
 set(NEK5000_GS_VERSION "1.0.5")
 set(NEK5000_GS_HASH "1b5b28de5b997c3b0893a3b4dcf5cee8614b9f27")
-set(PARRSB_VERSION "0.4")
+set(PARRSB_VERSION "0.6")
 
 if (${NEK5000_PPLIST} MATCHES "PARRSB")
   set(USE_PARRSB on)
@@ -34,6 +34,12 @@ endif()
 
 set(NEK5000_SOURCE_DIR ${nek5000_content_SOURCE_DIR})
 
+if (USE_PARRSB)
+  install(FILES ${NEK5000_SOURCE_DIR}/core/PARALLEL.dprocmap DESTINATION ${NEK5000_SOURCE_DIR}/core RENAME "PARALLEL")
+else()
+  install(FILES ${NEK5000_SOURCE_DIR}/core/PARALLEL.default DESTINATION ${NEK5000_SOURCE_DIR}/core RENAME "PARALLEL")
+endif()
+
 # blasLapack
 # ==========
 
@@ -110,7 +116,7 @@ ExternalProject_Add(
   BUILD_COMMAND 
     ${CMAKE_CURRENT_LIST_DIR}/run_nekconfig.sh 
     "CC=${CMAKE_C_COMPILER}" 
-    "CFLAGS=${EXTERNAL_C_FLAGS}"
+    "NEK5000_SOURCE_DIRCFLAGS=${EXTERNAL_C_FLAGS}"
     "FC=${CMAKE_Fortran_COMPILER}"
     "FFLAGS=${EXTERNAL_Fortran_FLAGS}"
     "NEK5000_SOURCE_DIR=${NEK5000_SOURCE_DIR}"

src/core/setup.cpp

@@ -12,7 +12,7 @@ static occa::memory o_scratch;
 
 static cds_t* cdsSetup(ins_t* ins, mesh_t* mesh, setupAide options, occa::properties &kernelInfoH);
 
-void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, int buildOnly, nrs_t *nrs)
+void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, nrs_t *nrs)
 {
   nrs->options = options;
   nrs->kernelInfo = new occa::properties();
@@ -24,7 +24,9 @@ void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, int buildO
   kernelInfo["include_paths"].asArray();
 
   int N, cubN;
+  int buildOnly = 0;
   string install_dir;
+  if(nrs->options.compareArgs("BUILD ONLY", "TRUE")) buildOnly = 1;
   nrs->options.getArgs("POLYNOMIAL DEGREE", N);
   nrs->options.getArgs("CUBATURE POLYNOMIAL DEGREE", cubN);
   nrs->options.getArgs("NUMBER OF SCALARS", nrs->Nscalar);
@@ -50,10 +52,13 @@ void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, int buildO
     string casename;
     nrs->options.getArgs("CASENAME", casename);
 
+    int err = 0;
     int npTarget = size;
     if (buildOnly) nrs->options.getArgs("NP TARGET", npTarget);
-    if (rank == 0) buildNekInterface(casename.c_str(), mymax(5, nrs->Nscalar), N, npTarget);
-    MPI_Barrier(comm);
+    if (rank == 0) err = buildNekInterface(casename.c_str(), mymax(5, nrs->Nscalar), N, npTarget);
+    MPI_Allreduce(MPI_IN_PLACE, &err, 1, MPI_INT, MPI_SUM, comm);
+    if (err) ABORT(EXIT_FAILURE);; 
+
     if (!buildOnly) {
       nek_setup(comm, nrs->options, nrs);
       nek_setic();
@@ -77,7 +82,7 @@ void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, int buildO
 
   if (nrs->cht && !nrs->options.compareArgs("SCALAR00 IS TEMPERATURE", "TRUE")) {
     if (mesh->rank == 0) cout << "Conjugate heat transfer requires solving for temperature!\n"; 
-    EXIT(1);
+    ABORT(EXIT_FAILURE);;
   } 
 
   { 
@@ -795,7 +800,7 @@ void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, int buildO
       if(nrs->pSolver->levels[0] > mesh->N || 
          nrs->pSolver->levels[nrs->pSolver->nLevels-1] < 1) {
         if(mesh->rank == 0) printf("ERROR: Invalid multigrid coarsening!\n");
-        EXIT(1);
+        ABORT(EXIT_FAILURE);;
       }
       nrs->pOptions.setArgs("MULTIGRID COARSENING","CUSTOM");
     } else if(nrs->pOptions.compareArgs("MULTIGRID DOWNWARD SMOOTHER","ASM") ||

src/core/setup.hpp

@@ -2,6 +2,6 @@
 #define nekrs_inssetup_hpp_
 
 #include "nrs.hpp"
-void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, int buildOnly, nrs_t *nrs);
+void nrsSetup(MPI_Comm comm, occa::device device, setupAide &options, nrs_t *nrs);
 
 #endif

src/core/setupAide.cpp

@@ -68,7 +68,7 @@ string setupAide::readFile(string filename)
   FILE* fh = fopen(filename.c_str(), "r");
   if (fh == 0) {
     printf("Failed to open: %s\n", filename.c_str());
-    exit(1);
+    ABORT(EXIT_FAILURE);
   }
 
   stat(filename.c_str(), &statbuf);

src/core/setupAide.hpp

@@ -37,6 +37,8 @@ SOFTWARE.
 #include <sys/types.h>
 #include <sys/stat.h>
 
+#include "nrssys.hpp"
+
 using std::stringstream;
 using std::fstream;
 using std::string;

src/elliptic/ellipticMultiGridLevelSetup.cpp

@@ -338,12 +338,11 @@ dfloat MGLevel::maxEigSmoothAx()
   const dlong N = Nrows;
   const dlong M = Ncols;
 
-  int k = 10;
-
   hlong Nlocal = (hlong) Nrows;
   hlong Ntotal = 0;
   MPI_Allreduce(&Nlocal, &Ntotal, 1, MPI_HLONG, MPI_SUM, mesh->comm);
-  if(k > Ntotal) k = (int) Ntotal;
+
+  const int k = std::min((hlong) 20, Ntotal); 
 
   // do an arnoldi
 

src/elliptic/ellipticMultiGridSchwarz.cpp

@@ -167,7 +167,7 @@ compute_element_lengths(ElementLengths* lengths, elliptic_t* elliptic)
       std::cout << "Encountered length of zero in middle for element e = " << e << "!\n";
       std::cout << "x,y,z = " << lengths->length_middle_x[e] << ", "
                 << lengths->length_middle_y[e] << ", " << lengths->length_middle_z[e] << "\n";
-      exit(-1);
+      ABORT(EXIT_FAILURE);;
     }
     bool negative = false;
     negative |= lengths->length_middle_x[e] < -tol;
@@ -177,7 +177,7 @@ compute_element_lengths(ElementLengths* lengths, elliptic_t* elliptic)
       std::cout << "Encountered negative length in middle for element e = " << e << "!\n";
       std::cout << "x,y,z = " << lengths->length_middle_x[e] << ", "
                 << lengths->length_middle_y[e] << ", " << lengths->length_middle_z[e] << "\n";
-      exit(-1);
+      ABORT(EXIT_FAILURE);;
     }
   }
 
@@ -506,7 +506,7 @@ void compute_1d_matrices(
     std::cout << "lbc = " << lbc << ", rbc = " << rbc << "\n";
     for(int iface = 0; iface < 6; ++iface)
       std::cout << "EToB[iface] = " << elliptic->EToB[6 * e + iface] << "\n";
-    exit(-1);
+    ABORT(EXIT_FAILURE);;
   }
   if(lbc > 0)
     row_zero(S,nl,0);
@@ -608,6 +608,9 @@ mesh_t* create_extended_mesh(elliptic_t* elliptic)
 {
   mesh_t* meshRoot = elliptic->mesh;
 
+  int buildOnly  = 0;
+  if(elliptic->options.compareArgs("BUILD ONLY", "TRUE")) buildOnly = 1;
+
   mesh_t* mesh = new mesh_t();
   mesh->rank = meshRoot->rank;
   mesh->size = meshRoot->size;
@@ -618,28 +621,31 @@ mesh_t* create_extended_mesh(elliptic_t* elliptic)
   mesh->Nelements = meshRoot->Nelements;
   mesh->Nverts = meshRoot->Nverts;
   mesh->Nfaces = meshRoot->Nfaces;
+  mesh->NfaceVertices = meshRoot->NfaceVertices;
+  mesh->Nnodes = meshRoot->Nnodes;
+
   mesh->EX = (dfloat*) calloc(mesh->Nverts * mesh->Nelements, sizeof(dfloat));
   mesh->EY = (dfloat*) calloc(mesh->Nverts * mesh->Nelements, sizeof(dfloat));
   mesh->EZ = (dfloat*) calloc(mesh->Nverts * mesh->Nelements, sizeof(dfloat));
   memcpy(mesh->EX, meshRoot->EX, mesh->Nverts * mesh->Nelements * sizeof(dfloat));
   memcpy(mesh->EY, meshRoot->EY, mesh->Nverts * mesh->Nelements * sizeof(dfloat));
   memcpy(mesh->EZ, meshRoot->EZ, mesh->Nverts * mesh->Nelements * sizeof(dfloat));
+
+  mesh->faceVertices = (int*) calloc(mesh->NfaceVertices * mesh->Nfaces, sizeof(int));
+  memcpy(mesh->faceVertices, meshRoot->faceVertices, mesh->NfaceVertices * mesh->Nfaces * sizeof(int));
+
   mesh->EToV = (hlong*) calloc(mesh->Nverts * mesh->Nelements, sizeof(hlong));
   memcpy(mesh->EToV, meshRoot->EToV, mesh->Nverts * mesh->Nelements * sizeof(hlong));
 
   meshParallelConnect(mesh);
   meshConnectBoundary(mesh);
-
   meshLoadReferenceNodesHex3D(mesh, mesh->N, 1);
-
-  int buildOnly  = 0;
-  if(elliptic->options.compareArgs("BUILD ONLY", "TRUE")) buildOnly = 1;
-
-  meshPhysicalNodesHex3D(mesh, buildOnly);
-
   meshHaloSetup(mesh);
+  meshPhysicalNodesHex3D(mesh, buildOnly);
+  meshHaloPhysicalNodes(mesh);
   meshConnectFaceNodes3D(mesh);
   meshParallelConnectNodes(mesh, buildOnly);
+
   mesh->ogs = ogsSetup(mesh->Nelements * mesh->Np, mesh->globalIds, mesh->comm, 1, mesh->device);
 
   const int bigNum = 1E9;
@@ -786,8 +792,8 @@ void MGLevel::build(
   elliptic_t* pSolver)
 {
   if(elliptic->elementType != HEXAHEDRA) {
-    printf("ERROR: Unsupported elements type!");
-    exit(-1);
+    printf("ERROR: Unsupported element type!");
+    ABORT(EXIT_FAILURE);
   }
 
   overlap = false;

src/elliptic/ellipticOperator.cpp

@@ -70,7 +70,7 @@ void ellipticAx(elliptic_t* elliptic,
         if(integrationType != 0)
           printf("Precision level (%s) does not support integrationType %d\n", precision, integrationType);
       }
-      exit(1);
+      ABORT(EXIT_FAILURE);
     }
   }
 
@@ -113,7 +113,7 @@ void ellipticAx(elliptic_t* elliptic,
         }
       }
     } else {
-      exit(1);
+      ABORT(EXIT_FAILURE);
     }
     return;
   }

src/elliptic/ellipticPreconditionerSetup.cpp

@@ -37,8 +37,7 @@ void ellipticPreconditionerSetup(elliptic_t* elliptic, ogs_t* ogs, occa::propert
   } else if(options.compareArgs("PRECONDITIONER", "SEMFEM")) {
     //ellipticSEMFEMSetup(elliptic,precon);
     printf("ERROR: SEMFEM does not work right now.\n");
-
-    exit(-1);
+    ABORT(EXIT_FAILURE);;
   } else if(options.compareArgs("PRECONDITIONER", "JACOBI")) {
     dfloat* invDiagA;
     ellipticBuildJacobi(elliptic,&invDiagA);
@@ -48,6 +47,6 @@ void ellipticPreconditionerSetup(elliptic_t* elliptic, ogs_t* ogs, occa::propert
     free(invDiagA);
   } else {
     printf("ERROR: Unknown preconditioner!\n");
-    exit(-1);
+    ABORT(EXIT_FAILURE);
   }
 }

src/elliptic/ellipticSolve.cpp

@@ -64,7 +64,7 @@ int ellipticSolve(elliptic_t* elliptic,
     Niter = pcg (elliptic, o_r, o_x, tol, maxIter);
   }else{
     printf("NONBLOCKING Krylov solvers currently not supported!");
-    exit(1);
+    ABORT(EXIT_FAILURE);
 /*
     if(!options.compareArgs("KRYLOV SOLVER", "FLEXIBLE"))
       Niter = nbpcg (elliptic, o_r, o_x, tol, maxIter);

src/elliptic/ellipticSolveSetup.cpp

@@ -43,17 +43,15 @@ void ellipticSolveSetup(elliptic_t* elliptic, occa::properties kernelInfo)
     if(mesh->rank == 0)
       printf("ERROR: Block solver is implemented for C0-HEXAHEDRA with Jacobi preconditioner only\n");
 
-    MPI_Finalize();
-    exit(-1);
+    ABORT(EXIT_FAILURE);
   }
 
   if (options.compareArgs("COEFFICIENT","VARIABLE") &&  elliptic->elementType != HEXAHEDRA &&
       !options.compareArgs("DISCRETIZATION", "CONTINUOUS")) {
     if(mesh->rank == 0)
       printf("ERROR: Varibale coefficient solver is implemented for C0-HEXAHEDRA only\n");
 
-    MPI_Finalize();
-    exit(-1);
+    ABORT(EXIT_FAILURE);
   }
 
   if (options.compareArgs("COEFFICIENT","VARIABLE")) {
@@ -63,8 +61,8 @@ void ellipticSolveSetup(elliptic_t* elliptic, occa::properties kernelInfo)
       if(mesh->rank == 0)
         printf(
           "ERROR: Varibale coefficient solver is implemented for constant multigrid preconditioner only\n");
-      MPI_Finalize();
-      exit(-1);
+
+      ABORT(EXIT_FAILURE);
     }
   }
 
@@ -828,7 +826,7 @@ void ellipticSolveSetup(elliptic_t* elliptic, occa::properties kernelInfo)
             if(elliptic->var_coeff || elliptic->blockSolver) {
               printf(
                 "ERROR: TRILINEAR form is not implemented for varibale coefficient and block solver yet \n");
-              exit(-1);
+              ABORT(EXIT_FAILURE);
             }
             kernelName = "ellipticPartialAxTrilinear" + suffix;
           }else {

src/io/io.hpp

@@ -1,5 +1,7 @@
 #include "nrs.hpp"
 
+void copyFile(const char *srcName, const char* destName);
+bool isFileNewer(const char *file1, const char* file2);
 void writeFld(nrs_t *nrs, dfloat t);
 void writeFld(nrs_t *nrs, dfloat t, int FP64);
 void writeFld(const char* suffix, dfloat t, int coords, int FP64,

src/io/utils.cpp

@@ -0,0 +1,35 @@
+#include <sys/stat.h>
+#include <assert.h>
+#include <fstream> 
+#include <cstdio>
+
+bool isFileNewer(const char *file1, const char* file2)
+{
+  struct stat s1, s2;
+  if (lstat(file1, &s1) != 0) assert(1);
+  if (lstat(file2, &s2) != 0) return true; 
+  if (s1.st_mtime > s2.st_mtime) 
+    return true;
+  else
+    return false;	  
+}
+
+void copyFile(const char *srcName, const char* destName)
+{
+  std::fstream src,dest;
+  src.open (srcName);
+  dest.open (destName);
+
+  std::filebuf* inbuf  = src.rdbuf();
+  std::filebuf* outbuf = dest.rdbuf();
+
+  char c = inbuf->sbumpc();
+  while (c != EOF)
+  {
+    outbuf->sputc (c);
+    c = inbuf->sbumpc();
+  }
+
+  dest.close();
+  src.close();
+}