Bugfixes, now restart appears to work(!)

grid.cpp

@@ -140,6 +140,7 @@ void initializeGrid(
       cerr << "(MAIN) ERROR: System boundary conditions were not set correctly." << endl;
       exit(1);
    }
+
    phiprof::stop("Classify cells (sys boundary conditions)");
 
    if (P::isRestart) {
@@ -151,8 +152,7 @@ void initializeGrid(
       }
       phiprof::stop("Read restart");
       const vector<CellID>& cells = getLocalCells();
-
-      // Set background field, FIXME should be read in from restart
+
       #pragma omp parallel for schedule(dynamic)
       for (size_t i=0; i<cells.size(); ++i) {
          SpatialCell* cell = mpiGrid[cells[i]];

ioread.cpp

@@ -780,9 +780,8 @@ bool exec_readGrid(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
       }
    }
 
-   //Do initial load balance based on pins. No need to transfer data now since
-   //we have not yet loaded any actual data.
-   SpatialCell::set_mpi_transfer_type(Transfer::NONE);
+   //Do initial load balance based on pins. Need to transfer at least sysboundaryflags
+   SpatialCell::set_mpi_transfer_type(Transfer::ALL_SPATIAL_DATA);
    mpiGrid.balance_load(false);
 
    //update list of local gridcells
@@ -832,9 +831,7 @@ bool exec_readGrid(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
    //todo, check file datatype, and do not just use double
    phiprof::start("readCellParameters");
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"perturbed_B",CellParams::PERBX,3,mpiGrid); }
-// This has to be set anyway, there are also the derivatives that should be written/read if we want to only read in background field
-//   if(success)
-//     success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"background_B",CellParams::BGBX,3,mpiGrid);
+// Backround B has to be set, there are also the derivatives that should be written/read if we wanted to only read in background field
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"moments",CellParams::RHO,4,mpiGrid); }
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"moments_dt2",CellParams::RHO_DT2,4,mpiGrid); }
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"moments_r",CellParams::RHO_R,4,mpiGrid); }
@@ -847,9 +844,9 @@ bool exec_readGrid(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"max_v_dt",CellParams::MAXVDT,1,mpiGrid); }
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"max_r_dt",CellParams::MAXRDT,1,mpiGrid); }
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"max_fields_dt",CellParams::MAXFDT,1,mpiGrid); }
-   // Read rho losses Note: vector size = 1 (In the older versions the rho loss wasn't recorded)
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"rho_loss_adjust",CellParams::RHOLOSSADJUST,1,mpiGrid); }
    if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"rho_loss_velocity_boundary",CellParams::RHOLOSSVELBOUNDARY,1,mpiGrid); }
+// Backround B has to be set, there are also the derivatives that should be written/read if we wanted to only read in background field
    phiprof::stop("readCellParameters");
 
    phiprof::start("readBlockData");

sysboundary/sysboundary.cpp

@@ -570,10 +570,17 @@ void SysBoundary::applySysBoundaryVlasovConditions(
 
 /*! Get a pointer to the SysBoundaryCondition of given index.
  * \param sysBoundaryType Type of the system boundary condition to return
- * \return Pointer to the instance of the SysBoundaryCondition.
+ * \return Pointer to the instance of the SysBoundaryCondition. NULL if sysBoundaryType is invalid.
  */
 SBC::SysBoundaryCondition* SysBoundary::getSysBoundary(cuint sysBoundaryType) const {
-   return indexToSysBoundary.find(sysBoundaryType)->second;
+   auto it = indexToSysBoundary.find(sysBoundaryType);
+   if(it != indexToSysBoundary.end()){
+      return it->second;
+   }
+   else{
+      std::cerr << "System boundary "<< sysBoundaryType << " is invalid  " << __FILE__ << ":" << __LINE__ << std::endl;
+      return NULL;
+   }
 }
 
 /*! Get the number of SysBoundaryConditions stored in SysBoundary.

velocity_mesh_old.h

@@ -589,7 +589,8 @@ namespace vmesh {
    bool VelocityMesh<GID,LID>::push_back(const std::vector<GID>& blocks) {
       if (size()+blocks.size() >= meshParameters[meshID].max_velocity_blocks) {
          std::cerr << "vmesh: too many blocks, current size is " << size();
-         std::cerr << " max " << meshParameters[meshID].max_velocity_blocks << std::endl;
+         std::cerr << ", adding " << blocks.size() << " blocks";
+         std::cerr << ", max is " << meshParameters[meshID].max_velocity_blocks << std::endl;
          return false;
       }