Merge pull request #569 from lkotipal/eVlasiator-dev

eVlasiator update to AMR / 6D runs

MAKE/Makefile.carrington_gcc_openmpi_ohpc

@@ -0,0 +1,73 @@
+# module load gnu9/9.3.0                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
+# module load openmpi4/4.0.5
+# module load pmix/3.1.4
+CMP = mpicxx
+LNK = mpicxx
+#======== Vectorization ==========
+#Set vector backend type for vlasov solvers, sets precision and length.  
+#Options:  
+# AVX:      VEC4D_AGNER, VEC4F_AGNER, VEC8F_AGNER
+# AVX512:   VEC8D_AGNER, VEC16F_AGNER
+# Fallback: VEC4D_FALLBACK, VEC4F_FALLBACK, VEC8F_FALLBACK
+ifeq ($(DISTRIBUTION_FP_PRECISION),SPF)
+#Single-precision         
+       VECTORCLASS = VEC8F_AGNER
+else
+#Double-precision
+       VECTORCLASS = VEC4D_AGNER
+endif
+#======= Compiler and compilation flags =========
+# NOTES on compiler flags:
+# CXXFLAGS is for compiler flags, they are always used
+# MATHFLAGS are for special math etc. flags, these are only applied on solver functions
+# LDFLAGS flags for linker
+#-DNO_WRITE_AT_ALL:  Define to disable write at all to  
+#                    avoid memleak (much slower IO)
+#-DMPICH_IGNORE_CXX_SEEK: Ignores some multiple definition  
+#                         errors that come up when using  
+#                         mpi.h in c++ on Cray
+#
+# CXXFLAGS = -DMPICH_IGNORE_CXX_SEEK
+FLAGS =  
+#GNU flags:
+CC_BRAND = gcc
+CC_BRAND_VERSION = 9.3.0
+CXXFLAGS += -g -O3 -fopenmp -funroll-loops -std=c++17 -W -Wall -Wno-unused -m64 -march=znver2 #-flto
+testpackage: CXXFLAGS = -g -O2 -fopenmp -funroll-loops -std=c++17 -m64 -march=znver2
+MATHFLAGS = -ffast-math
+LDFLAGS = -lrt -lgfortran -std=c++17 -lgomp
+LIB_MPI = -lgomp -lmpi
+#======== PAPI ==========
+#Add PAPI_MEM define to use papi to report memory consumption?
+CXXFLAGS += -DPAPI_MEM
+#======== Allocator =========
+#Use jemalloc instead of system malloc to reduce memory fragmentation? https://github.com/jemalloc/jemalloc
+#Configure jemalloc with  --with-jemalloc-prefix=je_ when installing it
+CXXFLAGS += -DUSE_JEMALLOC -DJEMALLOC_NO_DEMANGLE
+# BOOST_VERSION = current trilinos version
+# ZOLTAN_VERSION = current trilinos verson
+#
+#======== Libraries ===========
+MPT_BRAND = OpenMPI
+MPT_VERSION = 4.0.5
+JEMALLOC_VERSION = 5.2.1
+LIBRARY_PREFIX_B = /proj/ykempf/libraries
+LIBRARY_PREFIX = /proj/group/spacephysics/libraries
+#compiled libraries mostly in modules
+LIB_PROFILE = -L$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/phiprof/lib -lphiprof -Wl,-rpath=$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/phiprof/lib
+INC_PROFILE = -I$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/phiprof/include
+LIB_VLSV = -L$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/vlsv -lvlsv -Wl,-rpath=$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/vlsv
+INC_VLSV = -I$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/vlsv
+LIB_JEMALLOC = -L$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/jemalloc/lib -ljemalloc -Wl,-rpath=$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/jemalloc/lib
+INC_JEMALLOC = -I$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/jemalloc/include
+LIB_BOOST = -L$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/boost/lib -lboost_program_options -Wl,-rpath=$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/boost/lib
+INC_BOOST = -I$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/boost/include
+LIB_ZOLTAN = -L$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/zoltan/lib -lzoltan -Wl,-rpath=$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/zoltan/lib
+INC_ZOLTAN = -I$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/zoltan/include
+LIB_PAPI = -L$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/papi/lib -lpapi -Wl,-rpath=$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/papi/lib   
+INC_PAPI = -I$(LIBRARY_PREFIX_B)/$(CC_BRAND)/$(CC_BRAND_VERSION)/$(MPT_BRAND)/$(MPT_VERSION)/carrington/papi/include
+#header libraries
+INC_EIGEN = -I$(LIBRARY_PREFIX)/ -I$(LIBRARY_PREFIX)/Eigen/
+INC_FSGRID = -I$(LIBRARY_PREFIX)/fsgrid/
+INC_DCCRG = -I$(LIBRARY_PREFIX)/dccrg/
+INC_VECTORCLASS = -I/proj/ykempf/vectorclass/version1

Makefile

@@ -196,7 +196,7 @@ OBJS = 	version.o memoryallocation.o backgroundfield.o quadr.o dipole.o linedipo
 	IPShock.o object_wrapper.o\
 	verificationLarmor.o Shocktest.o grid.o ioread.o iowrite.o vlasiator.o logger.o\
 	common.o parameters.o readparameters.o spatial_cell.o mesh_data_container.o\
-	vlasovmover.o $(FIELDSOLVER).o fs_common.o fs_limiters.o gridGlue.o
+	vlasovmover.o $(FIELDSOLVER).o fs_common.o fs_limiters.o gridGlue.o vlsvreaderinterface.o
 
 # Add Vlasov solver objects (depend on mesh: AMR or non-AMR)
 ifeq ($(MESH),AMR)

fieldsolver/gridGlue.cpp

@@ -867,3 +867,77 @@ void getdBvolFieldsFromFsGrid(
   MPI_Waitall(sendRequests.size(), sendRequests.data(), MPI_STATUSES_IGNORE);
 
 }
+
+void feedBoundaryIntoFsGrid(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
+			const std::vector<CellID>& cells,
+			FsGrid< fsgrids::technical, 2> & technicalGrid) {
+
+  int ii;
+  //sorted list of dccrg cells. cells is typicall already sorted, but just to make sure....
+  std::vector<CellID> dccrgCells = cells;
+  std::sort(dccrgCells.begin(), dccrgCells.end());
+
+  //Datastructure for coupling
+  std::map<int, std::set<CellID> > onDccrgMapRemoteProcess; 
+  std::map<int, std::set<CellID> > onFsgridMapRemoteProcess; 
+  std::map<CellID, std::vector<int64_t> >  onFsgridMapCells;
+
+  // map receive process => receive buffers 
+  std::map<int, std::vector<int> > receivedData; 
+
+  // send buffers  to each process
+  std::map<int, std::vector<int> > sendData;
+
+  //list of requests
+  std::vector<MPI_Request> sendRequests;
+  std::vector<MPI_Request> receiveRequests;
+
+  //computeCoupling
+  computeCoupling(mpiGrid, cells, technicalGrid, onDccrgMapRemoteProcess, onFsgridMapRemoteProcess, onFsgridMapCells);
+
+  // Post receives
+  receiveRequests.resize(onFsgridMapRemoteProcess.size());  
+  ii=0;
+  for(auto const &receives: onFsgridMapRemoteProcess){
+    int process = receives.first;
+    int count = receives.second.size();
+    receivedData[process].resize(count);
+    MPI_Irecv(receivedData[process].data(), count * sizeof(int),
+	      MPI_BYTE, process, 1, MPI_COMM_WORLD,&(receiveRequests[ii++]));
+  }
+
+  // Launch sends
+  ii=0;
+  sendRequests.resize(onDccrgMapRemoteProcess.size());
+  for (auto const &snd : onDccrgMapRemoteProcess){
+    int targetProc = snd.first; 
+    auto& sendBuffer=sendData[targetProc];
+    for(CellID sendCell: snd.second){
+      //Collect data to send for this dccrg cell
+      sendBuffer.push_back(mpiGrid[sendCell]->sysBoundaryFlag);
+    }
+    int count = sendBuffer.size(); //note, compared to receive this includes all elements to be sent
+    MPI_Isend(sendBuffer.data(), sendBuffer.size() * sizeof(int),
+	      MPI_BYTE, targetProc, 1, MPI_COMM_WORLD,&(sendRequests[ii]));
+    ii++;
+  }
+
+  MPI_Waitall(receiveRequests.size(), receiveRequests.data(), MPI_STATUSES_IGNORE);
+
+  for(auto const &receives: onFsgridMapRemoteProcess){
+    int process = receives.first; //data received from this process
+    int* receiveBuffer = receivedData[process].data(); // data received from process
+    for(auto const &cell: receives.second){ //loop over cellids (dccrg) for receive
+      // this part heavily relies on both sender and receiver having cellids sorted!
+      for(auto lid: onFsgridMapCells[cell]){
+        // Now save the values to face-averages
+        technicalGrid.get(lid)->sysBoundaryFlag = receiveBuffer[0];
+      }
+
+      receiveBuffer++;
+    }
+  }
+
+  MPI_Waitall(sendRequests.size(), sendRequests.data(), MPI_STATUSES_IGNORE);
+
+}

fieldsolver/gridGlue.hpp

@@ -97,6 +97,10 @@ void getDerivativesFromFsGrid(
 int getNumberOfCellsOnMaxRefLvl(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
                                 const std::vector<CellID>& cells);
 
+void feedBoundaryIntoFsGrid(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
+			const std::vector<CellID>& cells,
+			FsGrid< fsgrids::technical, 2> & technicalGrid);
+
 
 /*! Transfer field data from an FsGrid back into the appropriate CellParams slot in DCCRG 
  * \param sourceGrid Fieldsolver grid for these quantities

grid.cpp

@@ -63,10 +63,6 @@ int globalflags::AMRstencilWidth = VLASOV_STENCIL_WIDTH;
 
 extern Logger logFile, diagnostic;
 
-void initVelocityGridGeometry(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
-void initSpatialCellCoordinates(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
-void initializeStencils(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
-
 void writeVelMesh(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid) {
    const vector<CellID>& cells = getLocalCells();
 
@@ -157,7 +153,10 @@ void initializeGrids(
 
 
    phiprof::start("Refine spatial cells");
-   if(P::amrMaxSpatialRefLevel > 0 && project.refineSpatialCells(mpiGrid)) {
+   recalculateLocalCellsCache();
+   // refineSpatialCells should be a nop if amrMaxSpatialRefLevel is 0. Make this better later
+   if(project.refineSpatialCells(mpiGrid)) {
+      mpiGrid.balance_load();
       recalculateLocalCellsCache();
    }
    phiprof::stop("Refine spatial cells");
@@ -184,6 +183,9 @@ void initializeGrids(
    phiprof::start("Set spatial cell coordinates");
    initSpatialCellCoordinates(mpiGrid);
    phiprof::stop("Set spatial cell coordinates");
+
+   SpatialCell::set_mpi_transfer_type(Transfer::CELL_DIMENSIONS);
+   mpiGrid.update_copies_of_remote_neighbors(SYSBOUNDARIES_NEIGHBORHOOD_ID);
 
    phiprof::start("Initialize system boundary conditions");
    if(sysBoundaries.initSysBoundaries(project, P::t_min) == false) {
@@ -204,8 +206,8 @@ void initializeGrids(
    // Check refined cells do not touch boundary cells
    phiprof::start("Check boundary refinement");
    if(!sysBoundaries.checkRefinement(mpiGrid)) {
-      cerr << "(MAIN) ERROR: Boundary cells must have identical refinement level " << endl;
-      exit(1);
+      cerr << "(MAIN) WARNING: Boundary cells don't have identical refinement level " << endl;
+      //exit(1);
    }
    phiprof::stop("Check boundary refinement");
 
@@ -251,7 +253,8 @@ void initializeGrids(
 
    // Update technicalGrid
    technicalGrid.updateGhostCells(); // This needs to be done at some point
-
+
+   bool needCurl = false;
    if (!P::isRestart) {
       //Initial state based on project, background field in all cells
       //and other initial values in non-sysboundary cells
@@ -262,12 +265,23 @@ void initializeGrids(
       // Each initialization has to be independent to avoid threading problems 
 
       // Allow the project to set up data structures for it's setCell calls
-      bool needCurl=false;
       project.setupBeforeSetCell(cells, mpiGrid, needCurl);
-      if (needCurl==true) {
-         // Communicate the perturbed B-fields and E-fileds read from the start file over to FSgrid
+      if (needCurl) {
          feedPerBIntoFsGrid(mpiGrid, cells, perBGrid);
          perBGrid.updateGhostCells();
+      }
+   }
+
+   phiprof::start("setProjectBField");
+   project.setProjectBField(perBGrid, BgBGrid, technicalGrid);
+   // Set E field here as well?
+   perBGrid.updateGhostCells();
+   BgBGrid.updateGhostCells();
+   EGrid.updateGhostCells();
+   phiprof::stop("setProjectBField");
+
+   if (!P::isRestart) {
+      if (needCurl) {
          // E is needed only because both volumetric fields are calculated in 1 call
          feedEIntoFsGrid(mpiGrid, cells, EGrid);
          EGrid.updateGhostCells();
@@ -367,13 +381,6 @@ void initializeGrids(
       phiprof::stop("Init moments");
    }
 
-   phiprof::start("setProjectBField");
-   project.setProjectBField(perBGrid, BgBGrid, technicalGrid);
-   perBGrid.updateGhostCells();
-   BgBGrid.updateGhostCells();
-   EGrid.updateGhostCells();
-   phiprof::stop("setProjectBField");
-
    phiprof::start("Finish fsgrid setup");
    feedMomentsIntoFsGrid(mpiGrid, cells, momentsGrid,technicalGrid, false);
    if(!P::isRestart) {

grid.h

@@ -115,5 +115,8 @@ void shrink_to_fit_grid_data(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>
 bool validateMesh(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,const uint popID);
 
 void setFaceNeighborRanks( dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid );
+void initVelocityGridGeometry(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
+void initSpatialCellCoordinates(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
+void initializeStencils(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
 
 #endif