Hypergraph

datareduction/datareducer.cpp

@@ -439,7 +439,7 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti
       }
       if(P::systemWriteAllDROs || lowercase == "lbweight" || lowercase == "vg_lbweight" || lowercase == "vg_loadbalanceweight" || lowercase == "vg_loadbalance_weight") {
          // Load balance metric for LB debugging
-         outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_loadbalance_weight",CellParams::LBWEIGHTCOUNTER,1));
+         outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_loadbalance_weight",CellParams::LBWEIGHTCOUNTER, 1));
          outputReducer->addMetadata(outputReducer->size()-1,"","","$\\mathrm{LB weight}$","");
          if(!P::systemWriteAllDROs) {
             continue;
@@ -3664,7 +3664,7 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti
          }
       }
       if(P::diagnosticWriteAllDROs || lowercase == "lbweight" || lowercase == "vg_lbweight" || lowercase == "vg_loadbalanceweight" || lowercase == "vg_loadbalance_weight" || lowercase == "loadbalance_weight") {
-         diagnosticReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_loadbalance_weight",CellParams::LBWEIGHTCOUNTER,1));
+         diagnosticReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_loadbalance_weight",CellParams::LBWEIGHTCOUNTER, 1));
          if(!P::diagnosticWriteAllDROs) {
             continue;
          }

grid.cpp

@@ -66,7 +66,6 @@ using namespace std;
 
 extern Logger logFile;
 
-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);
 
@@ -143,11 +142,24 @@ void initializeGrids(
                     sysBoundaries.isPeriodic(2))
       .initialize(comm)
       .set_geometry(geom_params);
+
+   mpiGrid.set_load_balance_norm(P::loadBalanceNorm);
+
+   for (const auto& [key, value] : P::loadBalanceOptions) {
+      mpiGrid.set_partitioning_option(key, value);
+   }
+
+   mpiGrid.set_partitioning_option("LB_APPROACH", "PARTITION");
+
+   // Hypergraph partitioning needs stencils initialized
+   initializeStencils(mpiGrid);
    dccrgTimer.stop();
 
    phiprof::Timer refineTimer {"Refine spatial cells"};
-   // We need this first as well
    recalculateLocalCellsCache(mpiGrid);
+
+   setPartitioningNeighborhoods(mpiGrid);
+
    if (!P::isRestart) {
       // Note call to project.refineSpatialCells below
       if (P::amrMaxSpatialRefLevel > 0 && project.refineSpatialCells(mpiGrid)) {
@@ -166,13 +178,17 @@ void initializeGrids(
       }
    }
    refineTimer.stop();
-   initializeStencils(mpiGrid);
+
+   mpiGrid.set_partitioning_option("OBJ_WEIGHTS_COMPARABLE", "1");
 
    for (const auto& [key, value] : P::loadBalanceOptions) {
       mpiGrid.set_partitioning_option(key, value);
    }
    phiprof::Timer initialLBTimer {"Initial load-balancing"};
    if (myRank == MASTER_RANK) logFile << "(INIT): Starting initial load balance." << endl << writeVerbose;
+
+   // TODO: do we really need two initial LB?
+   setPartitioningNeighborhoods(mpiGrid);
    mpiGrid.balance_load(); // Direct DCCRG call, recalculate cache afterwards
    recalculateLocalCellsCache(mpiGrid);
 
@@ -316,7 +332,7 @@ void initializeGrids(
    } else if (P::writeFullBGB) {
       // If, instead of starting a regular simulation, we are only writing out the background field, it is enough to set a dummy load balance value of 1 here.
       for (size_t i=0; i<cells.size(); ++i) {
-         mpiGrid[cells[i]]->parameters[CellParams::LBWEIGHTCOUNTER] = 1;
+         mpiGrid[cells[i]]->parameters[CellParams::LBWEIGHTCOUNTER] = 1.0;
       }
    }
 
@@ -398,6 +414,8 @@ void initializeGrids(
       P::dt = P::bailout_min_dt;
    }
 
+   mpiGrid.set_partitioning_option("LB_APPROACH", P::loadBalanceOptions.count("LB_APPROACH") ? P::loadBalanceOptions["LB_APPROACH"] : "REPARTITION");
+
    // With all cell data in place, make preparations for translation
    prepareAMRLists(mpiGrid);
    initialStateTimer.stop();
@@ -475,6 +493,24 @@ void setFaceNeighborRanks( dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>&
    }
 }
 
+void setPartitioningNeighborhoods(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid) {
+   const vector<CellID>& cells = getLocalCells();
+   for (auto& cell : cells){
+      mpiGrid.clear_partitioning_neighborhoods(cell);
+      for (auto neighborhood : P::partitioningNeighborhoods) {
+         // TODO: for now, layer 1 cells communicate in the extended neighborhood
+         // If this is ever fixed, SYSBOUNDRIES_NEIGHBORHOOD_ID needs a case
+         if (neighborhood == Neighborhoods::SYSBOUNDARIES_EXTENDED) {
+            if (mpiGrid[cell]->sysBoundaryLayer == 1 || mpiGrid[cell]->sysBoundaryLayer == 2) {
+               mpiGrid.add_partitioning_neighborhood(cell, neighborhood);
+            }
+         } else {
+            mpiGrid.add_partitioning_neighborhood(cell, neighborhood);
+         }
+      }
+   }
+}
+
 void balanceLoad(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid, SysBoundary& sysBoundaries, FsGrid<fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, bool doTranslationLists){
    // Invalidate cached cell lists
    Parameters::meshRepartitioned = true;
@@ -488,6 +524,7 @@ void balanceLoad(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid, S
    deallocTimer.stop();
 
    //set weights based on each cells LB weight counter
+   setPartitioningNeighborhoods(mpiGrid);
    const vector<CellID>& cells = getLocalCells();
    for (size_t i=0; i<cells.size(); ++i){
       // Set cell weight. We could use different counters or number of blocks if different solvers are active.

grid.h

@@ -51,6 +51,8 @@ void initializeGrids(
    Project& project
 );
 
+void setPartitioningNeighborhoods(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
+
 /*!
   \brief Balance load
 

ioread.cpp

@@ -1107,6 +1107,28 @@ bool readIonosphereNodeVariable(
    return true;
 }
 
+/*!
+\brief read length of SpatialGrid vector variable
+\param file The ParallelReader file to use
+\param variableName Name of the variable to check
+*/
+int readVectorSize(vlsv::ParallelReader& file,
+                   const std::string& variableName){
+   uint64_t arraySize;
+   uint64_t vectorSize;
+   vlsv::datatype::type dataType;
+   uint64_t byteSize;
+   list<pair<string,string> > attribs;
+   attribs.push_back(make_pair("name",variableName));
+   attribs.push_back(make_pair("mesh","SpatialGrid"));
+
+   if (file.getArrayInfo("VARIABLE",attribs,arraySize,vectorSize,dataType,byteSize) == false) {
+      logFile << "(RESTART)  ERROR: Failed to read " << endl << write;
+      return -1;
+   }
+   return (int)vectorSize;
+}
+
 /*!
 \brief Read in state from a vlsv file in order to restart simulations
 \param mpiGrid Vlasiator's grid

parameters.cpp

@@ -138,7 +138,9 @@ bool P::isRestart = false;
 int P::writeAsFloat = false;
 int P::writeRestartAsFloat = false;
 string P::loadBalanceAlgorithm = string("");
+std::vector<int> P::partitioningNeighborhoods {};
 std::map<std::string, std::string> P::loadBalanceOptions;
+int P::loadBalanceNorm {1};
 uint P::rebalanceInterval = numeric_limits<uint>::max();
 
 vector<string> P::outputVariableList;
@@ -404,8 +406,11 @@ bool P::addParameters() {
    RP::add("vlasovsolver.GhostTranslateExtent","Stencil size in all-local ghost translation (default: VLASOV_STENCIL_WIDTH+1",0);
 
    // Load balancing parameters
-   RP::add("loadBalance.algorithm", "Load balancing algorithm to be used", string("RCB"));
+   RP::add("loadBalance.algorithm", "Load balancing algorithm to be used", string("RIB"));
+   RP::addComposing("loadBalance.partitioning_neighborhood", "Neighborhood ID for (hyper)graph partitioning, see definitions.h");
    RP::add("loadBalance.tolerance", "Load imbalance tolerance", string("1.05"));
+   RP::add("loadBalance.weight_dim", "Dimension of object weight", string("1"));
+   RP::add("loadBalance.norm", "Norm to use for cell weights in load balance, default 1-norm (sum)", 1);
    RP::add("loadBalance.rebalanceInterval", "Load rebalance interval (steps)", 10);
 
    RP::addComposing("loadBalance.optionKey", "Zoltan option key. Has to be matched by loadBalance.optionValue.");
@@ -1058,10 +1063,20 @@ void Parameters::getParameters() {
    }
    // Get load balance parameters
    RP::get("loadBalance.algorithm", P::loadBalanceAlgorithm);
+   RP::get("loadBalance.partitioning_neighborhood", P::partitioningNeighborhoods);
+
+   if (P::partitioningNeighborhoods.empty()) {
+      P::partitioningNeighborhoods.push_back(Neighborhoods::FULL);
+   }
+
    loadBalanceOptions["IMBALANCE_TOL"] = "";
    RP::get("loadBalance.tolerance", loadBalanceOptions["IMBALANCE_TOL"]);
+   RP::get("loadBalance.norm", P::loadBalanceNorm);
    RP::get("loadBalance.rebalanceInterval", P::rebalanceInterval);
 
+   loadBalanceOptions["OBJ_WEIGHT_DIM"] = "";
+   RP::get("loadBalance.weight_dim", loadBalanceOptions["OBJ_WEIGHT_DIM"]);
+
    std::vector<std::string> loadBalanceKeys;
    std::vector<std::string> loadBalanceValues;
    RP::get("loadBalance.optionKey", loadBalanceKeys);

parameters.h

@@ -157,7 +157,9 @@ struct Parameters {
                                    in the Lorentz force and in the field solver.*/
 
    static std::string loadBalanceAlgorithm; /*!< Algorithm to be used for load balance.*/
+   static std::vector<int> partitioningNeighborhoods; /*!< Partitioning neighborhood for (hyper)graph.*/
    static std::map<std::string, std::string> loadBalanceOptions;  // Other Load balancing options
+   static int loadBalanceNorm;   // Norm to use for cell weights
    static uint rebalanceInterval;           /*!< Load rebalance interval (steps). */
    static bool prepareForRebalance; /**< If true, propagators should measure their time consumption in preparation
                                      * for mesh repartitioning.*/