Grid controller manages adaptation and load balancing in multi-solver environment. More...

#include <cartesiangridcontroller.h>

Collaboration diagram for maia::grid::Controller< nDim >:

Public Types
using	Grid = CartesianGrid< nDim >

using	Cell = typename CartesianGrid< nDim >::Cell

using	GridProxy = typename maia::grid::Proxy< nDim >

Public Member Functions
	Controller (Grid grid_, std::vector< std::unique_ptr< Solver > > solvers, std::vector< std::unique_ptr< Coupling > > *couplers)

	~Controller ()

MBool	balance (const MBool force=false, const MBool finalTimeStep=false, const MBool adaptation=false)

MBool	adaptation (const MBool force=false)
	performs mesh adaptation More...

void	writeRestartFile (const MBool, const MBool)

void	printDomainStatistics (const MString &status="")
	Print statistics regarding the cell distribution among domains. More...

void	savePartitionFile ()

MBool	updateGridPartitionWorkloads ()
	Update the partition cell workloads in the grid file using user specified or default weights for the solver load quantities. More...

void	castAdaptationIntervalToMultipleOfCoarsestTimeStep (MInt maxLevel, MInt maxUniformRefinementLevel)

void	logTimerStatistics (const MString &)

MBool	isDlbTimeStep ()

MBool	isAdaptationTimeStep ()

Private Member Functions
Grid &	gridb ()

const Grid &	gridb () const

MInt	domainId ()

MInt	noDomains ()

Solver &	solver (const MInt solverId)

const Solver &	solver (const MInt solverId) const

MInt	noSolvers () const

MInt	solverLocalRootDomain (Solver *const solver)
	Determine the global domain id of the solver local root domain. More...

Coupling &	coupler (const MInt couplerId)

const Coupling &	coupler (const MInt couplerId) const

MInt	noCouplers () const

void	initDlbProperties ()
	Read Dynamic Load Balancing properties. More...

void	initTimers ()
	Initialize timers. More...

MBool	needLoadBalancing (const MFloat localRunTime, MFloat *const loads, MFloat &imbalance)
	Return if dynamic load balancing is needed and compute domain loads. More...

void	partition (MLong partitionCellOffsets, MLong globalIdOffsets, const MBool onlyPartitionOffsets)
	Compute new domain decomposition. More...

MBool	loadBalancingPartition (const MFloat loads, const MFloat imbalance, MLong const partitionCellOffsets, MLong *const globalIdOffsets)
	Determine new partitioning for dynamic load balancing. More...

void	updateWeightsAndWorkloads (const std::vector< MFloat > &weights, const MBool restoreDefaultWeights)
	Determine the cell weights using the given weighting factors for the different solver load quantities and update the partition cell workloads. More...

void	getSpecifiedSolverWeights (std::vector< MFloat > &weights)
	Return the specified (or default) solver weights for all solvers. More...

void	accumulateCellWeights ()
	This function handels the setCellWeights-functionality of all solvers and writes the accumulated-cellweight of all solvers into the grid. More...

void	storeTimings ()
	Store timings of all timesteps for all domains for performance evaluations (i.e. runtime and idle time). More...

void	storeLoadsAndWeights (const MFloat const loads, const MInt noLoadTypes, const MInt const loadQuantities, const MFloat domainWeight, const MFloat *const weights)
	Store domain loads and additional information to file. More...

void	loadBalancingCalcNewGlobalOffsets (const MLong const oldPartitionCellOffsets, const MLong const newPartitionCellOffsets, MLong *const globalOffsets)
	Calculate new global domain offsets given the current and new global partition cell offsets. More...

void	loadBalancingCalcNoCellsToSend (const MLong const offsets, MInt const noCellsToSendByDomain, MInt const noCellsToReceiveByDomain, MInt const sortedCellId, MInt *const bufferIdToCellId)
	Based on new domain offsets calculate the number of cells to send/receive to/from each domain. More...

void	computeWeights (const MFloat *loads, const MFloat domainWeight, std::vector< MFloat > &weights)
	Compute computational weights for different components in the simulation based on the current distribution and the domain loads. More...

MInt	globalNoLoadTypes ()
	Return global number of load types of all solvers. More...

void	getLoadQuantities (MInt *const loadQuantities)
	Return load quantities of all solvers on this domain. More...

void	estimateParameters (MInt m, MInt n, const MFloat const A, const MFloat const b, MFloat *const x)
	Solve the parameter estimation problem A*x=b. More...

void	communicateGlobalSolverVars (Solver *const solver)

void	determineDataSizesDlb (const MInt solverId, const MInt mode, const MInt const noCellsToSend, const MInt const bufferIdToCellId, std::vector< std::vector< MInt > > &sendSizeVector, std::vector< std::vector< MInt > > &recvSizeVector)
	Determine the data sizes for a solver/coupler that need to be communicated during load balancing. More...

void	redistributeDataDlb (const MInt id, const MInt mode, std::vector< MInt > &sendSizeVector, std::vector< MInt > &recvSizeVector, const MInt const bufferIdToCellId, const MInt noCells, std::vector< MInt > &intDataRecv, std::vector< MLong * > &longDataRecv, std::vector< MFloat * > &floatDataRecv, std::vector< MInt > &dataTypes)
	Communicate all solver data for load balancing according to the send/recv sizes. More...

void	setDataDlb (const MInt solverId, const MInt mode, std::vector< MInt * > &intDataRecv, std::vector< MLong * > &longDataRecv, std::vector< MFloat * > &floatDataRecv, std::vector< MInt > &dataTypes, const MBool freeMemory)
	Set the solver/coupler data for load balancing. More...

void	resetAllTimer ()

std::vector< std::function< void(const MInt)> >	refineCellVec ()

std::vector< std::function< void(const MInt)> >	removeChildsVec ()

std::vector< std::function< void(const MInt, const MInt)> >	swapProxyVec ()

std::vector< std::function< MInt(const MFloat *, const MInt, const MInt)> >	cellOutsideVec ()

std::vector< std::function< void()> >	resizeGridMapVec ()

std::vector< std::function< void(const MInt)> >	removeCellVec ()

Private Attributes
Grid *	m_grid

const std::vector< std::unique_ptr< Solver > > *	m_solvers

const std::vector< std::function< void(const MInt)> >	m_refineCellSolver

const std::vector< std::function< void(const MInt)> >	m_removeChildsSolver

const std::vector< std::function< void(const MInt, const MInt)> >	m_swapProxySolver

const std::vector< std::function< MInt(const MFloat *, const MInt, const MInt)> >	m_cellOutside

const std::vector< std::function< void()> >	m_resizeGridMapSolver

const std::vector< std::function< void(const MInt)> >	m_removeCellSolver

const std::vector< std::unique_ptr< Coupling > > *	m_couplers

MBool	m_adaptation = false

MInt	m_adaptationInterval = 0

MInt	m_adaptationStart = -1

MInt	m_adaptationStop = -1

MInt	m_lastAdaptationTimeStep = 0

MInt	m_noAdaptations = 0

MInt	m_noMaxAdaptations = -1

MInt	m_noSensors {}

MBool	m_balance = false

MInt	m_forceBalance = 0

MBool	m_performanceOutput = false

MBool	m_debugBalance = false

MInt	m_loadBalancingMode = -1

MInt	m_loadBalancingTimerMode = 0

MInt	m_loadBalancingInterval = -1

MInt	m_loadBalancingOffset = 0

MInt	m_loadBalancingStartTimeStep = 0

MInt	m_loadBalancingStopTimeStep = -1

MInt	m_loadBalancingTimerStartOffset = 0

MBool	m_forceLoadBalancing = false

MBool	m_testDynamicLoadBalancing = false

MBool	m_testUpdatePartitionCells = false

MInt	m_testUpdatePartCellsOffspringThreshold = 10

MFloat	m_testUpdatePartCellsWorkloadThreshold = 100.0

MInt	m_nAdaptationsSinceBalance = 9999999

MBool	m_syncTimeSteps = false

MBool	m_syncTimerSteps = true

MInt	m_dlbStep = 0

MInt	m_lastLoadBalancingTimeStep = 0

MInt	m_dlbLastResetTimeStep = 0

MFloat	m_dlbImbalanceThreshold = 0.05

MFloat	m_maxPerformanceVarThreshold = 0.15

MInt	m_dlbPartitionMethod = DLB_PARTITION_DEFAULT

MBool	m_dlbUpdatePartitionCells = false

MBool	m_saveGridNewPartitionCells = false

MBool	m_dlbSmoothGlobalShifts = true

MInt	m_dlbNoLocalShifts = 0

MInt	m_dlbNoFinalLocalShifts = 0

MFloat	m_dlbMaxWorkloadLimit = 1.5

MFloat	m_dlbPreviousLocalRunTime = 0.0

MFloat	m_dlbPreviousLocalIdleTime = 0.0

std::vector< MFloat >	m_dlbTimings

MBool	m_outputDlbTimings = false

std::vector< MInt >	m_dlbTimeStepsAll

std::vector< MFloat >	m_dlbRunTimeAll

std::vector< MFloat >	m_dlbIdleTimeAll

std::vector< MInt >	m_previousLoadBins {}

std::array< MFloat, 5 >	m_previousLoadInfo {}

MInt	m_previousLoadInfoStep = -1

std::vector< MFloat >	m_domainWeights

MBool	m_useDomainFactor = false

std::vector< MInt >	m_lastOffsetShiftDirection {}

MFloat	m_timePerStepTotal = std::numeric_limits<MFloat>::max()

MLong	m_optPartitionCellOffsetTotal = -1

MFloat	m_imbalance = std::numeric_limits<MFloat>::max()

MBool	m_balanceAfterAdaptation = false

MInt	m_balanceAdaptationInterval = 1

MBool	m_dlbRestartWeights = false

MFloat *	m_dlbLastWeights = nullptr

MFloat *	m_dlbStaticWeights = nullptr

MInt	m_dlbStaticWeightMode = -1

MBool	m_timersInitialized = false

MInt	m_timerGroup = -1

std::array< MInt, Timers::_count >	m_timers {}

std::vector< MInt >	m_solverTimerGroups {}

std::vector< std::array< MInt, SolverTimers::_count > >	m_solverTimers {}

MString	m_currentGridFileName

MBool	m_useNonSpecifiedRestartFile = false

MString	m_outputDir

MInt *	m_recalcIds

Detailed Description

template<MInt nDim>
class maia::grid::Controller< nDim >

Definition at line 88 of file cartesiangridcontroller.h.

Member Typedef Documentation

◆ Cell

template<MInt nDim>

using maia::grid::Controller< nDim >::Cell = typename CartesianGrid<nDim>::Cell

Definition at line 91 of file cartesiangridcontroller.h.

◆ Grid

template<MInt nDim>

using maia::grid::Controller< nDim >::Grid = CartesianGrid<nDim>

Definition at line 90 of file cartesiangridcontroller.h.

◆ GridProxy

template<MInt nDim>

using maia::grid::Controller< nDim >::GridProxy = typename maia::grid::Proxy<nDim>

Definition at line 92 of file cartesiangridcontroller.h.

Constructor & Destructor Documentation

◆ Controller()

template<MInt nDim>

maia::grid::Controller< nDim >::Controller	(	Grid *	grid_,
		std::vector< std::unique_ptr< Solver > > *	solvers,
		std::vector< std::unique_ptr< Coupling > > *	couplers
	)

Definition at line 415 of file cartesiangridcontroller.h.

◆ ~Controller()

template<MInt nDim>

maia::grid::Controller< nDim >::~Controller ( )

inline

Definition at line 97 of file cartesiangridcontroller.h.

                {
    if(m_timersInitialized) { // Timers are not created if there is no Cartesian grid
      RECORD_TIMER_STOP(m_timers[Timers::Controller]);
    }
  }

Member Function Documentation

◆ accumulateCellWeights()

template<MInt nDim>

void maia::grid::Controller< nDim >::accumulateCellWeights

private

Author: Tim Wegmann

Definition at line 2219 of file cartesiangridcontroller.h.

                                             {
  TRACE();
 
  const MInt noGridCells = gridb().treeb().size();
 
  // Reset cell weights
  for(MInt gridCellId = 0; gridCellId < noGridCells; gridCellId++) {
    gridb().treeb().weight(gridCellId) = 0.0;
  }
 
  if(!m_dlbRestartWeights || m_dlbLastWeights == nullptr) {
    // use static cell weights
 
    // TODO labels:CONTROLLER why not use an array of size noGridCells, set to zero
    // and add cell weights of each solver?
    const MInt sizeTimesSolver = noGridCells * noSolvers();
    MFloatScratchSpace solverweight(sizeTimesSolver, FUN_, "solverweight");
    // might be necessary for dimensionless weight-comparison between different-solvers,
    // to devide the weight of ech solver-cell by its maximum-possible-value!
    // MFloatScratchSpace solvermaxweights(noSolvers(),FUN_, "solvermaxweights");
    solverweight.fill(F0);
 
    for(MInt i = 0; i < noSolvers(); i++) {
      if(solver(i).isActive()) { // Only for solvers active on this ranks
        solver(i).setCellWeights(&solverweight[0]);
      }
    }
 
    // Accumulate weights of all solvers
    for(MInt i = 0; i < noSolvers(); i++) {
      for(MInt gridCellId = 0; gridCellId < noGridCells; gridCellId++) {
        const MInt id = gridCellId + noGridCells * i;
        gridb().treeb().weight(gridCellId) += solverweight[id];
      }
    }
 
 
  } else {
    // use weights from last DLB
    MInt weightOffset = 0;
    for(MInt solverId = 0; solverId < noSolvers(); solverId++) {
      if(solver(solverId).isActive()) {
        for(MInt cellId = 0; cellId < noGridCells; cellId++) {
          if(gridb().a_hasProperty(cellId, Cell::IsHalo)) {
            continue;
          }
          gridb().a_weight(cellId) += solver(solverId).getCellLoad(cellId, &m_dlbLastWeights[weightOffset]);
        }
      }
      weightOffset += solver(solverId).noLoadTypes();
    }
  }
}

◆ adaptation()

template<MInt nDim>

MBool maia::grid::Controller< nDim >::adaptation ( const MBool force = false )

Author: Lennart Schneiders

Definition at line 1994 of file cartesiangridcontroller.h.

                                                    {
  if(m_grid == nullptr) {
    return false;
  }
 
  // Reset adaptation status
  gridb().m_wasAdapted = false;
 
  if(!m_adaptation
     || (!force
         && ((m_adaptationInterval <= 0 || globalTimeStep % m_adaptationInterval != 0)
             || globalTimeStep < m_adaptationStart))) {
    return false;
  }
 
  RECORD_TIMER_START(m_timers[Timers::Adaptation]);
 
  cerr0 << "=== Initialising solver-adaptation at time-step " << globalTimeStep << std::endl;
  // TODO labels:TIMERS,ADAPTATION adaptation timers for new split version -> initTimers()
  logTimerStatistics("before adaptation");
 
  // Disable all DLB timers
  const MInt noDlbTimers = maia::dlb::g_dlbTimerController.noDlbTimers();
  ScratchSpace<MBool> dlbTimersStatus(std::max(noDlbTimers, 1), AT_, "dlbTimersStatus");
  maia::dlb::g_dlbTimerController.disableAllDlbTimers(&dlbTimersStatus[0]);
 
  // TODO labels:CONTROLLER remove if all solvers have splitAdaptation
  MBool splitAdaptation = true;
  for(MInt i = 0; i < noSolvers(); i++) {
    if(!solver(i).m_splitAdaptation) {
      std::cerr << "Update Adaptation to splitAdaptation! solver " << i << std::endl;
      splitAdaptation = false;
    }
  }
 
  // Cancel any open MPI requests in the solver since these might conflict with communication
  // during adaptation
  for(MInt i = 0; i < noSolvers(); i++) {
    solver(i).cancelMpiRequests();
  }
 
  if(splitAdaptation) { // new split-Adaptation version
 
    // prepare Adaptations
    for(MInt i = 0; i < noSolvers(); i++) {
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Total]);
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Adaptation]);
      solver(i).prepareAdaptation();
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Adaptation]);
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Total]);
    }
 
    for(MInt i = 0; i < noCouplers(); i++) {
      const MInt timerId = noSolvers() + i;
      RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Total]);
      RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Adaptation]);
      coupler(i).prepareAdaptation();
      RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Adaptation]);
      RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Total]);
    }
 
    // loop over all available levels (and a specified additional number of steps)!
    // prefered over while-loop!
    MInt addedAdaptationSteps = 0;
    addedAdaptationSteps = Context::getBasicProperty<MInt>("addedAdaptationSteps", AT_, &addedAdaptationSteps);
    for(MInt level = gridb().maxUniformRefinementLevel(); level < gridb().maxRefinementLevel() + addedAdaptationSteps;
        level++) {
      const MLong oldLocalCnt = gridb().noInternalCells();
      const MLong oldCnt = gridb().noCellsGlobal();
      std::vector<std::vector<MFloat>> sensors;
      std::vector<std::bitset<CartesianGrid<nDim>::m_maxNoSensors>> sensorCellFlag(gridb().noInternalCells());
      std::vector<MFloat> sensorWeight;
      std::vector<MInt> sensorSolverId;
 
      // set sensors
      for(MInt i = 0; i < noSolvers(); i++) {
        solver(i).setSensors(sensors, sensorWeight, sensorCellFlag, sensorSolverId);
      }
 
      ASSERT(sensors.size() < CartesianGrid<nDim>::m_maxNoSensors, "Increase bitsetsize!");
 
      // saveSensorData
      // Determine of at least one solver wants to save its sensor data
      MBool saveSensorData = false;
      for(MInt i = 0; i < noSolvers(); i++) {
        saveSensorData = saveSensorData || solver(i).m_saveSensorData;
      }
      if(saveSensorData) {
        // Write grid file
        std::stringstream gridFileName;
        gridFileName << "sensorDataGrid_" << std::to_string(level) << "_" << globalTimeStep << ParallelIo::fileExt();
        gridb().saveGrid((m_outputDir + gridFileName.str()).c_str(), m_recalcIds);
        // Write sensor data
        for(MInt i = 0; i < noSolvers(); i++) {
          if(solver(i).m_saveSensorData) {
            solver(i).saveSensorData(sensors, level, gridFileName.str(), m_recalcIds);
          }
        }
      }
 
      RECORD_TIMER_START(m_timers[Timers::MeshAdaptation]);
      gridb().meshAdaptation(sensors, sensorWeight, sensorCellFlag, sensorSolverId, m_refineCellSolver,
                             m_removeChildsSolver, m_removeCellSolver, m_swapProxySolver, m_cellOutside,
                             m_resizeGridMapSolver);
      RECORD_TIMER_STOP(m_timers[Timers::MeshAdaptation]);
 
      // compacts cells, and emidiate after adaptation
      for(MInt i = 0; i < noSolvers(); i++) {
        RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Total]);
        RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Adaptation]);
        solver(i).postAdaptation();
        RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Adaptation]);
        RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Total]);
      }
 
      // call couplers only after all solvers have rebild their grid-structure!
      for(MInt j = 0; j < noCouplers(); j++) {
        const MInt timerId = noSolvers() + j;
        RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Total]);
        RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Adaptation]);
        coupler(j).postAdaptation();
        RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Adaptation]);
        RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Total]);
      }
 
      m_log << "Mesh adaptation: " << gridb().noCellsGlobal() - oldCnt << " cells generated"
            << " (before: " << oldCnt << ", now: " << gridb().noCellsGlobal() << ")." << std::endl;
 
      MBool skipLoop = false;
      for(MInt i = 0; i < noSolvers(); i++) {
        if(solver(i).m_singleAdaptation) {
          skipLoop = true;
        }
      }
 
      if(skipLoop && globalTimeStep > 0) break;
 
      if(globalTimeStep < 0 && m_loadBalancingInterval > 0) {
        MInt balanceTrigger = 0;
        MInt deltaCells = gridb().noInternalCells() - oldLocalCnt;
        if(gridb().noInternalCells() + deltaCells * gridb().m_maxNoChilds > gridb().maxNoCells()) {
          balanceTrigger = 1;
        }
        MPI_Allreduce(MPI_IN_PLACE, &balanceTrigger, 1, MPI_INT, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                      "balanceTrigger");
        if(balanceTrigger || level == (gridb().maxRefinementLevel() - 1)) {
          cerr0 << "Performing intermediate load balancing step." << std::endl;
#ifdef MAIA_GRID_SANITY_CHECKS
          gridb().gridSanityChecks();
          gridb().checkWindowHaloConsistency(true);
#endif
          // preserve loadBalancing mode and
          // switch to weight based mode for the balance during initial adaptation!
          const MInt backUp = m_loadBalancingMode;
          m_loadBalancingMode = 0;
          balance(true, false, true);
          m_loadBalancingMode = backUp;
        }
      }
    }
 
    // reinit solver after adaptation
    for(MInt i = 0; i < noSolvers(); i++) {
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Total]);
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Adaptation]);
      solver(i).finalizeAdaptation();
      for(MInt j = 0; j < noCouplers(); j++) {
        coupler(j).finalizeAdaptation(i);
      }
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Adaptation]);
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Total]);
    }
 
  } else { // TODO labels:CONTROLLER delete below
    const MLong oldCnt = gridb().noCellsGlobal();
    std::vector<std::vector<MFloat>> sensors;
    std::vector<std::bitset<CartesianGrid<nDim>::m_maxNoSensors>> sensorCellFlag(gridb().noInternalCells());
    std::vector<MFloat> sensorWeight;
    std::vector<MInt> sensorSolverId;
 
    for(MInt i = 0; i < noSolvers(); i++) {
      std::cerr << "preparing  adaptation for solver " << i << std::endl;
      solver(i).prepareAdaptation(sensors, sensorWeight, sensorCellFlag, sensorSolverId);
    }
    ASSERT(sensors.size() < CartesianGrid<nDim>::m_maxNoSensors, "Increase bitset size!");
 
    gridb().meshAdaptation(sensors, sensorWeight, sensorCellFlag, sensorSolverId, m_refineCellSolver,
                           m_removeChildsSolver, m_removeCellSolver, m_swapProxySolver, m_cellOutside,
                           m_resizeGridMapSolver);
 
    for(MInt i = 0; i < noSolvers(); i++) {
      // in multisolver case: translate tree ids back to solver-local ids
      solver(i).reinitAfterAdaptation();
    }
 
    m_log << "Mesh adaptation: " << gridb().noCellsGlobal() - oldCnt << " cells generated"
          << " (before: " << oldCnt << ", now: " << gridb().noCellsGlobal() << ")." << std::endl;
  }
 
  // Enable all previously active DLB timers again
  maia::dlb::g_dlbTimerController.enableAllDlbTimers(&dlbTimersStatus[0]);
 
  // reset all timers after an adaptation:
  // Tim version
  // resetAllTimer();
 
  // set the last adaptation timeStep
  m_lastAdaptationTimeStep = globalTimeStep;
  m_nAdaptationsSinceBalance = m_nAdaptationsSinceBalance + 1;
 
  printDomainStatistics("after adaptation");
 
  cerr0 << "=== Finished adaptation" << std::endl;
  m_log << "Finished adaptation at time-step " << globalTimeStep << std::endl;
  RECORD_TIMER_STOP(m_timers[Timers::Adaptation]);
 
  return true;
}

◆ balance()

template<MInt nDim>

MBool maia::grid::Controller< nDim >::balance	(	const MBool	force = `false`,
		const MBool	finalTimeStep = `false`,
		const MBool	adaptation = `false`
	)

Check if load balancing is necessary and if yes, call appropriate grid & solver methods.

Return true if actual load balancing took place.

Definition at line 1057 of file cartesiangridcontroller.h.

                                                                                                    {
  // Return if not enabled or in serial
  if(!m_balance || noDomains() == 1) {
    return false;
  }
 
  // Reset balance status
  gridb().m_wasBalanced = false;
 
  std::vector<std::pair<MFloat, MString>> durations{};
  auto logDuration = [&durations](const MFloat time, const MString comment, const MBool fromTimer = false) {
    const MFloat duration = (fromTimer) ? time : wallTime() - time;
    durations.push_back(std::make_pair(duration, comment));
  };
 
  RECORD_TIMER_START(m_timers[Timers::DLB]);
 
  const MFloat dlbStartTime = wallTime();
  const MInt oldAllocatedBytes = allocatedBytes();
 
  const MInt noDlbTimers = maia::dlb::g_dlbTimerController.noDlbTimers();
  if(noDlbTimers == 0 && m_loadBalancingMode == 1 && !m_testDynamicLoadBalancing) {
    std::cerr << "There are no DLB timers, but loadBalancingMode is 1 and testing is off; "
                 "switching to loadBalancingMode = 0."
              << std::endl;
    m_loadBalancingMode = 0;
  }
 
  // Accumulate timer records of all dlb timers
  MFloat localRunTime = 0.0;
  MFloat localIdleTime = 0.0;
 
  for(MInt i = 0; i < noDlbTimers; i++) {
    const MFloat loadRecord = maia::dlb::g_dlbTimerController.returnLoadRecord(i, m_loadBalancingTimerMode);
    const MFloat idleRecord = maia::dlb::g_dlbTimerController.returnIdleRecord(i, m_loadBalancingTimerMode);
 
    if(m_loadBalancingMode == 1 && !m_testDynamicLoadBalancing && (loadRecord < 0.0 || idleRecord < 0.0)) {
      TERMM(1, "Load/Idle record for dlb timer #" + std::to_string(i) + " is less than zero on global domain #"
                   + std::to_string(domainId()) + ": " + std::to_string(loadRecord) + ", "
                   + std::to_string(idleRecord));
    }
 
    localRunTime += loadRecord;
    localIdleTime += idleRecord;
  }
 
  const MInt timeStep = globalTimeStep;
  if(!adaptation) { // Skip storing timings for call from adaptation
    // Runtime and idle-time of last step
    const MFloat localRunTimeLastStep = localRunTime - m_dlbPreviousLocalRunTime;
    const MFloat localIdleTimeLastStep = localIdleTime - m_dlbPreviousLocalIdleTime;
 
    // Store current timings for next evaluation
    m_dlbPreviousLocalRunTime = localRunTime;
    m_dlbPreviousLocalIdleTime = localIdleTime;
 
    // Store timing for imbalance evaluation
    m_dlbTimings.push_back(localRunTimeLastStep);
 
    // Store runtime and idle time for performance evaluations
    if(m_outputDlbTimings) {
      m_dlbTimeStepsAll.push_back(timeStep);
      m_dlbRunTimeAll.push_back(localRunTimeLastStep);
      m_dlbIdleTimeAll.push_back(localIdleTimeLastStep);
 
      // Write timings to disk at final time step
      if(finalTimeStep) {
        storeTimings();
      }
    }
  }
 
  // Check if this is a timer reset time step
  MBool resetTimeStep = (m_lastLoadBalancingTimeStep + m_loadBalancingTimerStartOffset == timeStep
                         || timeStep == m_loadBalancingStartTimeStep + m_loadBalancingTimerStartOffset);
 
  // TODO labels:DLB @ansgar FIXME
  if(m_balanceAfterAdaptation) {
    // resetTimeStep = ((timeStep - m_loadBalancingTimerStartOffset) % m_loadBalancingInterval == 0
    resetTimeStep = resetTimeStep || (timeStep - m_lastAdaptationTimeStep == m_loadBalancingTimerStartOffset);
    // Reset timer at adaptation time step (if this is no DLB step, see below)
    resetTimeStep = resetTimeStep || (timeStep - m_lastAdaptationTimeStep == 0);
  }
 
  MBool dlbTimeStep = ((timeStep - m_loadBalancingOffset) % m_loadBalancingInterval == 0);
  if(m_balanceAfterAdaptation) {
    dlbTimeStep = dlbTimeStep
                  || (((timeStep - m_lastAdaptationTimeStep) == m_loadBalancingOffset)
                      && (m_nAdaptationsSinceBalance >= m_balanceAdaptationInterval || m_dlbStep < 2)
                      && (m_loadBalancingStopTimeStep < 0 || globalTimeStep < m_loadBalancingStopTimeStep));
  }
 
  if(resetTimeStep && !dlbTimeStep) {
    // reset all timer
    resetAllTimer();
 
    // ... and return
    RECORD_TIMER_STOP(m_timers[Timers::DLB]);
    return false;
  }
 
  // Quit early if balancing is not forced and interval does not match
  if(!force && !finalTimeStep && (m_loadBalancingInterval <= 0 || !dlbTimeStep)) {
    RECORD_TIMER_STOP(m_timers[Timers::DLB]);
    return false;
  }
 
  // const MBool atDlbInterval = (timeStep % m_loadBalancingInterval == 0);
  const MBool initialBalance = (timeStep == -1 && force);
  const MBool initialAdaptationOnly = (timeStep == -1 && m_performanceOutput);
  const MBool afterLastDlbStep = (timeStep > m_loadBalancingStopTimeStep && m_loadBalancingStopTimeStep != -1);
  const MBool beforeFirstDlbStep = (timeStep <= m_loadBalancingStartTimeStep);
  // Performance output step without load balancing
  const MBool performanceOutput =
      (m_performanceOutput || finalTimeStep || afterLastDlbStep || beforeFirstDlbStep) && !initialBalance;
 
  // Return if this is not a DLB time step and not the initial balancing
  if(!dlbTimeStep && (!initialBalance || initialAdaptationOnly)) {
    // Continue if forced or at final time step to determine imbalance before returning
    if((!force && !finalTimeStep) || initialAdaptationOnly) {
      RECORD_TIMER_STOP(m_timers[Timers::DLB]);
      return false;
    }
  }
 
  if(!performanceOutput) {
    m_log << "Dynamic load balancing at time step " << timeStep << std::endl;
    cerr0 << "=== Dynamic load balancing at time step " << timeStep << std::endl;
  }
 
  logTimerStatistics("before balance");
 
  // Determine if this is the last DLB step
  const MBool lastDlbStep =
      (timeStep + m_loadBalancingInterval > m_loadBalancingStopTimeStep && m_loadBalancingStopTimeStep != -1)
      && !m_balanceAfterAdaptation;
  // Determine if this is a DLB step at which the partitioning should be reverted to the best so far
  // to improve it by using local shifts only
  const MInt isDlbRevertStep =
      (m_loadBalancingMode == 1 && !lastDlbStep && m_loadBalancingStopTimeStep != -1 && m_dlbNoFinalLocalShifts > 0
       && timeStep == m_loadBalancingStopTimeStep - m_loadBalancingInterval * (m_dlbNoFinalLocalShifts + 1));
 
 
  // Determine the number of timesteps that are timed
  // MInt noTimeSteps = timeStep - m_lastLoadBalancingTimeStep - m_loadBalancingTimerStartOffset;
  MInt noTimeSteps = timeStep - m_dlbLastResetTimeStep;
  if(timeStep < 0) noTimeSteps = 0;
 
  if(m_loadBalancingMode == 1) {
    ASSERT(noTimeSteps > 0, "ERROR: noTimeSteps = " + std::to_string(noTimeSteps));
  }
 
  if(noTimeSteps != (MInt)m_dlbTimings.size()) {
    m_log << "DLB: Number of timings does not match: " + std::to_string(m_dlbTimings.size()) + " "
                 + std::to_string(noTimeSteps) + " (beforeFirstStep = " + std::to_string(beforeFirstDlbStep)
                 + "; afterLastDlbStep = " + std::to_string(afterLastDlbStep) + ")"
          << std::endl;
    std::cerr << "WARNING: number of timings mismatch! " << noTimeSteps << " " << m_dlbTimings.size() << std::endl;
    noTimeSteps = m_dlbTimings.size();
  }
 
 
  std::vector<MFloat> timings(m_dlbTimings.begin(), m_dlbTimings.end());
  // Sort timings for evaluation
  std::sort(timings.begin(), timings.end());
  const MInt noSamples = timings.size();
  if(noSamples < 1 && !initialBalance) {
    TERMM(1, "ERROR: number of samples is < 1");
  }
 
  const MBool truncatedMean = !m_balanceAfterAdaptation;
  // Compute 25% truncated mean for computing loads for static grids
  const MFloat trim = 0.25;
  const MInt lowerBound = std::floor(trim * noSamples);
  const MInt upperBound = std::max(MInt(std::floor((1 - trim) * noSamples)), std::min(noSamples, 1));
  const MInt noSamplesTruncated = upperBound - lowerBound;
 
  if(m_loadBalancingMode == 1 && noSamplesTruncated < 1 && !initialBalance) {
    std::cerr << "ERROR no samples truncated " << noSamplesTruncated << std::endl;
    /* return 0; */
  }
 
  const MFloat truncatedMeanRunTime =
      std::accumulate(&timings[lowerBound], &timings[upperBound], 0.0) / noSamplesTruncated;
  const MFloat meanRunTime = std::accumulate(timings.begin(), timings.end(), 0.0) / noSamples;
  const MFloat meanRunTimeTrigger = (truncatedMean) ? truncatedMeanRunTime : meanRunTime;
 
  /* const MFloat maxLocalRunTime = *std::max_element(timings.begin(), timings.end()); */
  /* const MFloat minLocalRunTime = *std::min_element(timings.begin(), timings.end()); */
 
  // Compute average time per step (all timings)
  MFloat timePerStep = (localRunTime + localIdleTime) / noTimeSteps;
  MFloat maxRunTime = localRunTime / noTimeSteps;
  const MFloat localTimePerStep = timePerStep;
  // Communicate and take maximum value -> assure same value on all domains
  MPI_Allreduce(MPI_IN_PLACE, &timePerStep, 1, MPI_DOUBLE, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "timePerStep");
  MPI_Allreduce(MPI_IN_PLACE, &maxRunTime, 1, MPI_DOUBLE, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "maxRunTime");
 
  MFloat maxIdleTime = localIdleTime / noTimeSteps;
  MFloat minIdleTime = localIdleTime / noTimeSteps;
  MPI_Allreduce(MPI_IN_PLACE, &maxIdleTime, 1, MPI_DOUBLE, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "maxIdleTime");
  MPI_Allreduce(MPI_IN_PLACE, &minIdleTime, 1, MPI_DOUBLE, MPI_MIN, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "minIdleTime");
 
  const MFloat minIdleTimeRel = minIdleTime / timePerStep;
 
  m_log << timeStep << " * Average time per step: " << timePerStep << " ; local: " << localTimePerStep
        << ", idle/comp = " << localIdleTime / localRunTime
        << ", idle/timePerStep = " << localIdleTime / (noTimeSteps * timePerStep) << std::endl;
  m_log << timeStep << " * Relative idle time: max = " << maxIdleTime / timePerStep << ", min "
        << minIdleTime / timePerStep << std::endl;
  m_log << timeStep << " * maxRunTime " << maxRunTime << std::endl;
 
  logDuration(dlbStartTime, "DLB preparation/timers");
 
  const MFloat imbalanceStartTime = wallTime();
  RECORD_TIMER_START(m_timers[Timers::Trigger]);
  MFloat imbalance = -1.0;
  ScratchSpace<MFloat> loads(noDomains(), FUN_, "loads");
  MBool loadBalance = initialBalance;
  if(!loadBalance) { // Skip this for initial balance since there are no timings
    loadBalance = needLoadBalancing(meanRunTimeTrigger, &loads[0], imbalance) || force || (m_forceBalance > 0);
  }
  if(m_forceBalance == 1) {
    m_forceBalance = 0;
  }
  RECORD_TIMER_STOP(m_timers[Timers::Trigger]);
  logDuration(imbalanceStartTime, "Imbalance evaluation");
 
  // Determine variance of timings to detect performance variations
  // TODO labels:CONTROLLER,DLB exclude ranks with low load?
  MFloat localRunTimeVariance = 0.0;
  for(MInt i = 0; i < noSamples; i++) {
    localRunTimeVariance += POW2(timings[i] - meanRunTime);
  }
  localRunTimeVariance /= noSamples;
  const MFloat localRunTimeStdev = std::sqrt(localRunTimeVariance);
 
  const MFloat performanceVariation = localRunTimeStdev / meanRunTime;
  MFloat perfVarMax[2];
  perfVarMax[0] = performanceVariation;
  // Note: use only ranks above a specific load
  perfVarMax[1] = (loads[domainId()] > 0.85) ? performanceVariation : 0.0;
 
  MPI_Allreduce(MPI_IN_PLACE, &perfVarMax, 2, type_traits<MFloat>::mpiType(), MPI_MAX, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "perfVarMax");
  const MFloat maxPerformanceVariation = perfVarMax[1];
 
  const MBool performanceVariationCheck = (maxPerformanceVariation > m_maxPerformanceVarThreshold && !m_adaptation);
 
  m_log << timeStep << " * DLB: Performance variation: load>0.85 " << maxPerformanceVariation << "; all "
        << perfVarMax[0] << std::endl;
 
  m_log << timeStep << " * Average timings: avgRunTime = " << localRunTime / noTimeSteps
        << ", truncated = " << truncatedMean << " " << truncatedMeanRunTime << ", diff "
        << meanRunTime - truncatedMeanRunTime << ", noTimeSteps = " << noTimeSteps << std::endl;
 
  m_log << timeStep << " * timePerStep = " << timePerStep << "; imbalance = " << imbalance
        << "; maxRunTime = " << maxRunTime << "; minIdleTimeRel = " << minIdleTimeRel << std::endl;
 
  if(m_loadBalancingMode == 1 && performanceVariationCheck && !lastDlbStep && !performanceOutput
     && !m_testDynamicLoadBalancing) {
    /* std::cerr << domainId() << " performance variation " << minLocalRunTime << " " */
    /*           << maxLocalRunTime << " " << localRunTimeVariance << " " << localRunTimeStdev */
    /*           << " " << localRunTimeStdev/meanRunTime << std::endl; */
 
    std::stringstream perfMessage;
    perfMessage << "DLB: Performance variation " << maxPerformanceVariation << ", skip DLB step." << std::endl;
    m_log << perfMessage.str();
    if(domainId() == 0) {
      std::cout << perfMessage.str();
    }
 
    resetAllTimer();
 
    m_lastLoadBalancingTimeStep = timeStep;
    m_nAdaptationsSinceBalance = 0;
    RECORD_TIMER_STOP(m_timers[Timers::DLB]);
    return false;
  }
 
  if(m_loadBalancingMode == 0) {
    m_lastLoadBalancingTimeStep = timeStep;
    m_nAdaptationsSinceBalance = 0;
  }
 
  // Check if the current time per timestep is the best one so far. If so, store the new best
  // time/timestep and the corresponding local domain offset. (Note: includes idle times, also the
  // imbalance should not be too large, accept also better imbalance with slightly increased time
  // per step as the latter might be due e.g. network latencies and slower communication)
  MBool newBestTimePerStep = false;
  if(m_loadBalancingMode == 1 && !performanceOutput && !performanceVariationCheck && !m_testDynamicLoadBalancing
     && ((timePerStep < m_timePerStepTotal && imbalance < 1.05 * m_imbalance)
         || (imbalance < m_imbalance && timePerStep < 1.05 * m_timePerStepTotal))) {
    if(m_optPartitionCellOffsetTotal == -1) {
      m_timePerStepTotal = -1.0;
      m_imbalance = -1.0;
    }
    m_log << timeStep << " * Storing new best domain partitioning: timePerStep = " << timePerStep << " ("
          << m_timePerStepTotal << "); imbalance = " << imbalance << " (" << m_imbalance << "), maxRunTime "
          << maxRunTime << " minIdleTimeRel " << minIdleTimeRel << std::endl;
    newBestTimePerStep = true;
    m_timePerStepTotal = timePerStep;
    m_imbalance = imbalance;
    m_optPartitionCellOffsetTotal = gridb().m_localPartitionCellOffsets[0];
 
    RECORD_TIMER_START(m_timers[Timers::IO]);
    // Store partition file for restarting with the same domain offsets on the
    // same number of domains
    // This will overwrite the partition_n[noDomains].[ext] file every time a
    // better partitioning is found.
    std::stringstream partitionFileName;
    partitionFileName << "partition_n" << noDomains();
    gridb().savePartitionFile(partitionFileName.str(), m_optPartitionCellOffsetTotal);
    RECORD_TIMER_STOP(m_timers[Timers::IO]);
  } else if(m_testDynamicLoadBalancing && finalTimeStep) {
    // Always write partition file at final time step for testing
    std::stringstream partitionFileName;
    partitionFileName << "partition_n" << noDomains();
    gridb().savePartitionFile(partitionFileName.str(), gridb().m_localPartitionCellOffsets[0]);
  }
 
  // Return if nothing needs to be done
  if(!loadBalance && !m_forceLoadBalancing && (isDlbRevertStep == 0)) {
    m_log << " * no load imbalance detected at timestep " << timeStep << "!" << std::endl;
    // Continue if this is the last DLB step and the time per step is worse than
    // the best one, else return here
    if(!(lastDlbStep && !newBestTimePerStep)) {
      RECORD_TIMER_STOP(m_timers[Timers::DLB]);
      return false;
    }
  } else {
    if(!performanceOutput) {
      m_log << " * load imbalance detected at timestep " << timeStep << "!" << std::endl;
    }
  }
 
  // Return here for final time step/or performance output after determining imbalance etc.
  if(performanceOutput) {
    // clear timings such that only the next x timesteps are used for the following performance
    // output
    m_log << "Performance evaluation: clear timings at timestep " << timeStep << "!" << std::endl;
 
    resetAllTimer();
 
    RECORD_TIMER_STOP(m_timers[Timers::DLB]);
    return false;
  }
 
  // Determine imbalance statistics *BEFORE* online restart
  printDomainStatistics("before load balancing");
 
  // reset newMinLeven for balancing, otherwise cells with zero noOffsprings occour!
  MInt backupLevel = -1;
  if(gridb().m_newMinLevel > 0) {
    backupLevel = gridb().m_newMinLevel;
    gridb().m_newMinLevel = -1;
  }
 
  gridb().computeGlobalIds();
  gridb().storeMinLevelCells();
 
  // Cancel any open MPI (receive) requests already opened to allow for interleaved execution, since
  // these might lead to conflicting messages
  for(MInt i = 0; i < noSolvers(); i++) {
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Total]);
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DLB]);
 
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DlbOther]);
    solver(i).cancelMpiRequests();
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DlbOther]);
 
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DLB]);
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Total]);
  }
 
  MBool partitionCellChange = false;
  const MFloat updatePartCellsStartTime = wallTime();
  // If enabled: Check if the partition cells need to be updated/changed, i.e., decrease/increase
  // the partition level shift
  //
  // Note: in case of a change in the partition cells the balancing needs to be performed (even if
  // the determined new partitioning is the same as after the partition cell update)
  // TODO labels:CONTROLLER,DLB last DLB or revert step!
  if(m_dlbUpdatePartitionCells && !lastDlbStep) {
    MInt offspringThreshold = gridb().m_partitionCellOffspringThreshold;
    MFloat workloadThreshold = gridb().m_partitionCellWorkloadThreshold;
 
    // For testing: set different thresholds to force partition cell changes
    if(m_testUpdatePartitionCells && m_dlbStep % 2 == 0) {
      offspringThreshold = m_testUpdatePartCellsOffspringThreshold;
      workloadThreshold = m_testUpdatePartCellsWorkloadThreshold;
    }
 
    // Update partition cells
    partitionCellChange = gridb().updatePartitionCells(offspringThreshold, workloadThreshold);
 
    if(partitionCellChange) {
      // Grid with new partition cells will be written with next restart files
      m_saveGridNewPartitionCells = true;
 
      // Reset best found partitioning since not useful anymore with changed partition cells
      m_optPartitionCellOffsetTotal = -1;
    }
  }
  logDuration(updatePartCellsStartTime, "Update partition cells");
 
  // Storage for new partitioning
  MLongScratchSpace partitionCellOffsets(noDomains() + 1, AT_, "partitionCellOffsets");
  MLongScratchSpace globalIdOffsets(noDomains() + 1, AT_, "globalIdOffsets");
 
  const MFloat partitionStartTime = wallTime();
  // Determine new partitioning based on load balancing mode
  if(m_loadBalancingMode == 0 || timeStep == -1) { // Use mode 0 always for initial balance
    RECORD_TIMER_START(m_timers[Timers::Partition]);
    accumulateCellWeights();
 
    // Compute new domain decomposition
    partition(&partitionCellOffsets[0], &globalIdOffsets[0], false);
 
    // Increase DLB step
    m_dlbStep++;
 
    RECORD_TIMER_STOP(m_timers[Timers::Partition]);
  } else {
    // Static mode (i.e. fixed number of DLB steps)
    RECORD_TIMER_START(m_timers[Timers::Partition]);
 
    MBool newPartition = false;
    // Check if this is the last DLB step
    if((!lastDlbStep || m_testDynamicLoadBalancing) && ((isDlbRevertStep == 0) || newBestTimePerStep)) {
      // Determine new partitioning, i.e. new partition cell offsets and domain offsets
      newPartition = loadBalancingPartition(&loads[0], imbalance, &partitionCellOffsets[0], &globalIdOffsets[0]);
    } else {
      // Check if the current configuration is not the best one (wrt. time per timestep)
      if(((!newBestTimePerStep && !m_adaptation) || (isDlbRevertStep != 0)) && m_optPartitionCellOffsetTotal != -1) {
        m_log << timeStep << " * DLB reverting to best configuration." << std::endl;
        // Last DLB step: since there is still some load imbalance, revert to the best partitioning
        // found
 
        // Gather current partition cell offsets
        MLongScratchSpace localPartitionCellOffsets(noDomains() + 1, AT_, "localPartitionCellOffsets");
        MPI_Allgather(&gridb().m_localPartitionCellOffsets[0], 1, type_traits<MLong>::mpiType(),
                      &localPartitionCellOffsets[0], 1, type_traits<MLong>::mpiType(), gridb().mpiComm(), AT_,
                      "gridb().m_localPartitionCellOffsets[0]", "localPartitionCellOffsets[0]");
        localPartitionCellOffsets[noDomains()] = gridb().m_noPartitionCellsGlobal;
 
        // Gather 'optimal' partition cell offsets
        MPI_Allgather(&m_optPartitionCellOffsetTotal, 1, type_traits<MLong>::mpiType(), &partitionCellOffsets[0], 1,
                      type_traits<MLong>::mpiType(), gridb().mpiComm(), AT_, "m_optPartitionCellOffsetTotal",
                      "partitionCellOffsets[0]");
        partitionCellOffsets[noDomains()] = gridb().m_noPartitionCellsGlobal;
 
        // Check if this is a new partitioning, i.e. different from the old one
        newPartition =
            (std::mismatch(partitionCellOffsets.begin(), partitionCellOffsets.end(), &localPartitionCellOffsets[0]))
                .first
            != partitionCellOffsets.end();
 
        loadBalancingCalcNewGlobalOffsets(&localPartitionCellOffsets[0], &partitionCellOffsets[0], &globalIdOffsets[0]);
      } else {
        // Current partitioning is accepted
        // Note: this will not work if m_forceLoadBalancing is activated, since
        // the number of cells to send, etc. are not determined.
        newPartition = false;
      }
    }
 
    RECORD_TIMER_STOP(m_timers[Timers::Partition]);
 
    // Return here if partitioning did not change
    if(!newPartition && !partitionCellChange && !m_forceLoadBalancing) {
      m_log << "Dynamic load balancing: load imbalance detected but partition did not change at "
               "timestep "
            << timeStep << "!" << std::endl;
      RECORD_TIMER_STOP(m_timers[Timers::DLB]);
      return false;
    }
  }
  logDuration(partitionStartTime, "New partitioning");
 
  RECORD_TIMER_START(m_timers[Timers::Prepare]);
  const MFloat prepareStartTime = wallTime();
 
  // Disable all running DLB timers and store status
  ScratchSpace<MBool> dlbTimersStatus(std::max(noDlbTimers, 1), AT_, "dlbTimersStatus");
  maia::dlb::g_dlbTimerController.disableAllDlbTimers(&dlbTimersStatus[0]);
 
  // Reset solvers before load balancing
  for(MInt i = 0; i < noSolvers(); i++) {
    const MFloat resetStartTime = wallTime();
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Total]);
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DLB]);
 
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Reset]);
    solver(i).resetSolver();
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Reset]);
 
 
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DlbOther]);
    // Communicate global solver variables from current local rank 0 to all domains
    communicateGlobalSolverVars(&solver(i));
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DlbOther]);
 
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DLB]);
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Total]);
    logDuration(resetStartTime, "  Reset solver #" + std::to_string(i));
  }
 
  const MFloat prepareGridStartTime = wallTime();
  gridb().deletePeriodicConnection(false);
 
  std::vector<MInt>().swap(gridb().m_minLevelCells);
 
  gridb().localToGlobalIds();
 
  const MInt oldNoCells = gridb().treeb().size();
 
  MIntScratchSpace noCellsToReceiveByDomain(noDomains() + 1, AT_, "noCellsToReceiveByDomain");
  MIntScratchSpace noCellsToSendByDomain(noDomains() + 1, AT_, "noCellsToSendByDomain");
  // Mapping: gridCellId -> sort index in buffer
  MIntScratchSpace sortedCellId(oldNoCells, AT_, "sortedCellId");
  // Reverse mapping: sort index in buffer -> gridCellId
  MIntScratchSpace bufferIdToCellId(oldNoCells, AT_, "bufferIdToCellId");
  bufferIdToCellId.fill(-1);
 
  loadBalancingCalcNoCellsToSend(&globalIdOffsets[0], &noCellsToSendByDomain[0], &noCellsToReceiveByDomain[0],
                                 &sortedCellId[0], &bufferIdToCellId[0]);
  logDuration(prepareGridStartTime, "  Prepare grid");
 
  std::vector<std::vector<MInt>> dataSendSize{};
  std::vector<std::vector<MInt>> dataRecvSize{};
 
  // Determine communication data sizes for each solver
  for(MInt i = 0; i < noSolvers(); i++) {
    const MFloat dataStartTime = wallTime();
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Total]);
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DLB]);
 
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::CalcDataSizes]);
    determineDataSizesDlb(i, 0, &noCellsToSendByDomain[0], &bufferIdToCellId[0], dataSendSize, dataRecvSize);
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::CalcDataSizes]);
 
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DlbOther]);
    // Change local to global ids here, since local ids previously needed to determine data sizes
    solver(i).localToGlobalIds();
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DlbOther]);
 
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DLB]);
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Total]);
    logDuration(dataStartTime, "  Data sizes solver #" + std::to_string(i));
  }
 
 
  std::vector<std::vector<MInt>> dataSendSizeCoupler{};
  std::vector<std::vector<MInt>> dataRecvSizeCoupler{};
 
  // Determine communication data sizes for each coupler
  for(MInt i = 0; i < noCouplers(); i++) {
    const MFloat dataStartTime = wallTime();
    const MInt timerId = noSolvers() + i;
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Total]);
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::DLB]);
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::CalcDataSizes]);
    determineDataSizesDlb(i, 1, &noCellsToSendByDomain[0], &bufferIdToCellId[0], dataSendSizeCoupler,
                          dataRecvSizeCoupler);
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::CalcDataSizes]);
 
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::DLB]);
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Total]);
    logDuration(dataStartTime, "  Data sizes coupler #" + std::to_string(i));
  }
 
 
  RECORD_TIMER_STOP(m_timers[Timers::Prepare]);
  logDuration(prepareStartTime, "Prepare balancing");
 
  const MFloat balanceGridStartTime = wallTime();
  RECORD_TIMER_START(m_timers[Timers::BalanceGrid]);
  // Rebalance grid
  gridb().balance(&noCellsToReceiveByDomain[0], &noCellsToSendByDomain[0], &sortedCellId[0], &partitionCellOffsets[0],
                  &globalIdOffsets[0]);
  RECORD_TIMER_STOP(m_timers[Timers::BalanceGrid]);
  logDuration(balanceGridStartTime, "Balance grid");
 
  RECORD_TIMER_START(m_timers[Timers::BalanceSolvers]);
  const MFloat balanceSolversStartTime = wallTime();
  // Rebalance solvers
  for(MInt i = 0; i < noSolvers(); i++) {
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Total]);
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DLB]);
 
    RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Balance]);
    const MFloat balanceSolverStartTime = wallTime();
 
    if(!solver(i).hasSplitBalancing()) {
      solver(i).balance(&noCellsToReceiveByDomain[0], &noCellsToSendByDomain[0], &sortedCellId[0], oldNoCells);
    } else {
      // Variables to store pointers to allocated memory
      std::vector<MInt*> intDataRecv{};
      std::vector<MLong*> longDataRecv{};
      std::vector<MFloat*> floatDataRecv{};
      std::vector<MInt> dataTypes{};
 
      // TODO labels:CONTROLLER,DLB,TIMERS create shorthand for: timer start, function call, timer stop, logDuration?
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::Redistribute]);
      // Communicate all solver data
      redistributeDataDlb(i, 0, dataSendSize[i], dataRecvSize[i], &bufferIdToCellId[0], oldNoCells, intDataRecv,
                          longDataRecv, floatDataRecv, dataTypes);
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Redistribute]);
      logDuration(RETURN_TIMER(m_solverTimers[i][SolverTimers::Redistribute]),
                  "Redistribute solver #" + std::to_string(i), true);
 
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::BalancePre]);
      // Perform reinitialization steps prior to setting solver data
      // i.e. apply balance to the proxy!
      solver(i).balancePre();
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::BalancePre]);
      logDuration(RETURN_TIMER(m_solverTimers[i][SolverTimers::BalancePre]), "BalancePre solver #" + std::to_string(i),
                  true);
 
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::SetData]);
      setDataDlb(i, 0, intDataRecv, longDataRecv, floatDataRecv, dataTypes, false);
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::SetData]);
      logDuration(RETURN_TIMER(m_solverTimers[i][SolverTimers::SetData]), "SetData1 solver #" + std::to_string(i),
                  true);
 
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::DlbOther]);
      // Change global to local ids in the solver
      solver(i).globalToLocalIds();
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DlbOther]);
 
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::BalancePost]);
      // Perform reinitialization steps after setting solver data the first time
      solver(i).balancePost();
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::BalancePost]);
      logDuration(RETURN_TIMER(m_solverTimers[i][SolverTimers::BalancePost]),
                  "BalancePost solver #" + std::to_string(i), true);
 
      RECORD_TIMER_START(m_solverTimers[i][SolverTimers::SetData]);
      // Set solver data again (if required by solver) and deallocate buffers
      setDataDlb(i, 0, intDataRecv, longDataRecv, floatDataRecv, dataTypes, true);
      RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::SetData]);
      logDuration(RETURN_TIMER(m_solverTimers[i][SolverTimers::SetData]), "SetData2 solver #" + std::to_string(i),
                  true);
 
      intDataRecv.clear();
      longDataRecv.clear();
      floatDataRecv.clear();
    }
 
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Balance]);
    logDuration(balanceSolverStartTime, "Balance solver #" + std::to_string(i));
 
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::DLB]);
    RECORD_TIMER_STOP(m_solverTimers[i][SolverTimers::Total]);
  }
  RECORD_TIMER_STOP(m_timers[Timers::BalanceSolvers]);
  logDuration(balanceSolversStartTime, "Balance solvers");
 
  RECORD_TIMER_START(m_timers[Timers::BalanceCouplers]);
  const MFloat balanceCouplersStartTime = wallTime();
  // Rebalance couplers (if required)
  for(MInt i = 0; i < noCouplers(); i++) {
    const MInt timerId = noSolvers() + i;
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Total]);
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::DLB]);
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Balance]);
 
    // Variables to store pointers to allocated memory
    std::vector<MInt*> intDataRecv{};
    std::vector<MLong*> longDataRecv{};
    std::vector<MFloat*> floatDataRecv{};
    std::vector<MInt> dataTypes{};
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::Redistribute]);
    // Communicate all coupler data
    redistributeDataDlb(i, 1, dataSendSizeCoupler[i], dataRecvSizeCoupler[i], &bufferIdToCellId[0], oldNoCells,
                        intDataRecv, longDataRecv, floatDataRecv, dataTypes);
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Redistribute]);
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::BalancePre]);
    // Perform reinitialization steps prior to setting solver data
    coupler(i).balancePre();
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::BalancePre]);
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::SetData]);
    setDataDlb(i, 1, intDataRecv, longDataRecv, floatDataRecv, dataTypes, false);
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::SetData]);
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::BalancePost]);
    // Perform reinitialization steps after setting solver data the first time
    coupler(i).balancePost();
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::BalancePost]);
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::SetData]);
    // Set coupler data again (if required by coupler) and deallocate buffers
    setDataDlb(i, 1, intDataRecv, longDataRecv, floatDataRecv, dataTypes, true);
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::SetData]);
 
    intDataRecv.clear();
    longDataRecv.clear();
    floatDataRecv.clear();
 
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Balance]);
 
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::DLB]);
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::Total]);
  }
  RECORD_TIMER_STOP(m_timers[Timers::BalanceCouplers]);
  logDuration(balanceCouplersStartTime, "Balance couplers");
 
  const MFloat finalizeBalanceStartTime = wallTime();
  // finalize balance for solvers and couplers
  for(MInt i = 0; i < noSolvers(); i++) {
    const MFloat finalizeSolverStartTime = wallTime();
    solver(i).finalizeBalance();
    for(MInt j = 0; j < noCouplers(); j++) {
      const MFloat finalizeCouplerStartTime = wallTime();
      coupler(j).finalizeBalance(i);
      logDuration(finalizeCouplerStartTime,
                  "Finalize balance #" + std::to_string(i) + " coupler #" + std::to_string(j));
    }
    logDuration(finalizeSolverStartTime, "Finalize balance solver #" + std::to_string(i));
  }
  logDuration(finalizeBalanceStartTime, "Finalize balance total");
 
  // Determine imbalance statistics *AFTER* online restart
  printDomainStatistics("after load balancing");
 
  // Enable all previously active DLB timers again
  maia::dlb::g_dlbTimerController.enableAllDlbTimers(&dlbTimersStatus[0]);
 
  // Reset DLB timer records
  maia::dlb::g_dlbTimerController.resetRecords();
 
  m_dlbPreviousLocalRunTime = 0.0;
  m_dlbPreviousLocalIdleTime = 0.0;
 
  // ... reset timings ...
  resetAllTimer();
 
  m_lastLoadBalancingTimeStep = timeStep;
  m_nAdaptationsSinceBalance = 0;
 
  // Write timings to file
  if(m_outputDlbTimings) {
    storeTimings();
  }
 
  logDuration(dlbStartTime, "Balance total");
  logDurations(durations, "DLB", gridb().mpiComm(), globalDomainId(), globalNoDomains());
 
  const MFloat dlbTimeTotal = wallTime() - dlbStartTime;
  std::stringstream dlbMessage;
  dlbMessage << "=== Dynamic load balancing performed at timestep " << timeStep << "! Duration: " << dlbTimeTotal
             << " s" << std::endl;
  m_log << dlbMessage.str();
  cerr0 << dlbMessage.str();
 
  printAllocatedMemory(oldAllocatedBytes, "gridcontroller::balance()", gridb().mpiComm());
 
  // restore backup value for next restartFile!
  if(backupLevel > 0) {
    gridb().m_newMinLevel = backupLevel;
  }
 
  RECORD_TIMER_STOP(m_timers[Timers::DLB]);
  return true;
}

◆ castAdaptationIntervalToMultipleOfCoarsestTimeStep()

template<MInt nDim>

void maia::grid::Controller< nDim >::castAdaptationIntervalToMultipleOfCoarsestTimeStep	(	MInt	maxLevel,
		MInt	maxUniformRefinementLevel
	)

Definition at line 3302 of file cartesiangridcontroller.h.

                                                                                                          {
  this->m_adaptationInterval = this->m_adaptationInterval / IPOW2(maxLevel - maxUniformRefinementLevel)
                               * IPOW2(maxLevel - maxUniformRefinementLevel);
  this->m_adaptationStart = this->m_adaptationStart / IPOW2(maxLevel - maxUniformRefinementLevel)
                            * IPOW2(maxLevel - maxUniformRefinementLevel);
  std::cout << "Set adaptationStart to: " << this->m_adaptationStart
            << " and adaptationInterval to: " << this->m_adaptationInterval << "\n";
}

◆ cellOutsideVec()

template<MInt nDim>

std::vector< std::function< MInt(const MFloat *, const MInt, const MInt)> > maia::grid::Controller< nDim >::cellOutsideVec ( )

inlineprivate

Definition at line 366 of file cartesiangridcontroller.h.

                                                                                       {
    std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)>> vec(noSolvers());
    for(MInt i = 0; i < noSolvers(); i++) {
      vec[i] = std::bind(&Solver::cellOutside, &solver(i), std::placeholders::_1, std::placeholders::_2,
                         std::placeholders::_3);
    }
    return vec;
  }

◆ communicateGlobalSolverVars()

template<MInt nDim>

void maia::grid::Controller< nDim >::communicateGlobalSolverVars ( Solver *const solver )

private

Communicate the solver variables which should be the same on all ranks from the solver local rank 0 to all other domains. This is needed if during load balancing the ranks of a solver change, i.e. ranks becoming active/inactive due to partitioning changes

Definition at line 2854 of file cartesiangridcontroller.h.

                                                                       {
  TRACE();
 
  // Determine solver local root domain
  const MInt localRootGlobalDomain = solverLocalRootDomain(solver);
 
  std::vector<MInt> globalIdVars(0);
  std::vector<MFloat> globalFloatVars(0);
 
  // Request global solver variables
  solver->getGlobalSolverVars(globalFloatVars, globalIdVars);
 
  const MInt noIdVars = globalIdVars.size();
  const MInt noFloatVars = globalFloatVars.size();
 
  // Distribute all global variables from local rank 0 to all domains
  MPI_Bcast(&globalIdVars[0], noIdVars, type_traits<MInt>::mpiType(), localRootGlobalDomain, gridb().mpiComm(), AT_,
            "globalIdVars[0]");
  MPI_Bcast(&globalFloatVars[0], noFloatVars, type_traits<MFloat>::mpiType(), localRootGlobalDomain, gridb().mpiComm(),
            AT_, "globalFloatVars[0]");
 
  // Set variables in solver
  solver->setGlobalSolverVars(globalFloatVars, globalIdVars);
}

◆ computeWeights()

template<MInt nDim>

void maia::grid::Controller< nDim >::computeWeights	(	const MFloat *	loads,
		const MFloat	domainWeight,
		std::vector< MFloat > &	weights
	)

private

Author: Ansgar Niemoeller (ansgar) a.nie.nosp@m.moel.nosp@m.ler@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Parameters

[in]	loads	Domain loads.
[in]	domainWeight	Local domain weight or processing capacity.
[out]	weights	Vector of computed weights.

Definition at line 40 of file cartesiangridcontroller.cpp.

                                                                              {
  TRACE();
  using namespace maia;
 
  // Number of load quantities (should be the same for every rank)
  const MInt noLoadTypes = globalNoLoadTypes();
  if(noLoadTypes < 1) {
    TERMM(1,
          "There are no load quantities present in the solver(s), but needed to determine a new "
          "partition for dynamic load balancing!");
  }
  weights.resize(noLoadTypes);
  std::fill(weights.begin(), weights.end(), 1.0);
 
  MIntScratchSpace loadQuantities(noLoadTypes, AT_, "loadQuantities");
  MFloatScratchSpace loadQuantitiesScaled(noLoadTypes, AT_, "loadQuantitiesScaled");
  MFloatScratchSpace globalLoadQuantities(noDomains(), noLoadTypes, AT_, "globalLoadQuantities");
 
  // Get the solver specific load quantities
  getLoadQuantities(&loadQuantities[0]);
 
  // Scale the load quantities with the domain weight
  for(MInt i = 0; i < noLoadTypes; i++) {
    loadQuantitiesScaled[i] = domainWeight * loadQuantities[i];
  }
 
  // Gather (scaled) load quantities on rank 0
  MPI_Gather(&loadQuantitiesScaled[0], noLoadTypes, type_traits<MFloat>::mpiType(), &globalLoadQuantities[0],
             noLoadTypes, type_traits<MFloat>::mpiType(), 0, gridb().mpiComm(), AT_, "loadQuantitiesScaled[0]",
             "globalLoadQuantities[0]");
 
  if(domainId() == 0) {
    // Estimate the parameters, i.e. the new weights
    estimateParameters(noDomains(), noLoadTypes, &globalLoadQuantities[0], &loads[0], &weights[0]);
    // limit weights for certain types
    if(!m_testDynamicLoadBalancing) {
      MInt offset = 0;
      for(MInt i = 0; i < noSolvers(); i++) {
        solver(i).limitWeights(&weights[offset]);
        offset += solver(i).noLoadTypes();
      }
      // scale weights
      const MFloat norm = std::accumulate(&weights[0], &weights[0] + noLoadTypes, 0.0);
      for(MInt i = 0; i < noLoadTypes; i++) {
        weights[i] /= norm;
      }
    }
  }
 
  // Broadcast the estimated parameters/weights
  MPI_Bcast(&weights[0], noLoadTypes, type_traits<MFloat>::mpiType(), 0, gridb().mpiComm(), AT_, "weights[0]");
 
  /* storeLoadsAndWeights(&loads[0], noLoadTypes, &loadQuantities[0], domainWeight, &weights[0]); */
 
  { // Output to m_log
    m_log << " * computed weights:";
    MInt count = 0;
    for(MInt i = 0; i < noSolvers(); i++) {
      const MInt solverCount = solver(i).noLoadTypes();
      std::vector<MString> names(solverCount);
      std::vector<MFloat> weightsTmp(solverCount);
 
      solver(i).getDefaultWeights(&weightsTmp[0], names);
      for(MInt j = 0; j < solverCount; j++) {
        m_log << " s" << i << "_" << names[j] << ":" << count;
        count++;
      }
    }
 
    m_log << std::endl << " * computed weights:";
    for(MInt i = 0; i < noLoadTypes; i++) {
      TERMM_IF_COND(std::isnan(weights[i]), "computed weight is NaN");
      m_log << " " << std::scientific << weights[i];
    }
    m_log << std::endl;
  }
 
  // use static pre-defined weights instead
  if(m_dlbStaticWeights != nullptr && m_dlbStaticWeightMode > -1) {
    if(m_dlbStaticWeightMode == 0
       || (m_dlbStaticWeightMode > 0 && m_dlbStep < m_dlbStaticWeightMode && m_dlbStaticWeightMode != 98)
       || (m_dlbStaticWeightMode == 98 && m_dlbStep % 2 == 0)) {
      m_log << "Using static weights at dlb-step " << m_dlbStep << " :" << std::endl;
      for(MInt i = 0; i < noLoadTypes; i++) {
        weights[i] = m_dlbStaticWeights[i];
        m_log << " " << std::scientific << weights[i];
      }
      m_log << std::endl;
    }
    if(m_dlbStaticWeightMode > 0 && m_dlbStaticWeightMode < 90 && m_dlbStep == m_dlbStaticWeightMode) {
      m_dlbStaticWeightMode = -1;
      m_log << "Switching to computed weights from now on!" << std::endl;
    }
  }
 
 
  // store weights and use last weights at restart
  if(m_dlbLastWeights == nullptr) {
    mAlloc(m_dlbLastWeights, noLoadTypes, "dlbLastWeights", AT_);
  }
  for(MInt i = 0; i < noLoadTypes; i++) {
    m_dlbLastWeights[i] = weights[i];
  }
}

◆ coupler() [1/2]

template<MInt nDim>

Coupling & maia::grid::Controller< nDim >::coupler ( const MInt couplerId )

inlineprivate

Definition at line 144 of file cartesiangridcontroller.h.

144{ return *m_couplers->at(couplerId); }

maia::grid::Controller::m_couplers

const std::vector< std::unique_ptr< Coupling > > * m_couplers

Definition: cartesiangridcontroller.h:210

◆ coupler() [2/2]

template<MInt nDim>

const Coupling & maia::grid::Controller< nDim >::coupler ( const MInt couplerId ) const

inlineprivate

Definition at line 145 of file cartesiangridcontroller.h.

145{ return *m_couplers->at(couplerId); }

◆ determineDataSizesDlb()

template<MInt nDim>

void maia::grid::Controller< nDim >::determineDataSizesDlb	(	const MInt	solverId,
		const MInt	mode,
		const MInt *const	noCellsToSend,
		const MInt *const	bufferIdToCellId,
		std::vector< std::vector< MInt > > &	sendSizeVector,
		std::vector< std::vector< MInt > > &	recvSizeVector
	)

private

Definition at line 2882 of file cartesiangridcontroller.h.

                                                                                         {
  TRACE();
 
  const MBool isSolver = (mode == 0); // mode 0: solver; mode 1: coupler
  const MInt timerId = (isSolver) ? id : noSolvers() + id;
  // Get the number of quantities that need to be communicated
  MInt dataCount = (isSolver) ? solver(id).noCellDataDlb() : coupler(id).noCellDataDlb();
 
  RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::CalcDataSizesMpiBlocking]);
  // Note: since inactive ranks might not return the correct count, take the max among all domains
  MPI_Allreduce(MPI_IN_PLACE, &dataCount, 1, type_traits<MInt>::mpiType(), MPI_MAX, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "dataCount");
  RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::CalcDataSizesMpiBlocking]);
 
  std::vector<MInt> dataSendSize(dataCount * (globalNoDomains() + 1));
  std::vector<MInt> dataRecvSize(dataCount * (globalNoDomains() + 1));
 
  // Determine data size to send/receive to/from each domain for all quantities
  // last entry holds the sum
  for(MInt dataId = 0; dataId < dataCount; dataId++) {
    const MInt offset = dataId * (globalNoDomains() + 1);
    MInt bufferId = 0;
 
    // Loop over all domains and determine the data send size for this domain
    for(MInt domain = 0; domain < globalNoDomains(); domain++) {
      MInt domainDataSize = 0;
 
      // Loop over all grid cells that are send to this domain
      for(MInt i = 0; i < noCellsToSend[domain]; i++) {
        // Inquire the data size for this cell
        // TODO labels:DLB @ansgar inquire all data sizes at once?
        const MInt cellId = bufferIdToCellId[bufferId];
        ASSERT(cellId > -1, "");
        const MInt dataSize =
            (isSolver) ? solver(id).cellDataSizeDlb(dataId, cellId) : coupler(id).cellDataSizeDlb(dataId, cellId);
        domainDataSize += dataSize;
        bufferId++;
      }
 
      // Store total data size to send to this domain for this quantity
      dataSendSize[offset + domain] = domainDataSize;
    }
 
    // TODO labels:DLB @ansgar perform just one alltoall for all data ids?
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::CalcDataSizesMpi]);
    // Exchange data send size with all other domains and store as receive size
    MPI_Alltoall(&dataSendSize[offset], 1, type_traits<MInt>::mpiType(), &dataRecvSize[offset], 1,
                 type_traits<MInt>::mpiType(), gridb().mpiComm(), AT_, "dataSendSize[offset]", "dataRecvSize[offset]");
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::CalcDataSizesMpi]);
 
    // Compute the total data send size and store as last entry
    dataSendSize[offset + globalNoDomains()] =
        std::accumulate(&dataSendSize[offset], &dataSendSize[offset + globalNoDomains()], 0);
 
    // Compute the total data receive size and store as last entry
    dataRecvSize[offset + globalNoDomains()] =
        std::accumulate(&dataRecvSize[offset], &dataRecvSize[offset + globalNoDomains()], 0);
  }
 
  sendSizeVector.push_back(dataSendSize);
  recvSizeVector.push_back(dataRecvSize);
}

◆ domainId()

template<MInt nDim>

MInt maia::grid::Controller< nDim >::domainId ( )

inlineprivate

Definition at line 134 of file cartesiangridcontroller.h.

134{ return gridb().domainId(); }

CartesianGrid::domainId

MInt domainId() const

Return the domainId (rank)

Definition: cartesiangrid.h:1057

◆ estimateParameters()

template<MInt nDim>

void maia::grid::Controller< nDim >::estimateParameters	(	MInt	m,
		MInt	n,
		const MFloat *const	A,
		const MFloat *const	b,
		MFloat *const	x
	)

private

The linear least squares problem A*x=b is regularized with a Tikhonov regularization and solved using Eigen. Parameters x are assumed to be positive. The regularization parameter is increased as long as the residual decreases. Estimated parameters are normalized in the 1-norm and any negative parameter is set to zero.

Parameters

[in]	m	Number of rows of the matrix A.
[in]	n	Number of columns of the matrix A (x=1 if n==1).
[in]	A	Matrix (size m x n).
[in]	b	Right hand side vector (size m).
[out]	x	Estimated parameter vector (size n).

Definition at line 159 of file cartesiangridcontroller.cpp.

                                                                       {
  TRACE();
 
  // Always return the same parameters for testing
  if(m_testDynamicLoadBalancing) {
    // For testing: switch weights in consecutive steps to force a partitioning change
    std::cerr << "estimateParameters: n = " << n << ", setting weights to 1.0" << std::endl;
    std::fill_n(x, n, 1.0);
    if(n > 1 && m_dlbStep % 2 == 0) {
      if(Context::propertyExists("solverWeights_0")) {
        std::vector<MFloat> weights;
        getSpecifiedSolverWeights(weights);
        std::copy(weights.begin(), weights.end(), &x[0]);
        std::cerr << "estimateParameters: using specified solver weights" << std::endl;
      } else {
        std::cerr << "estimateParameters: n = " << n << ", setting w[n-1] = 5.0" << std::endl;
        x[n - 1] = 5.0;
      }
    }
    return;
  }
 
  if(n == 1) {
    // Single parameter, return 1
    x[0] = 1.0;
  } else {
    // Solve the (overdetermined) system of linear equations (linear least squares)
 
    // Regularization parameter
    MFloat alpha = 0.01;
    // Maximum number of iterations
    const MInt maxNoIt = 30;
 
    // Determine column sums
    MFloatScratchSpace sum_j(n, AT_, "sum_j");
    for(MInt j = 0; j < n; j++) {
      sum_j[j] = 0.0;
      for(MInt i = 0; i < m; i++) {
        sum_j[j] += A[i * n + j];
      }
      if(m_debugBalance) {
        m_log << " * estimateParameters regularize " << j << " proportional to " << sum_j[j] << std::endl;
      }
    }
 
    MFloat oldResidual = 0.0;
    MFloatScratchSpace oldParam(n, AT_, "oldParam");
 
    MInt it = 0;
    while(it < maxNoIt) {
      // Work on copies of A and b since they get overwritten and additional
      // rows are added to allow a regularization of the overdetermined system
      MFloatScratchSpace A_work((m + n) * n, AT_, "A_work");
      MFloatScratchSpace b_work(m + n, AT_, "b_work");
 
      // Copy matrix A and fill additional rows with zeros
      std::copy_n(&A[0], m * n, &A_work[0]);
      std::fill_n(&A_work[m * n], n * n, 0.0);
 
      // Set 'alpha' on diagonal of additional (n x n) solver (Tikhonov regular.)
      for(MInt i = 0; i < n; i++) {
        // Set proportional to average column value or to 1 if column sum is zero
        const MFloat reg = (sum_j[i] > 0) ? alpha * sum_j[i] / noDomains() : 1.0;
        A_work[m * n + i * (n + 1)] = reg;
      }
 
      // Copy vector b and add zeros for regularization
      std::copy_n(&b[0], m, &b_work[0]);
      std::fill_n(&b_work[m], n, 0.0);
 
      // Compute the minimum-norm solution to the linear least squares problem:
      // minimize 2-norm(|b - A*x|) using singular value decomposition (SVD) of
      // A. A is an M-by-N matrix which can be rank deficient.
 
      // TODO labels:CONTROLLER,totest check if this works properly!
 
      std::vector<MFloat> singularValues = maia::math::svd(A_work.data(), b_work.data(), m + n, n, x);
 
      // Compute condition number of A (cond_2-norm = S(0)/S(min(m,n)-1))
      // TODO labels:CONTROLLER this is not correct it should be max(S)/min(S)
      const MFloat cond_2 = singularValues[0] / singularValues[std::min(m, n) - 1];
 
      // Compute residual
      // TODO labels:CONTROLLER could be replaced by
      // maia::math::MFloatVectorXd residual = A_eigen*x_eigen - b_eigen
      MFloat residual = 0.0;
      MFloat maxDiff = 0.0;
      MFloatScratchSpace mxv(m, AT_, "mxv");
      for(MInt i = 0; i < m; i++) {
        mxv[i] = 0.0;
        for(MInt j = 0; j < n; j++) {
          mxv[i] += A[i * n + j] * x[j];
        }
        residual += POW2(mxv[i] - b[i]);
        maxDiff = std::max(maxDiff, mxv[i] - b[i]);
      }
 
      if(m_debugBalance) {
        m_log << " * estimateParameters iteration " << it << ":";
        for(MInt i = 0; i < n; i++) {
          m_log << std::scientific << " " << x[i];
        }
        m_log << "; residual " << residual << "; maxDiff " << maxDiff << std::endl;
      }
 
      /* if (std::any_of(&x[0], &x[0] + n, [](MFloat p) { return p < 0.0; }) || cond_2 > maxCond)
       * { */
      /* const MBool hasNegativeParam = std::any_of(&x[0], &x[0] + n, [](MFloat p) { return p <
       * 0.0; }); */
 
      // If the residual decreased increase regularization constant and repeat
      if(it == 0 || residual <= oldResidual) {
        alpha *= 1.5;
        it++;
 
        oldResidual = residual;
        std::copy_n(&x[0], n, &oldParam[0]);
 
        if(m_debugBalance) {
          if(it == maxNoIt) {
            m_log << " * estimateParameters max iterations reached " << it << ", alpha " << alpha
                  << ", condition number " << cond_2 << std::endl;
          } else {
            m_log << " * estimateParameters increase alpha to " << alpha << ", condition number " << cond_2
                  << ", residual " << residual << std::endl;
          }
        }
      } else {
        // Take values of previous iteration
        std::copy_n(&oldParam[0], n, &x[0]);
 
        if(m_debugBalance) {
          m_log << " * estimateParameters found solution! iteration " << it << ", condition number " << cond_2
                << ", residual " << residual << ", previous " << oldResidual << ", maxDiff " << maxDiff;
          for(MInt i = 0; i < n; i++) {
            m_log << std::scientific << " " << x[i];
          }
          m_log << std::endl;
        }
        break;
      }
    } // alpha loop
  }
 
  // Prevent negative weights and Nans, i.e. set to zero
  for(MInt i = 0; i < n; i++) {
    if(x[i] < 0.0 || std::isnan(x[i])) {
      x[i] = 0.0;
    }
  }
  // Renormalize parameter vector using 1-norm, such that parameters are better suited for
  // comparison and will add up to one.
  const MFloat norm = std::accumulate(&x[0], &x[0] + n, 0.0);
  for(MInt i = 0; i < n; i++) {
    x[i] /= norm;
  }
}

◆ getLoadQuantities()

template<MInt nDim>

void maia::grid::Controller< nDim >::getLoadQuantities ( MInt *const loadQuantities )

private

Definition at line 2499 of file cartesiangridcontroller.h.

                                                                   {
  // Fill with zeros since inactive solvers will not write anything
  std::fill_n(&loadQuantities[0], globalNoLoadTypes(), 0);
 
  MInt offset = 0;
  for(MInt i = 0; i < noSolvers(); i++) {
    solver(i).getLoadQuantities(&loadQuantities[offset]);
    offset += solver(i).noLoadTypes();
  }
}

◆ getSpecifiedSolverWeights()

template<MInt nDim>

void maia::grid::Controller< nDim >::getSpecifiedSolverWeights ( std::vector< MFloat > & weights )

private

Definition at line 2355 of file cartesiangridcontroller.h.

                                                                           {
  TRACE();
 
  // Determine the weights for all solvers
  const MInt noLoadTypes = globalNoLoadTypes();
  weights.resize(noLoadTypes);
  std::fill(weights.begin(), weights.end(), 0.0);
 
  MInt offset = 0;
  for(MInt i = 0; i < noSolvers(); i++) {
    const MInt solverCount = solver(i).noLoadTypes();
    std::vector<MString> names(solverCount);
 
    // Get default weights (+names) for this solver
    solver(i).getDefaultWeights(&weights[offset], names);
 
    // Check for user specified solver weights
    const MString propName = "solverWeights_" + std::to_string(i);
    if(Context::propertyExists(propName)) {
      if(Context::propertyLength(propName) != solverCount) {
        TERMM(1, "wrong length of property '" + propName + "', should be of length " + std::to_string(solverCount));
      }
      for(MInt j = 0; j < solverCount; j++) {
        weights[offset + j] = Context::getBasicProperty<MFloat>(propName, AT_, j);
      }
    }
 
    for(MInt j = 0; j < solverCount; j++) {
      m_log << "Solver #" << i << ", weight #" << j << ": " << names[j] << ", " << weights[offset + j] << std::endl;
    }
 
    offset += solverCount;
  }
}

◆ globalNoLoadTypes()

template<MInt nDim>

MInt maia::grid::Controller< nDim >::globalNoLoadTypes

private

Definition at line 2488 of file cartesiangridcontroller.h.

                                         {
  MInt noLoadTypes = 0;
  for(MInt i = 0; i < noSolvers(); i++) {
    noLoadTypes += solver(i).noLoadTypes();
  }
  return noLoadTypes;
}

◆ gridb() [1/2]

template<MInt nDim>

Grid & maia::grid::Controller< nDim >::gridb ( )

inlineprivate

Definition at line 129 of file cartesiangridcontroller.h.

129{ return *m_grid; }

◆ gridb() [2/2]

template<MInt nDim>

const Grid & maia::grid::Controller< nDim >::gridb ( ) const

inlineprivate

Definition at line 130 of file cartesiangridcontroller.h.

130{ return *m_grid; }

◆ initDlbProperties()

template<MInt nDim>

void maia::grid::Controller< nDim >::initDlbProperties

private

Definition at line 532 of file cartesiangridcontroller.h.

◆ initTimers()

template<MInt nDim>

void maia::grid::Controller< nDim >::initTimers

private

Definition at line 977 of file cartesiangridcontroller.h.

                                  {
  NEW_TIMER_GROUP_NOCREATE(m_timerGroup, "GridController (noSolvers = " + std::to_string(noSolvers()) + ")");
  m_timers.fill(-1);
  NEW_TIMER_NOCREATE(m_timers[Timers::Controller], "total object lifetime", m_timerGroup);
  RECORD_TIMER_START(m_timers[Timers::Controller]);
 
  // Balancing timers
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::DLB], "DLB", m_timers[Timers::Controller]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Trigger], "Trigger", m_timers[Timers::DLB]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Partition], "Partition", m_timers[Timers::DLB]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Prepare], "Prepare", m_timers[Timers::DLB]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::BalanceGrid], "Balance grid", m_timers[Timers::DLB]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::BalanceSolvers], "Balance solvers", m_timers[Timers::DLB]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::BalanceCouplers], "Balance couplers", m_timers[Timers::DLB]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::IO], "I/O", m_timers[Timers::DLB]);
 
  // Adaptation timers
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Adaptation], "Adaptation", m_timers[Timers::Controller]);
  NEW_SUB_TIMER_NOCREATE(m_timers[Timers::MeshAdaptation], "Mesh Adaptation", m_timers[Timers::Adaptation]);
 
  const MInt noSolversAndCouplers = noSolvers() + noCouplers();
  m_solverTimerGroups.resize(noSolversAndCouplers);
  m_solverTimers.resize(noSolversAndCouplers);
 
  for(MInt b = 0; b < noSolversAndCouplers; b++) {
    m_solverTimerGroups[b] = -1;
    m_solverTimers[b].fill(-1);
 
    const MBool isSolver = (b < noSolvers());
    const MString groupName = (isSolver) ? "solverId = " + std::to_string(solver(b).solverId())
                                         : "couplerId = " + std::to_string(coupler(b - noSolvers()).couplerId());
 
    NEW_TIMER_GROUP_NOCREATE(m_solverTimerGroups[b], "GridController: " + groupName);
    NEW_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::Total], "total", m_solverTimerGroups[b]);
    const MInt solverTimer = m_solverTimers[b][SolverTimers::Total];
 
    // Balancing timers
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::DLB], "Dynamic load balancing", solverTimer);
    const MInt dlbTimer = m_solverTimers[b][SolverTimers::DLB];
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::Reset], "Reset solver", dlbTimer);
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::CalcDataSizes], "Determine data sizes", dlbTimer);
    const MInt dataSizeTimer = m_solverTimers[b][SolverTimers::CalcDataSizes];
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::CalcDataSizesMpiBlocking], "initial blocking MPI",
                           dataSizeTimer);
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::CalcDataSizesMpi], "remaining MPI", dataSizeTimer);
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::Balance], "Balance", dlbTimer);
    const MInt balanceTimer = m_solverTimers[b][SolverTimers::Balance];
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::Redistribute], "Redistribute data", balanceTimer);
    const MInt redistTimer = m_solverTimers[b][SolverTimers::Redistribute];
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::RedistributeMpiBlocking], "initial blocking MPI",
                           redistTimer);
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::RedistributeMpi], "remaining MPI", redistTimer);
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::BalancePre], "Balance pre", balanceTimer);
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::SetData], "Set data", balanceTimer);
    const MInt setDataTimer = m_solverTimers[b][SolverTimers::SetData];
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::SetDataMpiBlocking], "initial blocking MPI", setDataTimer);
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::SetDataMpi], "remaining MPI", setDataTimer);
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::BalancePost], "Balance post", balanceTimer);
 
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::DlbOther], "Other", dlbTimer);
 
    // Adaptation timers
    NEW_SUB_TIMER_NOCREATE(m_solverTimers[b][SolverTimers::Adaptation], "Adaptation", solverTimer);
  }
}

◆ isAdaptationTimeStep()

template<MInt nDim>

MBool maia::grid::Controller< nDim >::isAdaptationTimeStep

Definition at line 3406 of file cartesiangridcontroller.h.

                                             {
  TRACE();
 
  return (m_adaptationInterval > 0 && globalTimeStep % m_adaptationInterval == 0 && globalTimeStep > m_adaptationStart);
}

◆ isDlbTimeStep()

template<MInt nDim>

MBool maia::grid::Controller< nDim >::isDlbTimeStep

Definition at line 3381 of file cartesiangridcontroller.h.

                                      {
  TRACE();
  MBool dlbTimeStep = ((globalTimeStep - m_loadBalancingOffset) % m_loadBalancingInterval == 0)
                      && (m_loadBalancingStopTimeStep < 0 || globalTimeStep < m_loadBalancingStopTimeStep);
  if(m_balanceAfterAdaptation) {
    dlbTimeStep = dlbTimeStep
                  || (((globalTimeStep - m_lastAdaptationTimeStep) == m_loadBalancingOffset)
                      && (m_nAdaptationsSinceBalance >= m_balanceAdaptationInterval || m_dlbStep < 2)
                      && (m_loadBalancingStopTimeStep < 0 || globalTimeStep < m_loadBalancingStopTimeStep));
 
    if((m_syncTimerSteps && ((globalTimeStep - m_lastAdaptationTimeStep) < m_loadBalancingOffset)) || m_syncTimeSteps) {
      cerr0 << "Applying Barrier for correct timer!" << std::endl;
      if(m_syncTimeSteps && !((globalTimeStep - m_lastAdaptationTimeStep) < m_loadBalancingOffset)) {
        solver(0).startLoadTimer(AT_);
      }
      MPI_Barrier(gridb().mpiComm(), AT_);
      if(m_syncTimeSteps && !((globalTimeStep - m_lastAdaptationTimeStep) < m_loadBalancingOffset)) {
        solver(0).stopLoadTimer(AT_);
      }
    }
  }
 
  return dlbTimeStep;
}

◆ loadBalancingCalcNewGlobalOffsets()

template<MInt nDim>

void maia::grid::Controller< nDim >::loadBalancingCalcNewGlobalOffsets	(	const MLong *const	oldPartitionCellOffsets,
		const MLong *const	newPartitionCellOffsets,
		MLong *const	globalOffsets
	)

private

Author: Ansgar Niemoeller (ansgar) a.nie.nosp@m.moel.nosp@m.ler@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Parameters

[in]	oldPartitionCellOffsets	The current (old) partition cell offsets.
[in]	newPartitionCellOffsets	The new partition cell offsets.
[out]	globalOffsets	Pointer to storage of size noDomains()+1 that will hold the new global domain offsets on exit.

Definition at line 2687 of file cartesiangridcontroller.h.

                                                                                     {
  TRACE();
 
  const MInt oldNoPartitionCells = gridb().m_noPartitionCells;
  MLongScratchSpace partitionCellsGlobalId(oldNoPartitionCells, AT_, "partitionCellsGlobalId");
 
  for(MInt i = 0; i < oldNoPartitionCells; i++) {
    const MLong globalCellId = gridb().m_localPartitionCellGlobalIds[i]; // sorted by global id
 
    partitionCellsGlobalId(i) = globalCellId;
    ASSERT(i == 0 || partitionCellsGlobalId(i) > partitionCellsGlobalId(i - 1),
           "Partition cell global ids not in ascending order.");
  }
 
  // Fill with zeros
  std::fill_n(&globalOffsets[0], noDomains() + 1, 0);
 
  // Calculate global offsets for the given new partition cell offsets
  for(MInt i = 1; i < noDomains(); i++) {
    const MLong partitionCellId = newPartitionCellOffsets[i];
    // Determine if this partition cell is currently on this domain
    const MBool hasPartitionCell = (oldPartitionCellOffsets[domainId()] <= partitionCellId
                                    && partitionCellId < oldPartitionCellOffsets[domainId() + 1]);
 
    // Set global domain offset if the partition cell is present, else it remains zero
    if(hasPartitionCell) {
      // Determine local partition cell index and the global cell id
      const MLong partitionCellLocalId = partitionCellId - oldPartitionCellOffsets[domainId()];
      globalOffsets[i] = partitionCellsGlobalId[partitionCellLocalId];
 
      // Local cell id of the partition cell and its parent
      MInt currentId = gridb().m_localPartitionCellLocalIds[partitionCellLocalId];
      MInt parentId = gridb().a_parentId(currentId);
 
      MInt shift = 0;
      // Partition level shift: if the partition cell is not on the min-level check the global id of
      // its parent, if 'globalId(parent) == globalId(partitionCell) - 1' the parent is a partition
      // level ancestor that belongs to this domain, i.e., it has no offspring cells on the previous
      // domain. Continue to loop up the parent cells and check this condition to find the shift
      // that yields the correct domain offset.
      while(gridb().a_level(currentId) != gridb().minLevel()
            && gridb().a_globalId(parentId) == gridb().a_globalId(currentId) - 1) {
        shift++;
        currentId = parentId;
        parentId = gridb().a_parentId(currentId);
      }
 
      TERMM_IF_COND(shift > 0 && gridb().m_maxPartitionLevelShift == 0,
                    "Error: domain offset has a shift but the maximum partition level shift is 0.");
 
      // Correct global domain offset
      globalOffsets[i] -= shift;
    }
  }
 
  // Set only once since distributed via allreduce with sum to all domains
  if(domainId() == 0) {
    globalOffsets[noDomains()] = gridb().domainOffset(noDomains());
  }
 
  // Combine the new domain offset information of all ranks, i.e. each new
  // global domain offset is determined by the domain where the corresponding
  // partition cell is present, other domains will add a value of zero to this offset
  MPI_Allreduce(MPI_IN_PLACE, &globalOffsets[0], noDomains() + 1, type_traits<MLong>::mpiType(), MPI_SUM,
                gridb().mpiComm(), AT_, "MPI_IN_PLACE", "globalOffsets[0]");
}

◆ loadBalancingCalcNoCellsToSend()

template<MInt nDim>

void maia::grid::Controller< nDim >::loadBalancingCalcNoCellsToSend	(	const MLong *const	offsets,
		MInt *const	noCellsToSendByDomain,
		MInt *const	noCellsToReceiveByDomain,
		MInt *const	sortedCellId,
		MInt *const	bufferIdToCellId
	)

private

Author: Ansgar Niemoeller (ansgar) a.nie.nosp@m.moel.nosp@m.ler@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Parameters

[in]	offsets	New domain offsets.
[out]	noCellsToSendByDomains	Holds the number of cells to send to each domain.
[out]	noCellsToRecvByDomain	Holds the number of cells to receive from each domain.
[out]	sortedCellId	Buffer index for each cell, such that data is sorted by global id.
[out]	bufferIdToCellId	Reversed map of sortedCellId.

Definition at line 2768 of file cartesiangridcontroller.h.

                                                                                    {
  TRACE();
 
  const MInt noCells = gridb().treeb().size();
 
  MIntScratchSpace targetDomainsByCell(noCells, AT_, "targetDomainsByCell");
  std::fill_n(&targetDomainsByCell[0], noCells, -1);
  std::fill_n(&noCellsToSendByDomain[0], noDomains() + 1, 0);
 
  // Loop over all cells and determine to which domain they belong under the new domain
  // decomposition
  for(MInt cellId = 0; cellId < noCells; cellId++) {
    // Skip halo cells
    if(gridb().a_hasProperty(cellId, Cell::IsHalo)) {
      continue;
    }
 
    const MLong globalCellId = gridb().a_globalId(cellId);
    // Search for global cell id in domain offsets array
    auto lowerBound = std::lower_bound(&offsets[0], &offsets[0] + noDomains(), globalCellId);
    const MInt dist = std::distance(&offsets[0], lowerBound);
    // Check if this cell is a domain offset (i.e. in the offsets list)
    const MBool isDomainOffset = (*lowerBound == globalCellId);
    // Determine target domain id
    const MInt globalDomain = (isDomainOffset) ? dist : dist - 1;
 
    // Tag cell with target domain
    targetDomainsByCell[cellId] = globalDomain;
    // Increase number of cells to send to this domain
    noCellsToSendByDomain[globalDomain]++;
  }
 
  // Use all-to-all communication to determine how many cells to send/receive to/from each domain
  MPI_Alltoall(&noCellsToSendByDomain[0], 1, type_traits<MInt>::mpiType(), &noCellsToReceiveByDomain[0], 1,
               type_traits<MInt>::mpiType(), gridb().mpiComm(), AT_, "noCellsToSendByDomain[0]",
               "noCellsToReceiveByDomain[0]");
 
  // Determine total number of cells to send/recv (store as last entry)
  noCellsToSendByDomain[noDomains()] =
      std::accumulate(&noCellsToSendByDomain[0], &noCellsToSendByDomain[0] + noDomains(), 0);
  noCellsToReceiveByDomain[noDomains()] =
      std::accumulate(&noCellsToReceiveByDomain[0], &noCellsToReceiveByDomain[0] + noDomains(), 0);
 
  if(noCellsToSendByDomain[noDomains()] < 1) {
    TERMM(1, std::to_string(domainId()) + " noCellsToSend = " + std::to_string(noCellsToSendByDomain[noDomains()]));
  }
  if(noCellsToReceiveByDomain[noDomains()] < 1) {
    TERMM(1, std::to_string(domainId()) + " noCellsToRecv = " + std::to_string(noCellsToReceiveByDomain[noDomains()]));
  }
 
  // Sort cells depending on the domain to send to
  std::fill_n(&sortedCellId[0], noCells, -1);
 
  // @ansgar TODO labels:CONTROLLER,DLB this might be very inefficient for large number of cores
  MInt currentBufferId = 0;
  for(MInt dom = 0; dom < noDomains(); ++dom) {
    // Map from global cell id to local cell id for cells that will belong to the current domain
    std::map<MLong, MInt> cellMap;
    for(MInt cellId = 0; cellId < noCells; ++cellId) {
      // Halo cells are automatically skipped as they have tag[] of -1.
      // Halo cells will have a sortedCellId[] of -1.
      if(targetDomainsByCell[cellId] == dom) {
        cellMap[gridb().a_globalId(cellId)] = cellId;
      }
    }
 
    // Set buffer location for all cells in global id order -> received data is stored in
    // global id order and does not need to be resorted
    for(auto&& cell : cellMap) {
      const MInt cellId = cell.second;
      sortedCellId[cellId] = currentBufferId;
      bufferIdToCellId[currentBufferId] = cellId; // map from buffer id to cell id
      ++currentBufferId;
    }
  }
}

◆ loadBalancingPartition()

template<MInt nDim>

MBool maia::grid::Controller< nDim >::loadBalancingPartition	(	const MFloat *	loads,
		const MFloat	imbalance,
		MLong *const	partitionCellOffsets,
		MLong *const	globalIdOffsets
	)

private

Author: Ansgar Niemoeller (ansgar) a.nie.nosp@m.moel.nosp@m.ler@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Based on the domain loads and current domain decomposition determine a new partitioning that is expected to reduce imbalances among domains.

Parameters

[in]	loads	Vector of domain loads.
[in]	imbalance	Load imbalance among domains.
[out]	partitionCellOffsets	New partition-cell offsets (size noDomains()+1).
[out]	globalIdOffsets	New global domain offsets.

Available DLB partition methods:

DLB_PARTITION_WEIGHT
- compute new weights for the different load types
- determine the new workload of all min-cells and solve 1D partitioning problem
DLB_PARTITION_SHIFT_OFFSETS
- compute new weights and domain weights (quantifies processing capacities)
- determine summed load errors for each domain offset and shift offsets until the imbalance is predicted to be counterbalanced
DLB_PARTITION_DEFAULT
- combination of methods 1 and 2: 3 steps with pure weighting, the remaining steps use the offset shifting method
DLB_PARTITION_TEST
- for testing, use old partitioning to check if solver still works after reinitialization

Definition at line 345 of file cartesiangridcontroller.cpp.

                                                                                         {
  TRACE();
 
  using namespace maia;
 
  m_log << " * determine new partition" << std::endl;
 
  MLongScratchSpace localPartitionCellCounts(noDomains(), AT_, "localPartitionCellCounts");
  MLongScratchSpace localPartitionCellOffsets(noDomains() + 1, AT_, "localPartitionCellOffsets");
  // Determine the number of partition cells on all domains and the partition cell offsets
  gridb().determineNoPartitionCellsAndOffsets(&localPartitionCellCounts[0], &localPartitionCellOffsets[0]);
 
  const MInt noLocalPartitionCells = localPartitionCellCounts[domainId()];
 
  // Set the partition method
  MInt partitionMethod = m_dlbPartitionMethod;
 
  // Change the actual used partition method depending on the DLB step
  switch(m_dlbPartitionMethod) {
    case DLB_PARTITION_DEFAULT: {
      // Default partition method: combination of 3 steps with pure weighting,
      // the remaining steps use the offset shifting method
      const MInt noWeightPartitionSteps = 3;
      if(m_dlbStep < noWeightPartitionSteps) {
        partitionMethod = DLB_PARTITION_WEIGHT;
      } else {
        partitionMethod = DLB_PARTITION_SHIFT_OFFSETS;
      }
      break;
    }
    // Else: dont change the partition method
    case DLB_PARTITION_WEIGHT:
    case DLB_PARTITION_SHIFT_OFFSETS:
    case DLB_PARTITION_TEST:
      break;
    default:
      TERMM(1, "Unknown DLB partition method.");
      break;
  }
 
  // Initialize domain weights in first DLB step
  if(m_dlbStep == 0) {
    m_domainWeights.assign(noDomains(), 1.0);
    m_lastOffsetShiftDirection.assign(noDomains() + 1, 0);
  }
 
  // Increase DLB step
  m_dlbStep++;
 
  // Compute new partitioning with the activated partition method
  switch(partitionMethod) {
    case DLB_PARTITION_WEIGHT: {
      m_log << "Partition method #0 (compute weights, determine min cells workload and use the "
               "partition() algorithm)"
            << std::endl;
      // Use load quantities, i.e. number of DOF or number of active cells, and loads to estimate
      // new weights. The new domain distribution is then determined by the partition() algorithm.
      std::vector<MFloat> weights;
      computeWeights(loads, 1.0, weights);
 
      updateWeightsAndWorkloads(weights, false);
 
      // Note: in case of a partition level shift the global id offsets might not be correct here,
      // i.e., they are not corrected. This is done later in loadBalancingCalcNewGlobalOffsets().
      partition(&partitionCellOffsets[0], &globalIdOffsets[0], true);
      break;
    }
    case DLB_PARTITION_TEST: {
      // TESTING: use old offsets to check if solver still works after reinit
      for(MInt i = 0; i < noDomains() + 1; i++) {
        partitionCellOffsets[i] = localPartitionCellOffsets[i];
      }
      if(domainId() == 0) {
        std::cerr << "Using original partition!" << std::endl;
      }
      break;
    }
    case DLB_PARTITION_SHIFT_OFFSETS: {
      m_log << "Partition method #3 (compute weights and domain weights iteratively, determine "
               "summed load error at each domain offset and vary min cell offset such that "
               "predicted error is zero)"
            << std::endl;
 
      std::vector<MFloat> weights;
      std::vector<MFloat> oldWeights;
      std::vector<MFloat> oldDomainWeights;
      oldDomainWeights = m_domainWeights;
      std::fill(oldDomainWeights.begin(), oldDomainWeights.end(),
                1.0); // TODO labels:CONTROLLER,DLB reset in each step?
 
      MFloatScratchSpace domainWorkLoad(noDomains(), AT_, "domainWorkLoad");
      MFloat averageWorkLoad = -1.0;
 
      MFloatScratchSpace localPartitionCellsWorkload(noLocalPartitionCells, AT_, "localPartitionCellsWorkload");
      MFloatScratchSpace partitionCellsWorkload(gridb().m_noPartitionCellsGlobal, AT_, "partitionCellsWorkload");
 
      const MInt maxNoIt = 20;
      // Iterate until new weights and domain weights are found, i.e. converged
      for(MInt it = 0; it < maxNoIt; it++) {
        // Compute weights
        computeWeights(loads, m_domainWeights[domainId()], weights);
 
        // @ansgar_mb TODO labels:CONTROLLER,DLB test this
        updateWeightsAndWorkloads(weights, false);
 
        // Assemble local partition-cell workloads
        for(MInt i = 0; i < noLocalPartitionCells; i++) {
          const MLong globalPartitionCellId = gridb().m_localPartitionCellGlobalIds[i];
          const MInt partitionCellId = gridb().globalIdToLocalId(globalPartitionCellId, true);
          TERMM_IF_NOT_COND(gridb().a_hasProperty(partitionCellId, Cell::IsPartitionCell),
                            "Error: cell is not a partition cell.");
 
          const MFloat workload = gridb().a_workload(partitionCellId);
          TERMM_IF_NOT_COND(workload > 0.0, "Error: partition cell workload needs to be > 0.0");
          localPartitionCellsWorkload[i] = workload;
        }
 
        // @ansgar TODO labels:DLB Temporary fix; use MLong!? only an issue for more than 2billion partition cells...
        MIntScratchSpace localPartitionCellCounts_(noDomains(), AT_, "localPartitionCellCounts_");
        MIntScratchSpace localPartitionCellOffsets_(noDomains() + 1, AT_, "localPartitionCellOffsets_");
        for(MInt i = 0; i < noDomains(); i++) {
          localPartitionCellCounts_[i] = (MInt)localPartitionCellCounts[i];
          localPartitionCellOffsets_[i] = (MInt)localPartitionCellOffsets[i];
        }
        localPartitionCellOffsets_[noDomains()] = localPartitionCellOffsets[noDomains()];
 
        // Gather partition-cell workloads on root
        MPI_Gatherv(&localPartitionCellsWorkload[0], noLocalPartitionCells, type_traits<MFloat>::mpiType(),
                    &partitionCellsWorkload[0], &localPartitionCellCounts_[0], &localPartitionCellOffsets_[0],
                    type_traits<MFloat>::mpiType(), 0, gridb().mpiComm(), AT_, "localPartitionCellsWorkload[0]",
                    "partitionCellsWorkload[0]");
        MInt error = 0;
        if(domainId() == 0) {
          // Output some information on partition cell workloads
          const MFloat maxPartitionWorkload =
              *std::max_element(partitionCellsWorkload.begin(), partitionCellsWorkload.end());
          const MFloat minPartitionWorkload =
              *std::min_element(partitionCellsWorkload.begin(), partitionCellsWorkload.end());
          const MFloat avgPartitionWorkload =
              std::accumulate(partitionCellsWorkload.begin(), partitionCellsWorkload.end(), 0.0)
              / static_cast<MFloat>(gridb().m_noPartitionCellsGlobal);
          m_log << " * Maximum/minimum/average partition cell workload: " << maxPartitionWorkload << ", "
                << minPartitionWorkload << ", " << avgPartitionWorkload << std::endl;
          if(minPartitionWorkload <= 0) {
            m_log << globalTimeStep << " ERROR: minimum partition cell workload is " << minPartitionWorkload
                  << std::endl;
            cerr0 << globalTimeStep << " ERROR: minimum partition cell workload is " << minPartitionWorkload
                  << std::endl;
            error = 1;
          }
 
          // Compute domain workloads
          for(MInt i = 0; i < noDomains(); i++) {
            const MInt offset = localPartitionCellOffsets[i];
            const MInt count = localPartitionCellCounts[i];
            domainWorkLoad[i] =
                std::accumulate(&partitionCellsWorkload[offset], &partitionCellsWorkload[offset] + count, 0.0);
          }
          // Compute average workload
          averageWorkLoad = std::accumulate(domainWorkLoad.begin(), domainWorkLoad.end(), 0.0) / noDomains();
 
          m_log << " * Average domain workload: " << averageWorkLoad << std::endl;
 
          // Update domain weights
          for(MInt i = 0; i < noDomains(); i++) {
            m_domainWeights[i] = (m_useDomainFactor) ? loads[i] * averageWorkLoad / domainWorkLoad[i] : 1.0;
          }
 
          const MFloat domainWeightMinLoadThreshold = 0.2;
 
          MFloat averageDomainWeight = 0.0;
          MInt countDomainWeights = 0;
          for(MInt i = 0; i < noDomains(); i++) {
            if(loads[i] >= domainWeightMinLoadThreshold) {
              averageDomainWeight += m_domainWeights[i];
              countDomainWeights++;
            }
          }
          averageDomainWeight /= (MFloat)countDomainWeights;
 
          for(MInt i = 0; i < noDomains(); i++) {
            // For domains with low load reset domain weight since estimation might fail
            if(loads[i] < domainWeightMinLoadThreshold) {
              m_domainWeights[i] = 1.0;
            } else {
              // Scale domain weights with average such that the mean domain weight is one
              m_domainWeights[i] /= averageDomainWeight;
 
              // Weighting with previous values
              m_domainWeights[i] = 0.5 * (m_domainWeights[i] + oldDomainWeights[i]);
            }
          }
 
          const MFloat maxDomainWeight = *std::max_element(m_domainWeights.begin(), m_domainWeights.end());
          const MFloat minDomainWeight = *std::min_element(m_domainWeights.begin(), m_domainWeights.end());
          m_log << " * Domain weights: max = " << maxDomainWeight << ", min = " << minDomainWeight << std::endl;
        }
 
        // Broadcast error flag
        MPI_Bcast(&error, 1, type_traits<MInt>::mpiType(), 0, gridb().mpiComm(), AT_, "error");
        if(error > 0) {
          return false;
        }
 
        // Broadcast new domain weights
        MPI_Bcast(&m_domainWeights[0], noDomains(), type_traits<MFloat>::mpiType(), 0, gridb().mpiComm(), AT_,
                  "m_domainWeights[0]");
 
        // Compute residuals and check for convergence
        if(it > 0) {
          const MInt noWeights = weights.size();
          MFloat weightResidual = 0.0;
          for(MInt i = 0; i < noWeights; i++) {
            weightResidual += ABS(weights[i] - oldWeights[i]);
          }
          weightResidual /= noWeights;
 
          MFloat domainWeightResidual = 0.0;
          for(MInt i = 0; i < noDomains(); i++) {
            domainWeightResidual += ABS(m_domainWeights[i] - oldDomainWeights[i]);
          }
          domainWeightResidual /= noDomains();
 
          const MFloat weightResThreshold = 1e-3;
          const MFloat domainWeightResThreshold = 1e-2;
          if(weightResidual < weightResThreshold && domainWeightResidual < domainWeightResThreshold) {
            m_log << " * weights and domain weights converged, residuals: " << weightResidual << ", "
                  << domainWeightResidual << std::endl;
            break;
          } else if(it == maxNoIt - 1) {
            m_log << " * WARNING: Weights and domain weights not converged, but "
                     "maximum number of iterations reached. "
                  << weightResidual << " " << domainWeightResidual << std::endl;
          } else {
            m_log << " * Weights and domain weights iteration " << it << " " << weightResidual << " "
                  << domainWeightResidual << std::endl;
          }
        }
        oldWeights = weights;
        oldDomainWeights = m_domainWeights;
      }
 
      if(m_debugBalance) {
        const MInt noLoadTypes = globalNoLoadTypes();
        MIntScratchSpace loadQuantities(noLoadTypes, AT_, "loadQuantities");
        getLoadQuantities(&loadQuantities[0]);
        storeLoadsAndWeights(&loads[0], noLoadTypes, &loadQuantities[0], m_domainWeights[domainId()], &weights[0]);
      }
 
      // New Weights and domain weights are computed, rank 0 has already the partition
      // cell workload information, now determine the new partition cell offsets
      if(domainId() == 0) {
        MFloatScratchSpace errors(noDomains(), AT_, "errors");
        MFloatScratchSpace summedErrors(noDomains() + 1, AT_, "summedErrors");
        MFloatScratchSpace summedErrAbs(noDomains() + 1, AT_, "summedErrAbs");
        MIntScratchSpace summedErrDir(noDomains() + 1, AT_, "summedErrDir");
 
        // Compute load error for each domain ...
        for(MInt i = 0; i < noDomains(); i++) {
          errors[i] = loads[i] - 1.0;
        }
        // ... and summed load errors (and their absolut value) at each domain offset
        summedErrors[0] = 0.0;
        summedErrAbs[0] = 0.0;
        summedErrDir[0] = 0;
        for(MInt i = 1; i < noDomains() + 1; i++) {
          summedErrors[i] = summedErrors[i - 1] + errors[i - 1];
          summedErrAbs[i] = ABS(summedErrors[i]);
          summedErrDir[i] = -signum(summedErrors[i]);
        }
 
        const MFloat maxSummedErr = *std::max_element(summedErrors.begin(), summedErrors.end());
        const MFloat minSummedErr = *std::min_element(summedErrors.begin(), summedErrors.end());
        m_log << globalTimeStep << " * Summed load deviations: max = " << maxSummedErr << ", min = " << minSummedErr
              << std::endl;
 
        // Threshold below which an offset will not be changed
        const MFloat errorThreshold = m_dlbImbalanceThreshold;
 
        // Restrict the amount of imbalance to be shifted in a single step for higher degrees of
        // parallelism
        const MFloat restrictionFactor = 1.0 - std::min(0.5, 0.005 * noDomains());
 
        MFloat loadVariance = 0.0;
        for(MInt i = 0; i < noDomains(); i++) {
          loadVariance += POW2(loads[i] - 1.0);
        }
        loadVariance /= noDomains();
        const MFloat loadStd = std::sqrt(loadVariance);
        m_log << globalTimeStep << " * load distribution: variance = " << loadVariance << ", stdev = " << loadStd
              << std::endl;
 
        const MFloat imbalancePenaltyThreshold = 2.0 * errorThreshold;
        const MFloat imbPenaltyFactor =
            (imbalance < imbalancePenaltyThreshold)
                ? (1.0 + (imbalancePenaltyThreshold - imbalance) / imbalancePenaltyThreshold)
                : 1.0;
        m_log << globalTimeStep << " * Imbalance penalty factor " << imbPenaltyFactor << " (imbalance " << imbalance
              << "); restriction factor " << restrictionFactor << std::endl;
 
        // Final local shift for refining the best found partitioning
        const MBool finalLocalShift =
            m_dlbNoFinalLocalShifts > 0
            && globalTimeStep >= m_loadBalancingStopTimeStep - m_loadBalancingInterval * (m_dlbNoFinalLocalShifts + 1);
        // Intermediate local shift
        const MBool intermediateLocalShift =
            !finalLocalShift && (m_dlbNoLocalShifts > 0 && (m_dlbStep % (1 + m_dlbNoLocalShifts) != 1));
        // Check for global shift step
        const MBool globalShiftStep = !finalLocalShift && !intermediateLocalShift;
 
        m_log << globalTimeStep << " * Shift step: global " << globalShiftStep << "; intermediate local "
              << intermediateLocalShift << "; final local " << finalLocalShift << std::endl;
 
        // Flag to indicate if a global shift is required for each offset
        ScratchSpace<MInt> globalShiftFlag(noDomains() + 1, AT_, "globalShiftFlag");
        globalShiftFlag.fill(0);
        // Minimum distance of an offset to the next offset with a global shift
        ScratchSpace<MInt> minGlobalDist(noDomains() + 1, AT_, "minGlobalDist");
        minGlobalDist.fill(noDomains());
 
        { // Determine offsets with a global shift
          ScratchSpace<MInt> globalShiftFlagTmp(noDomains() + 1, AT_, "globalShiftFlagTmp");
          globalShiftFlagTmp.fill(0);
 
          MInt noGlobalShifts = 0;
          // 1. Determine offsets which require a global shift, i.e., the load is too imbalanced
          // among the neighboring domains
          for(MInt i = 1; i < noDomains(); i++) {
            const MInt dist = 4; // consider offsets in this distance to both sides
            const MInt firstOffset = std::max(0, i - dist);
            const MInt lastOffset = std::min(noDomains(), i + dist);
            // Difference in summed imbalance
            const MFloat summedErrorDiff = ABS(summedErrors[lastOffset] - summedErrors[firstOffset]);
            const MBool sameSign = (summedErrDir[lastOffset] == summedErrDir[firstOffset]);
            const MBool absErrCondition = (summedErrAbs[i] > errorThreshold);
 
            const MFloat smoothShiftThreshold = (lastOffset - firstOffset) * errorThreshold; // 0.5
            const MBool smoothShiftCondition =
                (absErrCondition && (sameSign && summedErrorDiff > smoothShiftThreshold));
            // Global shift required if the difference in the summed imbalance is too high
            MBool globalShiftCondition = false;
            if(globalShiftStep) {
              globalShiftCondition = (m_dlbSmoothGlobalShifts) ? smoothShiftCondition : absErrCondition;
            }
 
            globalShiftFlagTmp[i] = (globalShiftCondition) ? summedErrDir[i] : 0;
            if(globalShiftFlagTmp[i] != 0) {
              noGlobalShifts++;
              if(m_debugBalance) {
                m_log << globalTimeStep << " * global shift " << i << " diff=" << summedErrorDiff
                      << " tr=" << smoothShiftThreshold << " s_first=" << summedErrors[firstOffset]
                      << " s_last=" << summedErrors[lastOffset] << " s_abs=" << summedErrAbs[i] << std::endl;
              }
            }
          }
          m_log << globalTimeStep << " * number of global shifts " << noGlobalShifts << std::endl;
 
          // 2. For each offset marked with a global shift, mark also the nearby offsets up to a
          // certain distance while reducing the size of the shift such that global shifts are
          // smoothed out over some neighborhood. For this the minimum distance to a global shift is
          // stored, which later reduces the summed imbalance that needs to be counterbalanced by
          // the shift.
          for(MInt i = 1; i < noDomains(); i++) {
            // Check for offset with global shift
            if(globalShiftFlagTmp[i] != 0) {
              minGlobalDist[i] = 0;
              globalShiftFlag[i] = globalShiftFlagTmp[i];
 
              // Distance over which to smooth out the global shift
              const MInt dist = std::ceil(summedErrAbs[i] / errorThreshold);
              if(m_debugBalance) {
                m_log << globalTimeStep << " * global shift " << i << " marking distance " << dist << " "
                      << summedErrAbs[i] << std::endl;
              }
 
              // Loop over lower offsets
              for(MInt nId = i - 1; nId >= std::max(i - dist, 1); nId--) {
                const MInt distance = ABS(nId - i);
 
                // Stop marking offsets if shift direction is not the same anymore or if the
                // summed imbalance is too small for this distance, which would result in a summed
                // imbalance with opposite sign
                if(summedErrDir[nId] != summedErrDir[i] || errorThreshold * distance > summedErrAbs[nId]) {
                  break;
                }
 
                minGlobalDist[nId] = std::min(minGlobalDist[nId], distance);
                globalShiftFlag[nId] = summedErrDir[nId];
              }
              // Loop over higher offsets
              for(MInt nId = i + 1; nId <= std::min(i + dist, noDomains() - 1); nId++) {
                const MInt distance = ABS(nId - i);
 
                if(summedErrDir[nId] != summedErrDir[i] || errorThreshold * distance > summedErrAbs[nId]) {
                  break;
                }
 
                minGlobalDist[nId] = std::min(minGlobalDist[nId], distance);
                globalShiftFlag[nId] = summedErrDir[nId];
              }
            }
          }
 
          MInt noGlobalShiftsTotal = 0;
          for(MInt i = 1; i < noDomains(); i++) {
            if(globalShiftFlag[i] != 0) {
              noGlobalShiftsTotal++;
            }
          }
          m_log << globalTimeStep << " * number of total global shifts " << noGlobalShiftsTotal << std::endl;
        }
 
        // Variables for some user output about the shifts
        MInt noShifts = 0;
        MInt noShiftsWithoutChange = 0;
        MFloat maxWeight = 0.0;
        MFloat maxDiff = 0.0;
        // Storage for shift information that is written to file in debug mode
        ScratchSpace<MInt> shifts(noDomains() + 1, AT_, "shifts");
        shifts.fill(0);
        ScratchSpace<MFloat> shiftedWorkload(noDomains() + 1, AT_, "shiftedWorkload");
        shiftedWorkload.fill(0.0);
 
        partitionCellOffsets[0] = 0;
        // vary min cell offsets individually
        for(MInt offsetId = 1; offsetId < noDomains(); offsetId++) {
          MFloat summedWorkload = 0.0;
          // Determine current workload of the domain left to this offset (with new previous offset)
          for(MInt i = partitionCellOffsets[offsetId - 1]; i < localPartitionCellOffsets[offsetId]; i++) {
            summedWorkload += (m_useDomainFactor) ? m_domainWeights[offsetId - 1] * partitionCellsWorkload[i]
                                                  : partitionCellsWorkload[i];
          }
 
          // Check for a global or local shift
          const MBool globalShiftCondition = (globalShiftFlag[offsetId] != 0);
          MBool localShiftCondition = false;
          MInt localShiftDir = 0;
          MFloat localShiftDiff = 0.0;
 
          // No global shift for this offset and the neighboring ones, use a local shift to
          // distribute load among neighboring domains
          if(!globalShiftCondition && globalShiftFlag[offsetId - 1] == 0 && globalShiftFlag[offsetId + 1] == 0) {
            const MFloat errLeft = errors[offsetId - 1];
            const MFloat errRight = errors[offsetId];
            localShiftDiff = 0.5 * (errLeft - errRight);
 
            if(errLeft > errorThreshold && errRight < errLeft) {
              // overload left
              localShiftDir = -1;
            } else if(errRight > errorThreshold && errRight > errLeft) {
              // overload right
              localShiftDir = 1;
            } else if(errLeft < -errorThreshold && errLeft < errRight) {
              // underload left (l/r no overload)
              localShiftDir = 1;
            } else if(errRight < -errorThreshold && errRight < errLeft) {
              // underload right (l/r no overload)
              localShiftDir = -1;
            }
            localShiftCondition = (localShiftDir != 0);
          }
 
          const MBool shift = globalShiftCondition || localShiftCondition;
 
          if(m_debugBalance) {
            m_log << globalTimeStep << " * shift conditions " << offsetId << " shift=" << shift
                  << " g=" << globalShiftCondition << " l=" << localShiftCondition << " err=" << summedErrors[offsetId]
                  << " minDist=" << minGlobalDist[offsetId] << " lDir=" << localShiftDir << " diff=" << localShiftDiff
                  << " eL=" << errors[offsetId - 1] << " eR=" << errors[offsetId] << std::endl;
          }
 
          if(shift) {
            // Set shift direction depending on global or local shift
            const MInt dir = (globalShiftCondition) ? summedErrDir[offsetId] : localShiftDir;
            MFloat penaltyFactor = 1.25;
 
            const MInt lastDir = (globalShiftCondition) ? m_lastOffsetShiftDirection[offsetId] : 0;
            // Check for current shift in opposite direction of last shift (only if shift is global)
            if(lastDir != dir && lastDir != 0) {
              penaltyFactor = 2.0;
              if(m_debugBalance) {
                m_log << "Opposite shift " << offsetId << ", lastDir " << lastDir << ", dir " << dir << std::endl;
              }
            }
            penaltyFactor *= imbPenaltyFactor;
 
            // Determine 'new' and 'old' domain id depending on the direction in
            // which the partition cell offset will be shifted
            const MInt domainOld = (dir < 0) ? offsetId - 1 : offsetId;
            const MInt domainNew = (dir < 0) ? offsetId : offsetId - 1;
 
            // Initialize estimated summed load error
            // To smooth out the globally enforced shifts, the error threshold times
            // the minimum distance to the next domain which requires a global shift is subtracted,
            // i.e., the shift of an offset further away from a global shift is reduced
            const MFloat summedErrorTmp =
                restrictionFactor * (summedErrors[offsetId] + dir * errorThreshold * minGlobalDist[offsetId]);
            //(dir > 0) ? std::max(-2.0, summedErrorTmp) : std::min(2.0, summedErrorTmp);
            const MFloat summedError = summedErrorTmp;
            MFloat oldEstimate = (globalShiftCondition) ? summedError : localShiftDiff;
            const MFloat initialEstimate = oldEstimate;
 
            MFloat penalty = 1.0;
 
            // Domain id of the currently considered partition cell
            MInt partitionCellDomain = domainOld;
 
            if(m_debugBalance) {
              m_log << "DLB_DEBUG: Start offset shift " << offsetId << ", domainOld " << domainOld << ", domainNew "
                    << domainNew << ", oldEstimate " << oldEstimate << ", summedError " << summedErrors[offsetId]
                    << ", penaltyFactor " << penaltyFactor << std::endl;
            }
 
            // Shift the partition cell offset such that the new estimated summed load error is
            // minimized
            for(MInt i = 1;; i++) {
              // TODO labels:DLB check if first/last partition cell is reached and stop!
              // Currently considered partition cell id
              const MInt newPartitionCellId = localPartitionCellOffsets[offsetId] + dir * i;
 
              // Update domain id on which the current partition cell belongs
              if(dir > 0 && newPartitionCellId == localPartitionCellOffsets[partitionCellDomain + 1]) {
                // Next partition cell offsets reached, increase domain id
                partitionCellDomain++;
              } else if(dir < 0 && newPartitionCellId == localPartitionCellOffsets[partitionCellDomain] - 1) {
                // Previous partition cell offset exceeded, decrease domain id
                partitionCellDomain--;
              }
 
              const MBool domainFactorShift = true; //(oldEstimate < 1.0); // TODO labels:CONTROLLER,DLB
              // Performance factor of domains
              const MFloat domainFactor = (m_useDomainFactor && domainFactorShift)
                                              ? (m_domainWeights[domainNew] / m_domainWeights[partitionCellDomain])
                                              : 1.0;
 
              // Compute new estimate of summed load error
              const MFloat estimate = oldEstimate
                                      + dir * penaltyFactor * penalty * domainFactor * loads[partitionCellDomain]
                                            * partitionCellsWorkload[newPartitionCellId]
                                            / domainWorkLoad[partitionCellDomain];
              const MFloat diff = estimate - oldEstimate;
 
              // Reset summed workload if the previous new offset is still larger than the current
              // one
              if(newPartitionCellId <= partitionCellOffsets[offsetId - 1]) {
                summedWorkload = 0;
              }
              const MFloat workloadIncrement =
                  (m_useDomainFactor) ? m_domainWeights[offsetId - 1] * partitionCellsWorkload[newPartitionCellId]
                                      : partitionCellsWorkload[newPartitionCellId];
              const MFloat summedWorkloadNew = summedWorkload + dir * workloadIncrement;
              const MFloat summedWorkloadRel = summedWorkload / averageWorkLoad;
 
#ifndef NDEBUG
              // Only useful for debugging small cases
              if(m_debugBalance) {
                m_log << "DLB_DEBUG: offset " << offsetId << ", partCellId " << newPartitionCellId
                      << ", partCellDomain " << partitionCellDomain << ", domainFactor " << domainFactor
                      << ", oldEstimate " << oldEstimate << ", estimate " << estimate << ", diff " << diff
                      << ", penalty " << penalty << ", partCellWorkload " << partitionCellsWorkload[newPartitionCellId]
                      << ", domainWorkload " << domainWorkLoad[partitionCellDomain] << ", load "
                      << loads[partitionCellDomain] << std::endl;
              }
#endif
 
              // Accept the previous new partition cell offset if the predicted summed load error
              // reached its minimum or if the previous domain has only a single parititon cell left
              // (when decreasing the current partition cell offset).
              const MBool prevOffsetReached = (dir < 0 && newPartitionCellId == partitionCellOffsets[offsetId - 1]);
              const MBool lastCellReached = (newPartitionCellId == gridb().m_noPartitionCellsGlobal - 1);
 
              // Continue shift if the workload is too large
              const MBool workloadCheck = (m_dlbMaxWorkloadLimit > 1.0 && globalShiftCondition)
                                              ? (summedWorkloadRel < m_dlbMaxWorkloadLimit)
                                              : true;
 
              if(prevOffsetReached || lastCellReached || (ABS(estimate) >= ABS(oldEstimate) && workloadCheck)
                 || (dir == 1 && !workloadCheck)) {
                // New domain offset
                partitionCellOffsets[offsetId] = localPartitionCellOffsets[offsetId] + dir * (i - 1);
 
                if(m_debugBalance) {
                  m_log << globalTimeStep << " * shifted offset " << offsetId << " " << dir * (i - 1) << " "
                        << localPartitionCellOffsets[offsetId] << " " << partitionCellOffsets[offsetId] << " "
                        << summedErrors[offsetId] << " " << estimate << " " << oldEstimate << " " << penalty << " "
                        << summedWorkloadRel << std::endl;
                }
                shifts[offsetId] = dir * (i - 1);
 
                // Set direction of shift (if it was a global shift)
                m_lastOffsetShiftDirection[offsetId] = (i > 1 && globalShiftCondition) ? dir : 0;
 
                // Make sure the current domain has at least a single min cell
                // TODO labels:DLB handle last domain offset!
                if(partitionCellOffsets[offsetId] <= partitionCellOffsets[offsetId - 1]) {
                  partitionCellOffsets[offsetId] = partitionCellOffsets[offsetId - 1] + 1;
                  shifts[offsetId] = partitionCellOffsets[offsetId] - localPartitionCellOffsets[offsetId];
                  m_log << "WARNING: setting domain offset #" << offsetId
                        << " with just a single partition cell on this domain. " << partitionCellOffsets[offsetId - 1]
                        << " " << partitionCellOffsets[offsetId] << std::endl;
                } else if(partitionCellOffsets[offsetId] == partitionCellOffsets[offsetId - 1] + 1) {
                  m_log << "WARNING: domain #" << offsetId - 1 << " has just a single partition cell. "
                        << partitionCellOffsets[offsetId - 1] << " " << partitionCellOffsets[offsetId] << std::endl;
                }
 
                if(shifts[offsetId] != 0) {
                  noShifts++;
                } else {
                  noShiftsWithoutChange++;
                }
 
                // Continue with next domain offset
                break;
              } else {
                // Update estimate and continue with next partition cell
                oldEstimate = estimate;
                // Update summed workload for the domain left to the offset
                summedWorkload = summedWorkloadNew;
 
                // Determine max difference and weight for some user output
                maxDiff = std::max(maxDiff, ABS(diff));
                maxWeight = std::max(maxWeight, partitionCellsWorkload[newPartitionCellId]);
                // Increase shifted workload
                shiftedWorkload[offsetId] += dir * partitionCellsWorkload[newPartitionCellId];
 
                // Compute new penalty factor: 1.0/cos(pi/2*((s-e)/s))
                // Goes to infinity as the estimate goes to zero and is intended to prevent
                // 'overshooting' of the domain offsets, i.e. shifting the domain offsets too far
                const MFloat s = initialEstimate;
                const MFloat x = ABS((ABS(s) - ABS(estimate)) / ABS(s));
                penalty = 1.0 / cos(M_PI / 2.0 * x);
              }
            }
          } else {
            // Accept old partition cell offset
            partitionCellOffsets[offsetId] = localPartitionCellOffsets[offsetId];
            m_lastOffsetShiftDirection[offsetId] = 0;
            shifts[offsetId] = 0;
 
            if(m_debugBalance) {
              m_log << globalTimeStep << " * keep offset " << offsetId << " " << summedErrors[offsetId] << std::endl;
            }
          }
        }
 
        // Correct offsets which were not shifted to prevent domains with an invalid number of cells!
        for(MInt offsetId = 1; offsetId < noDomains(); offsetId++) {
          if(partitionCellOffsets[offsetId] <= partitionCellOffsets[offsetId - 1]) {
            partitionCellOffsets[offsetId] = partitionCellOffsets[offsetId - 1] + 1;
            shifts[offsetId] = partitionCellOffsets[offsetId] - localPartitionCellOffsets[offsetId];
            m_log << "WARNING: correcting domain offset #" << offsetId
                  << " with now just a single partition cell on this domain. " << partitionCellOffsets[offsetId - 1]
                  << " " << partitionCellOffsets[offsetId] << std::endl;
          }
        }
 
        m_log << " * number of shifted offsets: " << noShifts
              << "; offsets to shift without change: " << noShiftsWithoutChange << std::endl;
        m_log << " * max weight during shifting: " << maxWeight << "; max difference: " << maxDiff << std::endl;
 
        // Debug: store partition offset shifts to file
        if(m_debugBalance) {
          FILE* shiftsFile = nullptr;
          const MString fileName = "shifts_" + std::to_string(globalTimeStep) + ".dat";
          shiftsFile = fopen(fileName.c_str(), "w");
 
          fprintf(shiftsFile, "# 1:offsetId 2:partitionCellShift 3:shiftedWorkload "
                              "4:summedImbalance 5:load 6:domainWorkload 7:domainWeight\n");
          for(MInt i = 0; i < noDomains(); i++) {
            fprintf(shiftsFile, "%d %d %.8e %.8e %.8e %.8e %.8e\n", i, shifts[i], shiftedWorkload[i], summedErrors[i],
                    loads[i], domainWorkLoad[i], m_domainWeights[i]);
          }
          fclose(shiftsFile);
        }
      }
 
      // Broadcast new partition cell offsets
      MPI_Bcast(&partitionCellOffsets[0], noDomains(), type_traits<MLong>::mpiType(), 0, gridb().mpiComm(), AT_,
                "partitionCellOffsets[0]");
      partitionCellOffsets[noDomains()] = gridb().m_noPartitionCellsGlobal;
 
      break;
    }
    default: {
      TERMM(1, "Partition method not known.");
      break;
    }
  }
 
  // Check if new distribution is different to the old one
  const MBool newPartition =
      !std::equal(&localPartitionCellOffsets[0], &localPartitionCellOffsets[0] + noDomains(), &partitionCellOffsets[0]);
 
  m_log << " * checking new partition" << std::endl;
  MBool validPartition =
      (partitionCellOffsets[0] == 0) && (partitionCellOffsets[noDomains() - 1] < gridb().m_noPartitionCellsGlobal);
  if(!validPartition) {
    m_log << "   * invalid first/last domain offset:" << partitionCellOffsets[0] << " "
          << partitionCellOffsets[noDomains() - 1] << "; noPartitionCellsGlobal = " << gridb().m_noPartitionCellsGlobal
          << std::endl;
  }
  // Check if new partitioning is valid, i.e. the domain offsets are in ascending order without duplicates
  for(MInt i = 1; i < noDomains(); i++) {
    if(partitionCellOffsets[i] <= partitionCellOffsets[i - 1]) {
      validPartition = false;
      m_log << "   * invalid offsets " << i - 1 << " " << partitionCellOffsets[i - 1] << " and " << i << " "
            << partitionCellOffsets[i] << std::endl;
    }
  }
 
  // Return if partitioning is not valid
  // Note: this may not work if m_forceLoadBalancing is active
  if(!validPartition) {
    m_log << " * new partition is not valid!" << std::endl;
    return false;
  }
 
  // Return if partitioning did not change and DLB is not forced
  if(!newPartition && !m_forceLoadBalancing) {
    m_log << " * partitioning did not change, return." << std::endl;
    return false;
  } else {
    m_log << " * new partition (" << newPartition << "), forced balance (" << m_forceLoadBalancing << ")" << std::endl;
  }
 
  // Calculate new global domain offsets from the new partition cell offsets
  loadBalancingCalcNewGlobalOffsets(&localPartitionCellOffsets[0], &partitionCellOffsets[0], &globalIdOffsets[0]);
 
  return true;
}

◆ logTimerStatistics()

template<MInt nDim>

void maia::grid::Controller< nDim >::logTimerStatistics ( const MString & status )

Definition at line 3330 of file cartesiangridcontroller.h.

                                                               {
  m_log << "DLB: Timer statistics ";
  if(!status.empty()) m_log << status << " ";
  m_log << "at global time step " << globalTimeStep << std::endl;
 
  // Accumulate timer records of all dlb timers
  MFloat localRunTime = 0.0;
  MFloat localIdleTime = 0.0;
  const MInt noDlbTimers = maia::dlb::g_dlbTimerController.noDlbTimers();
  for(MInt i = 0; i < noDlbTimers; i++) {
    const MFloat loadRecord = maia::dlb::g_dlbTimerController.returnLoadRecord(i, m_loadBalancingTimerMode);
    const MFloat idleRecord = maia::dlb::g_dlbTimerController.returnIdleRecord(i, m_loadBalancingTimerMode);
 
    if(m_loadBalancingMode == 1 && !m_testDynamicLoadBalancing && (loadRecord < 0.0 || idleRecord < 0.0)) {
      TERMM(1, "Load/Idle record for dlb timer #" + std::to_string(i) + " is less than zero on global domain #"
                   + std::to_string(domainId()));
    }
 
    localRunTime += loadRecord;
    localIdleTime += idleRecord;
  }
 
  const MInt noTimeSteps = globalTimeStep - m_dlbLastResetTimeStep;
 
  MFloat timePerStep = (localRunTime + localIdleTime) / noTimeSteps;
  MFloat maxRunTime = localRunTime / noTimeSteps;
  const MFloat localTimePerStep = timePerStep;
 
  // Communicate and take maximum value -> assure same value on all domains
  MPI_Allreduce(MPI_IN_PLACE, &timePerStep, 1, MPI_DOUBLE, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "timePerStep");
  MPI_Allreduce(MPI_IN_PLACE, &maxRunTime, 1, MPI_DOUBLE, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "maxRunTime");
 
  MFloat maxIdleTime = localIdleTime / noTimeSteps;
  MFloat minIdleTime = localIdleTime / noTimeSteps;
  MPI_Allreduce(MPI_IN_PLACE, &maxIdleTime, 1, MPI_DOUBLE, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "maxIdleTime");
  MPI_Allreduce(MPI_IN_PLACE, &minIdleTime, 1, MPI_DOUBLE, MPI_MIN, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "minIdleTime");
 
 
  m_log << globalTimeStep << " Average time per step: " << timePerStep << " ; local: " << localTimePerStep
        << ", idle/comp = " << localIdleTime / localRunTime
        << ", idle/timePerStep = " << localIdleTime / (noTimeSteps * timePerStep) << std::endl;
  m_log << globalTimeStep << " Relative idle time: max = " << maxIdleTime / timePerStep << ", min "
        << minIdleTime / timePerStep << std::endl;
  m_log << globalTimeStep << " maxRunTime " << maxRunTime << std::endl;
}

◆ needLoadBalancing()

template<MInt nDim>

MBool maia::grid::Controller< nDim >::needLoadBalancing	(	const MFloat	localRunTime,
		MFloat *const	loads,
		MFloat &	imbalance
	)

private

Author: Ansgar Niemoeller (ansgar) a.nie.nosp@m.moel.nosp@m.ler@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Parameters

[in]	localRunTime	The local run time, i.e. computation time.
[out]	load	Computed domain loads.
[out]	imbalance	Computed imbalance.

Definition at line 2519 of file cartesiangridcontroller.h.

                                                                                                           {
  TRACE();
 
  // Nothing to do in serial
  if(noDomains() == 1) {
    return false;
  }
 
  MFloatScratchSpace localRunTimes(noDomains(), FUN_, "localRunTimes");
  MFloatScratchSpace localIdleTimes(noDomains(), FUN_, "localIdleTimes");
 
  // Gather the local runtimes on all ranks
  MPI_Allgather(&localRunTime, 1, type_traits<MFloat>::mpiType(), &localRunTimes[0], 1, type_traits<MFloat>::mpiType(),
                gridb().mpiComm(), AT_, "localRunTime", "localRunTimes[0]");
 
  // Calculate average and maximum runtime of all ranks
  const MFloat averageRunTime = std::accumulate(localRunTimes.begin(), localRunTimes.end(), 0.0) / noDomains();
  const MFloat maxRunTime = *std::max_element(localRunTimes.begin(), localRunTimes.end());
 
  // Imbalance percentage: I_% = (t_max - t_avg)/t_max * N/(N-1)
  // (DeRose 2007, Detecting Application Load Imbalance on High End Massively Parallel Systems)
  imbalance = (maxRunTime - averageRunTime) / maxRunTime * noDomains() / (noDomains() - 1.0);
  const MBool loadBalance = (imbalance > m_dlbImbalanceThreshold);
 
  // Compute domain loads as fraction of local run time to average run time
  for(MInt i = 0; i < noDomains(); i++) {
    const MFloat load = localRunTimes[i] / averageRunTime;
    loads[i] = load;
  }
 
  const MFloat maxLoad = *std::max_element(&loads[0], &loads[0] + noDomains());
  const MFloat minLoad = *std::min_element(&loads[0], &loads[0] + noDomains());
 
  char imb[10];
  std::sprintf(imb, "%.2f", imbalance * 100.0);
  std::stringstream message;
  message << globalTimeStep << " * Imbalance percentage: " << imb << "%, t_avg = " << averageRunTime
          << ", t_max = " << maxRunTime << std::endl;
  message << globalTimeStep << " * Loads: max = " << maxLoad << ", min = " << minLoad << std::endl;
  m_log << message.str();
  cerr0 << message.str();
 
  // TODO labels:CONTROLLER,DLB generalize histogram output for other purposes?
  const MFloat binWidth = 0.05;
  const MInt noBins = std::ceil(maxLoad / binWidth) + 1;
  std::vector<MInt> loadBins(noBins);
  std::fill(loadBins.begin(), loadBins.end(), 0);
  for(MInt i = 0; i < noDomains(); i++) {
    const MInt binId = std::floor(loads[i] / binWidth);
    loadBins[binId] += 1;
  }
 
  const MInt sumBins = std::accumulate(loadBins.begin(), loadBins.end(), 0);
  if(sumBins != noDomains()) {
    m_log << "ERROR in load bins, count " << sumBins << " != " << noDomains() << std::endl;
  }
 
  const MInt noBinsPrev = m_previousLoadBins.size();
  const MInt maxNoBins = std::max(noBins, noBinsPrev);
  const MInt maxBinCountCurr = *std::max_element(loadBins.begin(), loadBins.end());
  const MInt maxBinCountPrev =
      (noBinsPrev > 0) ? *std::max_element(m_previousLoadBins.begin(), m_previousLoadBins.end()) : 0;
  const MInt maxBinCount = std::max(maxBinCountCurr, maxBinCountPrev);
 
  const MInt firstBin =
      (std::find_if(loadBins.begin(), loadBins.end(), [](MInt i) { return i > 0; }) - loadBins.begin()) - 1;
  const MInt firstBinPrev =
      (std::find_if(m_previousLoadBins.begin(), m_previousLoadBins.end(), [](MInt i) { return i > 0; })
       - m_previousLoadBins.begin())
      - 1;
 
  const MInt minFirstBin = (noBinsPrev > 0) ? std::max(std::min(firstBin, firstBinPrev), 0) : std::max(firstBin, 0);
 
  // Set maximum line length
  const MInt maxLineLength = 256;
  const MString dashes(maxLineLength, '-');
  // Create sprintf buffer and string stream
  MChar b[maxLineLength];
  std::stringstream hist;
 
  snprintf(b, maxLineLength, " |%.126s|\n", dashes.c_str());
  MString separatorTop(b);
 
  snprintf(b, maxLineLength, " |%.18s|%.42s|%.42s|%.21s|\n", dashes.c_str(), dashes.c_str(), dashes.c_str(),
           dashes.c_str());
  MString separator(b);
 
  // Create header
  snprintf(b, maxLineLength,
           " | Load distribution at global timestep %-8d%*s - imbalance %5.2f%% - t_avg %5.3e - t_max %5.3e - loads "
           "[%5.3f,%5.3f] |\n",
           globalTimeStep, 2, " ", imbalance * 100.0, averageRunTime, maxRunTime, minLoad, maxLoad);
  hist << "\n" << separatorTop << b;
 
  if(noBinsPrev > 0) {
    snprintf(b, maxLineLength,
             " | Load distribution at global timestep %-8d%*s - imbalance %5.2f%% - t_avg %5.3e - t_max %5.3e - loads "
             "[%5.3f,%5.3f] |\n",
             m_previousLoadInfoStep, 2, " ", m_previousLoadInfo[0], m_previousLoadInfo[1], m_previousLoadInfo[2],
             m_previousLoadInfo[3], m_previousLoadInfo[4]);
    hist << separatorTop << b;
  }
 
  snprintf(b, maxLineLength, " |     load bin     | %-32s%8d | %-32s%8d |   curr/prev count   |\n",
           "current distribution - timestep", globalTimeStep, "previous distribution - timestep",
           m_previousLoadInfoStep);
  hist << separator << b << separator;
 
  MInt checksum = 0, checksumPrev = 0;
  // Output horizontal histogram
  for(MInt i = maxNoBins - 1; i >= minFirstBin; i--) {
    const MInt maxBarWidth = 40;
    const MInt count = (i < noBins) ? loadBins[i] : 0;
    const MInt countPrev = (i < noBinsPrev) ? m_previousLoadBins[i] : 0;
    // Display at least one char if the count is > 0
    const MInt width = (count > 0) ? std::max((MInt)(maxBarWidth * (MFloat)count / (MFloat)maxBinCount), 1) : 0;
    const MInt widthPrev =
        (countPrev > 0) ? std::max((MInt)(maxBarWidth * (MFloat)countPrev / (MFloat)maxBinCount), 1) : 0;
    const MString bar(width, '#');
    const MString barPrev(widthPrev, '@');
 
    if(approx(i * binWidth, 1.0, 0.0001)) {
      // Add horizontal separator marks between overloaded/underloaded ranks
      snprintf(b, maxLineLength, " | [%6.3f, %-6.3f) |_%-40s_|_%-40s_| %8d | %8d |\n", i * binWidth, (i + 1) * binWidth,
               bar.c_str(), barPrev.c_str(), count, countPrev);
    } else {
      snprintf(b, maxLineLength, " | [%6.3f, %-6.3f) | %-40s | %-40s | %8d | %8d |\n", i * binWidth, (i + 1) * binWidth,
               bar.c_str(), barPrev.c_str(), count, countPrev);
    }
    hist << b;
 
    // Make sure all ranks appear in the histogram
    checksum += count;
    checksumPrev += countPrev;
  }
  hist << separator;
 
  if(checksum != noDomains() || (checksumPrev != noDomains() && noBinsPrev > 0)) {
    m_log << "ERROR in load histogram" << std::endl;
    cerr0 << "ERROR in load histogram" << std::endl;
  } else {
    m_log << hist.str() << std::endl;
    cerr0 << hist.str() << std::endl;
  }
 
  // Store current information for next output
  m_previousLoadBins = loadBins;
  m_previousLoadInfoStep = globalTimeStep;
  m_previousLoadInfo[0] = imbalance * 100.0;
  m_previousLoadInfo[1] = averageRunTime;
  m_previousLoadInfo[2] = maxRunTime;
  m_previousLoadInfo[3] = minLoad;
  m_previousLoadInfo[4] = maxLoad;
 
  return loadBalance;
}

◆ noCouplers()

template<MInt nDim>

MInt maia::grid::Controller< nDim >::noCouplers ( ) const

inlineprivate

Definition at line 146 of file cartesiangridcontroller.h.

146{ return m_couplers->size(); }

◆ noDomains()

template<MInt nDim>

MInt maia::grid::Controller< nDim >::noDomains ( )

inlineprivate

Definition at line 135 of file cartesiangridcontroller.h.

135{ return gridb().noDomains(); }

CartesianGrid::noDomains

MInt noDomains() const

Definition: cartesiangrid.h:1054

◆ noSolvers()

template<MInt nDim>

MInt maia::grid::Controller< nDim >::noSolvers ( ) const

inlineprivate

Definition at line 140 of file cartesiangridcontroller.h.

140{ return m_solvers->size(); }

maia::grid::Controller::m_solvers

const std::vector< std::unique_ptr< Solver > > * m_solvers

Definition: cartesiangridcontroller.h:202

◆ partition()

template<MInt nDim>

void maia::grid::Controller< nDim >::partition	(	MLong *	partitionCellOffsets,
		MLong *	globalIdOffsets,
		const MBool	onlyPartitionOffsets
	)

private

Definition at line 1836 of file cartesiangridcontroller.h.

                                                                   {
  TRACE();
  if(gridb().m_maxPartitionLevelShift > 0) {
    cerr0 << "WARNING: Load balancing with partition level shifts might not fully work!" << std::endl;
  }
 
  // Note: noPartitionCells = 0 allowed after updatePartitionCells()!
  const MLong noPartitionCells = gridb().m_noPartitionCells;
 
  // Update noOffsprings and workLoad for all cells and calculate accumulated workload
  MFloat totalWorkload = 0.0;
  MFloatScratchSpace partitionCellsWorkload(std::max(noPartitionCells, 1L), AT_, "partitionCellsWorkload");
  MLongScratchSpace partitionCellsGlobalId(std::max(noPartitionCells, 1L), AT_, "partitionCellsGlobalId");
 
  gridb().calculateNoOffspringsAndWorkload(static_cast<Collector<void>*>(nullptr), gridb().treeb().size());
 
  for(MUint i = 0; i < noPartitionCells; i++) {
    const MLong globalCellId = gridb().m_localPartitionCellGlobalIds[i]; // sorted by global id
    const MInt cellId = gridb().globalIdToLocalId(globalCellId, true);
 
    partitionCellsWorkload(i) = gridb().a_workload(cellId);
    partitionCellsGlobalId(i) = gridb().a_globalId(cellId);
    ASSERT(globalCellId == partitionCellsGlobalId(i),
           "global cell id mismatch! " + std::to_string(globalCellId) + " " + std::to_string(partitionCellsGlobalId(i))
               + " " + std::to_string(cellId));
 
    ASSERT(i == 0 || partitionCellsGlobalId(i) > partitionCellsGlobalId(i - 1),
           "Partition cells not sorted by global id: " + std::to_string(partitionCellsGlobalId(i))
               + " <= " + std::to_string(partitionCellsGlobalId(i - 1)));
    totalWorkload += gridb().a_workload(cellId);
  }
 
  MPI_Allreduce(MPI_IN_PLACE, &totalWorkload, 1, MPI_DOUBLE, MPI_SUM, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "totalWorkload");
 
  MBool calcGlobalOffsets = false;
 
  if(gridb().m_partitionParallelSplit) {
    maia::grid::partitionParallelSplit(&partitionCellsWorkload[0], noPartitionCells,
                                       gridb().m_localPartitionCellOffsets[0], static_cast<MLong>(noDomains()),
                                       static_cast<MLong>(domainId()), gridb().mpiComm(), &partitionCellOffsets[0]);
    partitionCellOffsets[noDomains()] = gridb().m_noPartitionCellsGlobal;
 
    if(!onlyPartitionOffsets) {
      // Determine the globalIdOffsets later (required for load balancing mode 0)
      calcGlobalOffsets = true;
    }
  } else {
    gridb().partitionParallel(gridb().m_noPartitionCells, gridb().m_localPartitionCellOffsets[0],
                              &partitionCellsWorkload[0], &partitionCellsGlobalId[0], totalWorkload,
                              partitionCellOffsets, globalIdOffsets, true);
 
    if(!onlyPartitionOffsets && gridb().m_maxPartitionLevelShift > 0) {
      // Redetermine the global offsets below in case of a partition level shift since
      // partitionParallel() does not include the correction of the offsets as done in
      // correctDomainOffsetsAtPartitionLevelShifts()
      calcGlobalOffsets = true;
    }
  }
 
  // Determine the globalIdOffsets
  if(calcGlobalOffsets) {
    MLongScratchSpace localPartitionCellCounts(noDomains(), AT_, "localPartitionCellCounts");
    MLongScratchSpace localPartitionCellOffsets(noDomains() + 1, AT_, "localPartitionCellOffsets");
    // Determine the number of partition cells on all domains and the partition cell offsets
    gridb().determineNoPartitionCellsAndOffsets(&localPartitionCellCounts[0], &localPartitionCellOffsets[0]);
 
    // Calculate new global domain offsets from the new partition cell offsets
    loadBalancingCalcNewGlobalOffsets(&localPartitionCellOffsets[0], &partitionCellOffsets[0], &globalIdOffsets[0]);
  }
 
  if(domainId() == 0) std::cerr << std::endl;
}

◆ printDomainStatistics()

template<MInt nDim>

void maia::grid::Controller< nDim >::printDomainStatistics ( const MString & status = "" )

Definition at line 1914 of file cartesiangridcontroller.h.

                                                                  {
  TRACE();
 
  MIntScratchSpace minNoSolverCells(noSolvers(), FUN_, "minNoSolverCells");
  MIntScratchSpace maxNoSolverCells(noSolvers(), FUN_, "maxNoSolverCells");
  MIntScratchSpace avgNoSolverCells(noSolvers(), FUN_, "avgNoSolverCells");
  MLongScratchSpace globalNoSolverCells(noSolvers(), FUN_, "globalNoSolverCells");
 
  MLong noGridLeafCells = 0;
  minNoSolverCells.fill(0);
  maxNoSolverCells.fill(0);
  avgNoSolverCells.fill(0);
  globalNoSolverCells.fill(0);
 
  for(MInt cellId = 0; cellId < gridb().treeb().size(); cellId++) {
    if(!gridb().a_hasProperty(cellId, Cell::IsHalo) && gridb().a_noChildren(cellId) == 0) {
      noGridLeafCells++;
    }
  }
 
  MLong globalGridNoLeafCells = noGridLeafCells;
  MLong globalGridNoCells = (MLong)gridb().noInternalCells();
  MInt maxGridNoCells = gridb().noInternalCells();
  MInt minGridNoCells = gridb().noInternalCells();
 
  for(MInt i = 0; i < noSolvers(); i++) {
    minNoSolverCells(i) = solver(i).noInternalCells();
    maxNoSolverCells(i) = solver(i).noInternalCells();
    globalNoSolverCells(i) = solver(i).noInternalCells();
  }
 
  MPI_Allreduce(MPI_IN_PLACE, &globalGridNoLeafCells, 1, MPI_LONG, MPI_SUM, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "globalGridNoLeafCells");
  MPI_Allreduce(MPI_IN_PLACE, &globalGridNoCells, 1, MPI_LONG, MPI_SUM, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "globalGridNoCells");
  MPI_Allreduce(MPI_IN_PLACE, &maxGridNoCells, 1, MPI_INT, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "maxGridNoCells");
  MPI_Allreduce(MPI_IN_PLACE, &minGridNoCells, 1, MPI_INT, MPI_MIN, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "minGridNoCells");
  MPI_Allreduce(MPI_IN_PLACE, &maxNoSolverCells[0], noSolvers(), MPI_INT, MPI_MAX, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "maxNoSolverCells[0]");
  MPI_Allreduce(MPI_IN_PLACE, &minNoSolverCells[0], noSolvers(), MPI_INT, MPI_MIN, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "minNoSolverCells[0]");
  MPI_Allreduce(MPI_IN_PLACE, &globalNoSolverCells[0], noSolvers(), MPI_LONG, MPI_SUM, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "globalNoSolverCells[0]");
 
  for(MInt i = 0; i < noSolvers(); i++) {
    avgNoSolverCells(i) = (MInt)(((MFloat)globalNoSolverCells(i)) / ((MFloat)noDomains()));
  }
 
  const MInt avgGridNoLeafCells = (MInt)(((MFloat)globalGridNoLeafCells) / ((MFloat)noDomains()));
  const MInt avgGridNoCells = (MInt)(((MFloat)globalGridNoCells) / ((MFloat)noDomains()));
  MInt devNoGridCells = std::abs(gridb().noInternalCells() - avgGridNoCells);
  MInt devNoGridLeafCells = std::abs(noGridLeafCells - avgGridNoLeafCells);
 
 
  MPI_Allreduce(MPI_IN_PLACE, &devNoGridCells, 1, MPI_INT, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "devNoGridCells");
  MPI_Allreduce(MPI_IN_PLACE, &devNoGridLeafCells, 1, MPI_INT, MPI_MAX, gridb().mpiComm(), AT_, "MPI_IN_PLACE",
                "devNoGridLeafCells");
  if(domainId() == 0) {
    std::cerr << "Domain statistics ";
    if(!status.empty()) std::cerr << status << " ";
    std::cerr << "at global time step " << globalTimeStep << ": avg no cells=" << avgGridNoCells
              << ", min=" << minGridNoCells << ", max=" << maxGridNoCells << ", max deviation=" << devNoGridCells
              << ", max deviation leaf=" << devNoGridLeafCells << ", total=" << globalGridNoCells << std::endl;
    for(MInt i = 0; i < noSolvers(); i++) {
      std::cerr << "Solver-statistics for solver " << i << " min= " << minNoSolverCells(i)
                << " max= " << maxNoSolverCells(i) << " total= " << globalNoSolverCells(i)
                << " avg= " << avgNoSolverCells(i) << std::endl;
    }
  }
}

◆ redistributeDataDlb()

template<MInt nDim>

void maia::grid::Controller< nDim >::redistributeDataDlb	(	const MInt	id,
		const MInt	mode,
		std::vector< MInt > &	sendSizeVector,
		std::vector< MInt > &	recvSizeVector,
		const MInt *const	bufferIdToCellId,
		const MInt	noCells,
		std::vector< MInt * > &	intDataRecv,
		std::vector< MLong * > &	longDataRecv,
		std::vector< MFloat * > &	floatDataRecv,
		std::vector< MInt > &	dataTypes
	)

private

Definition at line 2951 of file cartesiangridcontroller.h.

                                                                       {
  TRACE();
 
  const MBool isSolver = (mode == 0); // mode 0: solver; mode 1: coupler
  const MInt timerId = (isSolver) ? id : noSolvers() + id;
  // Solver pointer for easier access
  Solver* const solverP = (isSolver) ? &solver(id) : nullptr;
  Coupling* const couplerP = (isSolver) ? nullptr : &coupler(id);
 
  RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::RedistributeMpiBlocking]);
 
  // Determine solver local root domain
  // TODO labels:CONTROLLER,DLB @ansgar_coupler check this!
  const MInt localRootGlobalDomain = (isSolver) ? solverLocalRootDomain(solverP) : 0;
 
  // Get the number of quantities that need to be communicated
  MInt dataCount = (isSolver) ? solverP->noCellDataDlb() : couplerP->noCellDataDlb();
 
  // Note: since inactive ranks might not return the correct count, take the max among all domains
  MPI_Allreduce(MPI_IN_PLACE, &dataCount, 1, type_traits<MInt>::mpiType(), MPI_MAX, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "dataCount");
  RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::RedistributeMpiBlocking]);
 
  MInt intDataCount = 0;
  MInt longDataCount = 0;
  MInt floatDataCount = 0;
 
  // Get solver/solver data and exchange
  for(MInt dataId = 0; dataId < dataCount; dataId++) {
    const MInt offset = dataId * (globalNoDomains() + 1);
    // Inquire data type (communicate from local rank 0 to avoid errors with inactive ranks) and
    // total data send size for this quantity
    // MInt dataType = (solverP->domainId() == 0) ? solverP->cellDataTypeDlb(dataId) : -1;
    MInt dataType = -1;
    if(isSolver) {
      dataType = (solverP->domainId() == 0) ? solverP->cellDataTypeDlb(dataId) : -1;
    } else {
      // TODO labels:CONTROLLER,DLB @ansgar_coupler
      /* dataType = (couplerP->domainId() == 0) ? couplerP->cellDataTypeDlb(dataId) : -1; */
      dataType = couplerP->cellDataTypeDlb(dataId);
    }
 
    RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
    MPI_Bcast(&dataType, 1, type_traits<MInt>::mpiType(), localRootGlobalDomain, gridb().mpiComm(), AT_, "dataType");
    RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
 
    const MInt dataSize = dataSendSize[offset + globalNoDomains()];
    const MInt recvSize = dataRecvSize[offset + globalNoDomains()];
 
    // Check data type (MInt, MLong or MFloat for now)
    switch(dataType) {
      case MINT: {
        dataTypes.push_back(MINT);
        // Allocate (persistent) data receive buffer and store pointer
        MInt* newIntData = nullptr;
        mAlloc(newIntData, std::max(recvSize, 1), "newIntData", AT_);
        intDataRecv.push_back(newIntData);
 
        // Data send buffer, temporary since not needed anymore after exchange
        ScratchSpace<MInt> intDataSend(std::max(dataSize, 1), AT_, "intDataSend");
 
        if(dataSize > 0) {
          // Get the solver/solver data and store in send buffer (sorted according to buffer mapping)
          if(isSolver) {
            solverP->getCellDataDlb(dataId, oldNoCells, &bufferIdToCellId[0], &intDataSend[0]);
          } else {
            couplerP->getCellDataDlb(dataId, oldNoCells, &bufferIdToCellId[0], &intDataSend[0]);
          }
        }
 
        RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
        // Exchange
        maia::mpi::exchangeData(&intDataSend[0], globalDomainId(), globalNoDomains(), gridb().mpiComm(), 1,
                                &dataSendSize[offset], &dataRecvSize[offset], &intDataRecv[intDataCount][0]);
        RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
 
        intDataCount++;
        break;
      }
      case MLONG: {
        dataTypes.push_back(MLONG);
        // Allocate (persistent) data receive buffer and store pointer
        MLong* newLongData = nullptr;
        mAlloc(newLongData, std::max(recvSize, 1), "newLongData", AT_);
        longDataRecv.push_back(newLongData);
 
        // Data send buffer, temporary since not needed anymore after exchange
        ScratchSpace<MLong> longDataSend(std::max(dataSize, 1), AT_, "longDataSend");
 
        if(dataSize > 0) {
          // Get the solver/solver data and store in send buffer (sorted according to buffer mapping)
          if(isSolver) {
            solverP->getCellDataDlb(dataId, oldNoCells, &bufferIdToCellId[0], &longDataSend[0]);
          } else {
            couplerP->getCellDataDlb(dataId, oldNoCells, &bufferIdToCellId[0], &longDataSend[0]);
          }
        }
 
        RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
        // Exchange
        maia::mpi::exchangeData(&longDataSend[0], globalDomainId(), globalNoDomains(), gridb().mpiComm(), 1,
                                &dataSendSize[offset], &dataRecvSize[offset], &longDataRecv[longDataCount][0]);
        RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
 
        longDataCount++;
        break;
      }
      case MFLOAT: {
        dataTypes.push_back(MFLOAT);
        // Allocate (persistent) data receive buffer and store pointer
        MFloat* newFloatData = nullptr;
        mAlloc(newFloatData, std::max(recvSize, 1), "newFloatData", AT_);
        floatDataRecv.push_back(newFloatData);
 
        // Data send buffer, temporary since not needed anymore after exchange
        ScratchSpace<MFloat> floatDataSend(std::max(dataSize, 1), AT_, "floatDataSend");
 
        if(dataSize > 0) {
          // Get the solver/solver data and store in send buffer
          if(isSolver) {
            solverP->getCellDataDlb(dataId, oldNoCells, &bufferIdToCellId[0], &floatDataSend[0]);
          } else {
            couplerP->getCellDataDlb(dataId, oldNoCells, &bufferIdToCellId[0], &floatDataSend[0]);
          }
        }
 
        RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
        // Exchange
        maia::mpi::exchangeData(&floatDataSend[0], globalDomainId(), globalNoDomains(), gridb().mpiComm(), 1,
                                &dataSendSize[offset], &dataRecvSize[offset], &floatDataRecv[floatDataCount][0]);
        RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::RedistributeMpi]);
 
        floatDataCount++;
        break;
      }
      default:
        TERMM(1, "Unknown data type: " + std::to_string(dataType) + ".");
        break;
    }
  }
}

◆ refineCellVec()

template<MInt nDim>

std::vector< std::function< void(const MInt)> > maia::grid::Controller< nDim >::refineCellVec ( )

inlineprivate

Definition at line 345 of file cartesiangridcontroller.h.

                                                           {
    std::vector<std::function<void(const MInt)>> vec(noSolvers());
    for(MInt i = 0; i < noSolvers(); i++) {
      vec[i] = std::bind(&Solver::refineCell, &solver(i), std::placeholders::_1);
    }
    return vec;
  }

◆ removeCellVec()

template<MInt nDim>

std::vector< std::function< void(const MInt)> > maia::grid::Controller< nDim >::removeCellVec ( )

inlineprivate

Definition at line 381 of file cartesiangridcontroller.h.

                                                           {
    std::vector<std::function<void(const MInt)>> vec(noSolvers());
    for(MInt i = 0; i < noSolvers(); i++) {
      vec[i] = std::bind(&Solver::removeCell, &solver(i), std::placeholders::_1);
    }
    return vec;
  }

◆ removeChildsVec()

template<MInt nDim>

std::vector< std::function< void(const MInt)> > maia::grid::Controller< nDim >::removeChildsVec ( )

inlineprivate

Definition at line 352 of file cartesiangridcontroller.h.

                                                             {
    std::vector<std::function<void(const MInt)>> vec(noSolvers());
    for(MInt i = 0; i < noSolvers(); i++) {
      vec[i] = std::bind(&Solver::removeChilds, &solver(i), std::placeholders::_1);
    }
    return vec;
  }

◆ resetAllTimer()

template<MInt nDim>

void maia::grid::Controller< nDim >::resetAllTimer

private

Definition at line 3313 of file cartesiangridcontroller.h.

                                     {
  // Reset timer records
  maia::dlb::g_dlbTimerController.resetRecords();
 
  m_dlbPreviousLocalRunTime = 0.0;
  m_dlbPreviousLocalIdleTime = 0.0;
 
  // ... clear timings ...
  std::vector<MFloat>().swap(m_dlbTimings);
 
  m_dlbLastResetTimeStep = globalTimeStep;
 
  m_log << "Resetting DLB timers at timestep " << globalTimeStep << std::endl;
  cerr0 << "Resetting DLB timers at timestep " << globalTimeStep << std::endl;
}

◆ resizeGridMapVec()

template<MInt nDim>

std::vector< std::function< void()> > maia::grid::Controller< nDim >::resizeGridMapVec ( )

inlineprivate

Definition at line 374 of file cartesiangridcontroller.h.

                                                    {
    std::vector<std::function<void()>> vec(noSolvers());
    for(MInt i = 0; i < noSolvers(); i++) {
      vec[i] = std::bind(&Solver::resizeGridMap, &solver(i));
    }
    return vec;
  }

◆ savePartitionFile()

template<MInt nDim>

void maia::grid::Controller< nDim >::savePartitionFile ( )

inline

Definition at line 108 of file cartesiangridcontroller.h.

                           {
    if(!m_balance || !gridb().wasBalancedAtLeastOnce()) {
      return;
    }
    gridb().savePartitionFile();
  }

◆ setDataDlb()

template<MInt nDim>

void maia::grid::Controller< nDim >::setDataDlb	(	const MInt	solverId,
		const MInt	mode,
		std::vector< MInt * > &	intDataRecv,
		std::vector< MLong * > &	longDataRecv,
		std::vector< MFloat * > &	floatDataRecv,
		std::vector< MInt > &	dataTypes,
		const MBool	freeMemory
	)

private

Definition at line 3100 of file cartesiangridcontroller.h.

                                                                                      {
  TRACE();
 
  const MBool isSolver = (mode == 0); // mode 0: solver; mode 1: coupler
  const MInt timerId = (isSolver) ? id : noSolvers() + id;
  // Pointer for easier access
  Solver* const solverP = (isSolver) ? &solver(id) : nullptr;
  Coupling* const couplerP = (isSolver) ? nullptr : &coupler(id);
 
  // Get the number of quantities that need to be communicated
  MInt dataCount = (isSolver) ? solverP->noCellDataDlb() : couplerP->noCellDataDlb();
 
  RECORD_TIMER_START(m_solverTimers[timerId][SolverTimers::SetDataMpiBlocking]);
  // Note: since inactive ranks might return the correct count, take the max among all domains
  MPI_Allreduce(MPI_IN_PLACE, &dataCount, 1, type_traits<MInt>::mpiType(), MPI_MAX, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "dataCount");
  RECORD_TIMER_STOP(m_solverTimers[timerId][SolverTimers::SetDataMpiBlocking]);
 
  MInt intDataCount = 0;
  MInt longDataCount = 0;
  MInt floatDataCount = 0;
 
  // Set solver/solver data
  for(MInt dataId = 0; dataId < dataCount; dataId++) {
    const MInt dataType = dataTypes[dataId];
 
    switch(dataType) {
      case MINT: {
        if(isSolver) {
          solverP->setCellDataDlb(dataId, &intDataRecv[intDataCount][0]);
        } else {
          couplerP->setCellDataDlb(dataId, &intDataRecv[intDataCount][0]);
        }
        if(freeMemory) {
          mDeallocate(intDataRecv[intDataCount]);
        }
        intDataCount++;
        break;
      }
      case MLONG: {
        if(isSolver) {
          solverP->setCellDataDlb(dataId, &longDataRecv[longDataCount][0]);
        } else {
          couplerP->setCellDataDlb(dataId, &longDataRecv[longDataCount][0]);
        }
        if(freeMemory) {
          mDeallocate(longDataRecv[longDataCount]);
        }
        longDataCount++;
        break;
      }
      case MFLOAT: {
        if(isSolver) {
          solverP->setCellDataDlb(dataId, &floatDataRecv[floatDataCount][0]);
        } else {
          couplerP->setCellDataDlb(dataId, &floatDataRecv[floatDataCount][0]);
        }
        if(freeMemory) {
          mDeallocate(floatDataRecv[floatDataCount]);
        }
        floatDataCount++;
        break;
      }
      default:
        TERMM(1, "Unknown data type.");
        break;
    }
  }
}

◆ solver() [1/2]

template<MInt nDim>

Solver & maia::grid::Controller< nDim >::solver ( const MInt solverId )

inlineprivate

Definition at line 138 of file cartesiangridcontroller.h.

138{ return *m_solvers->at(solverId); }

◆ solver() [2/2]

template<MInt nDim>

const Solver & maia::grid::Controller< nDim >::solver ( const MInt solverId ) const

inlineprivate

Definition at line 139 of file cartesiangridcontroller.h.

139{ return (const Solver&)m_solvers->at(solverId); }

◆ solverLocalRootDomain()

template<MInt nDim>

MInt maia::grid::Controller< nDim >::solverLocalRootDomain ( Solver *const solver )

private

Definition at line 3289 of file cartesiangridcontroller.h.

                                                                 {
  MInt localRootGlobalDomain = (solver->domainId() == 0) ? globalDomainId() : 0;
  MPI_Allreduce(MPI_IN_PLACE, &localRootGlobalDomain, 1, type_traits<MInt>::mpiType(), MPI_SUM, gridb().mpiComm(), AT_,
                "MPI_IN_PLACE", "localRootGlobalDomain");
  return localRootGlobalDomain;
}

◆ storeLoadsAndWeights()

template<MInt nDim>

void maia::grid::Controller< nDim >::storeLoadsAndWeights	(	const MFloat *const	loads,
		const MInt	noLoadTypes,
		const MInt *const	loadQuantities,
		const MFloat	domainWeight,
		const MFloat *const	weights
	)

private

Definition at line 3242 of file cartesiangridcontroller.h.

                                                                         {
  TRACE();
  using namespace maia::parallel_io;
 
  RECORD_TIMER_START(m_timers[Timers::IO]);
 
  if(domainId() == 0) {
    std::cerr << "Store loads and weights" << std::endl;
  }
 
  std::stringstream fileName;
  fileName << m_outputDir << "loads_" << noDomains() << "_" << std::setw(8) << std::setfill('0') << globalTimeStep;
  fileName << ParallelIo::fileExt();
 
  ParallelIo file(fileName.str(), PIO_REPLACE, gridb().mpiComm());
 
  file.defineArray(PIO_FLOAT, "weights", noLoadTypes);
  file.defineArray(PIO_FLOAT, "loads", noDomains());
  file.defineArray(PIO_FLOAT, "domainWeights", noDomains());
 
  // Dimensions of load quantities: noDomains x noLoadTypes
  ParallelIo::size_type dimSizes[] = {noDomains(), noLoadTypes};
  file.defineArray(PIO_INT, "loadQuantities", 2, &dimSizes[0]);
 
  // root writes estimated weights
  if(domainId() == 0) {
    file.setOffset(noLoadTypes, 0);
  } else {
    file.setOffset(0, 0);
  }
  file.writeArray(&weights[0], "weights");
 
  file.setOffset(1, domainId());
  file.writeArray(&loads[domainId()], "loads");
  file.writeArray(&domainWeight, "domainWeights");
 
  // Write load quantities of all domains
  file.setOffset(1, domainId(), 2);
  file.writeArray(&loadQuantities[0], "loadQuantities");
 
  RECORD_TIMER_STOP(m_timers[Timers::IO]);
}

◆ storeTimings()

template<MInt nDim>

void maia::grid::Controller< nDim >::storeTimings

private

Author: Ansgar Niemoeller (ansgar) a.nie.nosp@m.moel.nosp@m.ler@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Definition at line 3178 of file cartesiangridcontroller.h.

                                    {
  TRACE();
  using namespace maia::parallel_io;
 
  RECORD_TIMER_START(m_timers[Timers::IO]);
 
  // Number of timings
  const MInt noTimings = m_dlbTimeStepsAll.size();
  const MInt noValues = 2; // runTime and idleTime
 
  if(noTimings == 0) {
    RECORD_TIMER_STOP(m_timers[Timers::IO]);
    return;
  }
 
  // Timings output file name
  std::stringstream fileName;
  fileName << m_outputDir << "timings_" << noDomains() << "_" << globalTimeStep << ParallelIo::fileExt();
 
  m_log << "Write timings to file: " << fileName.str() << std::endl;
 
  ParallelIo file(fileName.str(), PIO_REPLACE, gridb().mpiComm());
 
  file.defineArray(PIO_INT, "timeStep", noTimings);
 
  // Dimensions of timings: noDomains x noTimings x noValues
  ParallelIo::size_type dimSizes[] = {noDomains(), noTimings, noValues};
 
  file.defineArray(PIO_FLOAT, "timings", 3, &dimSizes[0]);
 
  file.setAttribute("domain index", "dim_0", "timings");
  file.setAttribute("time step index", "dim_1", "timings");
  file.setAttribute("timings index", "dim_2", "timings");
 
  file.setAttribute("run time", "var_0", "timings");
  file.setAttribute("idle time", "var_1", "timings");
 
  // Assemble timesteps and timings
  MIntScratchSpace timeSteps(noTimings, AT_, "timeSteps");
  MFloatScratchSpace data(noTimings, noValues, AT_, "data");
  for(MInt i = 0; i < noTimings; i++) {
    timeSteps[i] = m_dlbTimeStepsAll[i];
    data(i, 0) = m_dlbRunTimeAll[i];
    data(i, 1) = m_dlbIdleTimeAll[i];
  }
 
  // root writes timestep data
  if(domainId() == 0) {
    file.setOffset(noTimings, 0);
  } else {
    file.setOffset(0, 0);
  }
  file.writeArray(&timeSteps[0], "timeStep");
 
  // Write timings of all domains
  file.setOffset(1, domainId(), 3);
  file.writeArray(&data[0], "timings");
 
  RECORD_TIMER_STOP(m_timers[Timers::IO]);
}

◆ swapProxyVec()

template<MInt nDim>

std::vector< std::function< void(const MInt, const MInt)> > maia::grid::Controller< nDim >::swapProxyVec ( )

inlineprivate

Definition at line 359 of file cartesiangridcontroller.h.

                                                                      {
    std::vector<std::function<void(const MInt, const MInt)>> vec(noSolvers());
    for(MInt i = 0; i < noSolvers(); i++) {
      vec[i] = std::bind(&Solver::swapProxy, &solver(i), std::placeholders::_1, std::placeholders::_2);
    }
    return vec;
  }

◆ updateGridPartitionWorkloads()

template<MInt nDim>

MBool maia::grid::Controller< nDim >::updateGridPartitionWorkloads

Definition at line 2312 of file cartesiangridcontroller.h.

                                                     {
  TRACE();
 
  // Check if update of partition workloads is enabled
  MBool updateGridPartitionWorkloads = false;
  updateGridPartitionWorkloads =
      Context::getBasicProperty<MBool>("updateGridPartitionWorkloads", AT_, &updateGridPartitionWorkloads);
  if(!updateGridPartitionWorkloads) {
    return false;
  }
 
  m_log << "Updating partition cell workloads and save them to grid... " << std::endl;
 
  MBool restore = false;
  restore = Context::getBasicProperty<MBool>("restoreDefaultWorkloads", AT_, &restore);
 
  // Number of load quantities (should be the same for every rank)
  const MInt noLoadTypes = globalNoLoadTypes();
  std::vector<MFloat> weights{};
 
  if(noLoadTypes < 1 || restore) {
    restore = true;
    m_log << "Restoring default weights" << std::endl;
  } else {
    m_log << "Using specified solver weights" << std::endl;
    // Get the weights for all solvers
    getSpecifiedSolverWeights(weights);
  }
 
  // Use the determined weights to update the partition workloads (or restore default weighting)
  updateWeightsAndWorkloads(weights, restore);
  // Save new partition cell workloads to grid file
  gridb().savePartitionCellWorkloadsGridFile();
  m_log << "done" << std::endl;
 
  if(domainId() == 0) {
    std::cout << "Updated cell weights. Restart solver with 'updateGridPartitionWorkloads = false'" << std::endl;
  }
  return true;
}

◆ updateWeightsAndWorkloads()

template<MInt nDim>

void maia::grid::Controller< nDim >::updateWeightsAndWorkloads	(	const std::vector< MFloat > &	weights,
		const MBool	restoreDefaultWeights
	)

private

Definition at line 2277 of file cartesiangridcontroller.h.

                                                                                    {
  TRACE();
 
  const MInt noGridCells = gridb().treeb().size();
 
  // Restore default cell weights: each cell has weight 1.
  // else reset all weights to 0
  const MFloat initWeight = (restoreDefaultWeights) ? 1.0 : 0.0;
  for(MInt gridCellId = 0; gridCellId < noGridCells; gridCellId++) {
    gridb().a_weight(gridCellId) = initWeight;
  }
 
  if(!restoreDefaultWeights) {
    // Update the cell weights
    MInt weightOffset = 0;
    for(MInt solverId = 0; solverId < noSolvers(); solverId++) {
      if(solver(solverId).isActive()) { // Only add cell loads for active solvers
        for(MInt cellId = 0; cellId < noGridCells; cellId++) {
          if(gridb().a_hasProperty(cellId, Cell::IsHalo)) continue;
          gridb().a_weight(cellId) += solver(solverId).getCellLoad(cellId, &weights[weightOffset]);
        }
      }
      weightOffset += solver(solverId).noLoadTypes();
    }
  }
 
  // Calculate the partition-cell workloads
  gridb().calculateNoOffspringsAndWorkload(static_cast<Collector<void>*>(nullptr), noGridCells);
}

◆ writeRestartFile()

template<MInt nDim>

void maia::grid::Controller< nDim >::writeRestartFile	(	const MBool	forceRestart,
		const MBool	writeBackup
	)

Definition at line 2392 of file cartesiangridcontroller.h.

                                                                                         {
  if(m_grid == nullptr) {
    return;
  }
  MBoolScratchSpace writeRestart(noSolvers(), FUN_, "writeRestart");
  MBoolScratchSpace writeGridRestart(noSolvers(), FUN_, "writeGridRestart");
 
  // disable load timers:
  std::vector<MInt> dlbTimersStatus(noSolvers());
 
  for(MInt i = 0; i < noSolvers(); i++) {
    dlbTimersStatus[i] = solver(i).dlbTimersEnabled();
    if(dlbTimersStatus[i] != 0) {
      solver(i).disableDlbTimers(); // Disable DLB timers if active
    }
  }
 
  for(MInt i = 0; i < noSolvers(); i++) {
    // prepareRestart() and reInitAfterRestart must be called for all ranks!
    writeRestart[i] = solver(i).prepareRestart(forceRestart, writeGridRestart[i]);
  }
 
  MBool allsolversrestart = false;
  MBool gridrestart = false;
 
  for(MInt i = 0; i < noSolvers(); i++) {
    if(writeRestart[i]) allsolversrestart = true;
    if(writeRestart[i] && writeGridRestart[i]) gridrestart = true;
  }
 
  gridb().deletePeriodicConnection(true);
 
  if(gridrestart || (allsolversrestart && m_saveGridNewPartitionCells)) {
    accumulateCellWeights();
 
    {
      std::stringstream s;
      s << "restartGrid";
      if(writeBackup) s << "Backup";
      if(!m_useNonSpecifiedRestartFile || writeBackup) s << "_00" << globalTimeStep;
      // TODO labels:CONTROLLER,IO,totest change to the version below, fix testcases
      // if (!m_useNonSpecifiedRestartFile || writeBackup) {
      //  s << "_" << setw(8) << setfill('0') << globalTimeStep;
      //}
      s << ParallelIo::fileExt();
      m_currentGridFileName = s.str();
 
      // Cancel any open MPI requests in the solver since these might conflict with communication
      // necessary when writing the grid file
      for(MInt i = 0; i < noSolvers(); i++) {
        solver(i).cancelMpiRequests();
      }
 
      cerr0 << "Saving adapted grid file for all solvers " << s.str() << "... ";
 
      gridb().saveGrid((m_outputDir + m_currentGridFileName).c_str(), m_recalcIds);
 
      cerr0 << "ok." << std::endl;
 
      m_saveGridNewPartitionCells = false;
 
    }
  }
 
  for(MInt c = 0; c < noCouplers(); c++) {
    coupler(c).writeRestartFile(allsolversrestart);
  }
 
  if(allsolversrestart) {
    cerr0 << "Writing restart files at time step " << globalTimeStep << " :" << std::endl;
 
    for(MInt i = 0; i < noSolvers(); i++) {
      if(solver(i).isActive()) {
        solver(i).writeRestartFile(writeRestart[i], writeBackup, m_currentGridFileName, m_recalcIds);
      }
    }
    // savePartitionFile();
  }
 
  for(MInt i = 0; i < noSolvers(); i++) {
    solver(i).reIntAfterRestart(writeRestart[i]);
  }
 
  gridb().restorePeriodicConnection();
 
  // reanable all previously running timers:
  for(MInt i = 0; i < noSolvers(); i++) {
    if(dlbTimersStatus[i] != 0) {
      solver(i).enableDlbTimers(); // Reenable previously active DLB timers
    }
  }
}

Member Data Documentation

◆ m_adaptation

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_adaptation = false

private

Definition at line 213 of file cartesiangridcontroller.h.

◆ m_adaptationInterval

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_adaptationInterval = 0

private

Definition at line 214 of file cartesiangridcontroller.h.

◆ m_adaptationStart

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_adaptationStart = -1

private

Definition at line 215 of file cartesiangridcontroller.h.

◆ m_adaptationStop

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_adaptationStop = -1

private

Definition at line 216 of file cartesiangridcontroller.h.

◆ m_balance

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_balance = false

private

Definition at line 225 of file cartesiangridcontroller.h.

◆ m_balanceAdaptationInterval

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_balanceAdaptationInterval = 1

private

Definition at line 324 of file cartesiangridcontroller.h.

◆ m_balanceAfterAdaptation

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_balanceAfterAdaptation = false

private

Definition at line 323 of file cartesiangridcontroller.h.

◆ m_cellOutside

template<MInt nDim>

const std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)> > maia::grid::Controller< nDim >::m_cellOutside

private

Definition at line 206 of file cartesiangridcontroller.h.

◆ m_couplers

template<MInt nDim>

const std::vector<std::unique_ptr<Coupling> >* maia::grid::Controller< nDim >::m_couplers

private

Definition at line 210 of file cartesiangridcontroller.h.

◆ m_currentGridFileName

template<MInt nDim>

MString maia::grid::Controller< nDim >::m_currentGridFileName

private

Definition at line 390 of file cartesiangridcontroller.h.

◆ m_debugBalance

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_debugBalance = false

private

Definition at line 231 of file cartesiangridcontroller.h.

◆ m_dlbIdleTimeAll

template<MInt nDim>

std::vector<MFloat> maia::grid::Controller< nDim >::m_dlbIdleTimeAll

private

Definition at line 300 of file cartesiangridcontroller.h.

◆ m_dlbImbalanceThreshold

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_dlbImbalanceThreshold = 0.05

private

Definition at line 266 of file cartesiangridcontroller.h.

◆ m_dlbLastResetTimeStep

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_dlbLastResetTimeStep = 0

private

Definition at line 263 of file cartesiangridcontroller.h.

◆ m_dlbLastWeights

template<MInt nDim>

MFloat* maia::grid::Controller< nDim >::m_dlbLastWeights = nullptr

private

Definition at line 328 of file cartesiangridcontroller.h.

◆ m_dlbMaxWorkloadLimit

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_dlbMaxWorkloadLimit = 1.5

private

Definition at line 286 of file cartesiangridcontroller.h.

◆ m_dlbNoFinalLocalShifts

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_dlbNoFinalLocalShifts = 0

private

Definition at line 284 of file cartesiangridcontroller.h.

◆ m_dlbNoLocalShifts

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_dlbNoLocalShifts = 0

private

Definition at line 281 of file cartesiangridcontroller.h.

◆ m_dlbPartitionMethod

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_dlbPartitionMethod = DLB_PARTITION_DEFAULT

private

Definition at line 270 of file cartesiangridcontroller.h.

◆ m_dlbPreviousLocalIdleTime

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_dlbPreviousLocalIdleTime = 0.0

private

Definition at line 290 of file cartesiangridcontroller.h.

◆ m_dlbPreviousLocalRunTime

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_dlbPreviousLocalRunTime = 0.0

private

Definition at line 289 of file cartesiangridcontroller.h.

◆ m_dlbRestartWeights

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_dlbRestartWeights = false

private

Definition at line 327 of file cartesiangridcontroller.h.

◆ m_dlbRunTimeAll

template<MInt nDim>

std::vector<MFloat> maia::grid::Controller< nDim >::m_dlbRunTimeAll

private

Definition at line 299 of file cartesiangridcontroller.h.

◆ m_dlbSmoothGlobalShifts

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_dlbSmoothGlobalShifts = true

private

Definition at line 279 of file cartesiangridcontroller.h.

◆ m_dlbStaticWeightMode

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_dlbStaticWeightMode = -1

private

Definition at line 330 of file cartesiangridcontroller.h.

◆ m_dlbStaticWeights

template<MInt nDim>

MFloat* maia::grid::Controller< nDim >::m_dlbStaticWeights = nullptr

private

Definition at line 329 of file cartesiangridcontroller.h.

◆ m_dlbStep

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_dlbStep = 0

private

Definition at line 258 of file cartesiangridcontroller.h.

◆ m_dlbTimeStepsAll

template<MInt nDim>

std::vector<MInt> maia::grid::Controller< nDim >::m_dlbTimeStepsAll

private

Definition at line 298 of file cartesiangridcontroller.h.

◆ m_dlbTimings

template<MInt nDim>

std::vector<MFloat> maia::grid::Controller< nDim >::m_dlbTimings

private

Definition at line 293 of file cartesiangridcontroller.h.

◆ m_dlbUpdatePartitionCells

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_dlbUpdatePartitionCells = false

private

Definition at line 274 of file cartesiangridcontroller.h.

◆ m_domainWeights

template<MInt nDim>

std::vector<MFloat> maia::grid::Controller< nDim >::m_domainWeights

private

Definition at line 308 of file cartesiangridcontroller.h.

◆ m_forceBalance

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_forceBalance = 0

private

Definition at line 227 of file cartesiangridcontroller.h.

◆ m_forceLoadBalancing

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_forceLoadBalancing = false

private

Definition at line 244 of file cartesiangridcontroller.h.

◆ m_grid

template<MInt nDim>

Grid* maia::grid::Controller< nDim >::m_grid

private

Definition at line 200 of file cartesiangridcontroller.h.

◆ m_imbalance

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_imbalance = std::numeric_limits<MFloat>::max()

private

Definition at line 320 of file cartesiangridcontroller.h.

◆ m_lastAdaptationTimeStep

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_lastAdaptationTimeStep = 0

private

Definition at line 217 of file cartesiangridcontroller.h.

◆ m_lastLoadBalancingTimeStep

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_lastLoadBalancingTimeStep = 0

private

Definition at line 260 of file cartesiangridcontroller.h.

◆ m_lastOffsetShiftDirection

template<MInt nDim>

std::vector<MInt> maia::grid::Controller< nDim >::m_lastOffsetShiftDirection {}

private

Definition at line 313 of file cartesiangridcontroller.h.

◆ m_loadBalancingInterval

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_loadBalancingInterval = -1

private

Definition at line 239 of file cartesiangridcontroller.h.

◆ m_loadBalancingMode

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_loadBalancingMode = -1

private

Definition at line 234 of file cartesiangridcontroller.h.

◆ m_loadBalancingOffset

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_loadBalancingOffset = 0

private

Definition at line 240 of file cartesiangridcontroller.h.

◆ m_loadBalancingStartTimeStep

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_loadBalancingStartTimeStep = 0

private

Definition at line 241 of file cartesiangridcontroller.h.

◆ m_loadBalancingStopTimeStep

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_loadBalancingStopTimeStep = -1

private

Definition at line 242 of file cartesiangridcontroller.h.

◆ m_loadBalancingTimerMode

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_loadBalancingTimerMode = 0

private

Definition at line 236 of file cartesiangridcontroller.h.

◆ m_loadBalancingTimerStartOffset

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_loadBalancingTimerStartOffset = 0

private

Definition at line 243 of file cartesiangridcontroller.h.

◆ m_maxPerformanceVarThreshold

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_maxPerformanceVarThreshold = 0.15

private

Definition at line 268 of file cartesiangridcontroller.h.

◆ m_nAdaptationsSinceBalance

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_nAdaptationsSinceBalance = 9999999

private

Definition at line 252 of file cartesiangridcontroller.h.

◆ m_noAdaptations

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_noAdaptations = 0

private

Definition at line 218 of file cartesiangridcontroller.h.

◆ m_noMaxAdaptations

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_noMaxAdaptations = -1

private

Definition at line 219 of file cartesiangridcontroller.h.

◆ m_noSensors

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_noSensors {}

private

Definition at line 222 of file cartesiangridcontroller.h.

◆ m_optPartitionCellOffsetTotal

template<MInt nDim>

MLong maia::grid::Controller< nDim >::m_optPartitionCellOffsetTotal = -1

private

Definition at line 319 of file cartesiangridcontroller.h.

◆ m_outputDir

template<MInt nDim>

MString maia::grid::Controller< nDim >::m_outputDir

private

Definition at line 392 of file cartesiangridcontroller.h.

◆ m_outputDlbTimings

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_outputDlbTimings = false

private

Definition at line 296 of file cartesiangridcontroller.h.

◆ m_performanceOutput

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_performanceOutput = false

private

Definition at line 229 of file cartesiangridcontroller.h.

◆ m_previousLoadBins

template<MInt nDim>

std::vector<MInt> maia::grid::Controller< nDim >::m_previousLoadBins {}

private

Definition at line 303 of file cartesiangridcontroller.h.

◆ m_previousLoadInfo

template<MInt nDim>

std::array<MFloat, 5> maia::grid::Controller< nDim >::m_previousLoadInfo {}

private

Definition at line 304 of file cartesiangridcontroller.h.

◆ m_previousLoadInfoStep

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_previousLoadInfoStep = -1

private

Definition at line 305 of file cartesiangridcontroller.h.

◆ m_recalcIds

template<MInt nDim>

MInt* maia::grid::Controller< nDim >::m_recalcIds

private

Definition at line 393 of file cartesiangridcontroller.h.

◆ m_refineCellSolver

template<MInt nDim>

const std::vector<std::function<void(const MInt)> > maia::grid::Controller< nDim >::m_refineCellSolver

private

Definition at line 203 of file cartesiangridcontroller.h.

◆ m_removeCellSolver

template<MInt nDim>

const std::vector<std::function<void(const MInt)> > maia::grid::Controller< nDim >::m_removeCellSolver

private

Definition at line 208 of file cartesiangridcontroller.h.

◆ m_removeChildsSolver

template<MInt nDim>

const std::vector<std::function<void(const MInt)> > maia::grid::Controller< nDim >::m_removeChildsSolver

private

Definition at line 204 of file cartesiangridcontroller.h.

◆ m_resizeGridMapSolver

template<MInt nDim>

const std::vector<std::function<void()> > maia::grid::Controller< nDim >::m_resizeGridMapSolver

private

Definition at line 207 of file cartesiangridcontroller.h.

◆ m_saveGridNewPartitionCells

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_saveGridNewPartitionCells = false

private

Definition at line 276 of file cartesiangridcontroller.h.

◆ m_solvers

template<MInt nDim>

const std::vector<std::unique_ptr<Solver> >* maia::grid::Controller< nDim >::m_solvers

private

Definition at line 202 of file cartesiangridcontroller.h.

◆ m_solverTimerGroups

template<MInt nDim>

std::vector<MInt> maia::grid::Controller< nDim >::m_solverTimerGroups {}

private

Definition at line 340 of file cartesiangridcontroller.h.

◆ m_solverTimers

template<MInt nDim>

std::vector<std::array<MInt, SolverTimers::_count> > maia::grid::Controller< nDim >::m_solverTimers {}

private

Definition at line 342 of file cartesiangridcontroller.h.

◆ m_swapProxySolver

template<MInt nDim>

const std::vector<std::function<void(const MInt, const MInt)> > maia::grid::Controller< nDim >::m_swapProxySolver

private

Definition at line 205 of file cartesiangridcontroller.h.

◆ m_syncTimerSteps

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_syncTimerSteps = true

private

Definition at line 255 of file cartesiangridcontroller.h.

◆ m_syncTimeSteps

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_syncTimeSteps = false

private

Definition at line 254 of file cartesiangridcontroller.h.

◆ m_testDynamicLoadBalancing

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_testDynamicLoadBalancing = false

private

Definition at line 245 of file cartesiangridcontroller.h.

◆ m_testUpdatePartCellsOffspringThreshold

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_testUpdatePartCellsOffspringThreshold = 10

private

Definition at line 249 of file cartesiangridcontroller.h.

◆ m_testUpdatePartCellsWorkloadThreshold

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_testUpdatePartCellsWorkloadThreshold = 100.0

private

Definition at line 250 of file cartesiangridcontroller.h.

◆ m_testUpdatePartitionCells

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_testUpdatePartitionCells = false

private

Definition at line 247 of file cartesiangridcontroller.h.

◆ m_timePerStepTotal

template<MInt nDim>

MFloat maia::grid::Controller< nDim >::m_timePerStepTotal = std::numeric_limits<MFloat>::max()

private

Definition at line 318 of file cartesiangridcontroller.h.

◆ m_timerGroup

template<MInt nDim>

MInt maia::grid::Controller< nDim >::m_timerGroup = -1

private

Definition at line 336 of file cartesiangridcontroller.h.

◆ m_timers

template<MInt nDim>

std::array<MInt, Timers::_count> maia::grid::Controller< nDim >::m_timers {}

private

Definition at line 338 of file cartesiangridcontroller.h.

◆ m_timersInitialized

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_timersInitialized = false

private

Definition at line 334 of file cartesiangridcontroller.h.

◆ m_useDomainFactor

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_useDomainFactor = false

private

Definition at line 310 of file cartesiangridcontroller.h.

◆ m_useNonSpecifiedRestartFile

template<MInt nDim>

MBool maia::grid::Controller< nDim >::m_useNonSpecifiedRestartFile = false

private

Definition at line 391 of file cartesiangridcontroller.h.

The documentation for this class was generated from the following files:

/home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/GRID/cartesiangrid.h
/home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/GRID/cartesiangridcontroller.h
/home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/GRID/cartesiangridcontroller.cpp

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ Cell

◆ Grid

◆ GridProxy

Constructor & Destructor Documentation

◆ Controller()

◆ ~Controller()

Member Function Documentation

◆ accumulateCellWeights()

◆ adaptation()

◆ balance()

◆ castAdaptationIntervalToMultipleOfCoarsestTimeStep()

◆ cellOutsideVec()

◆ communicateGlobalSolverVars()

◆ computeWeights()

◆ coupler() [1/2]

◆ coupler() [2/2]

◆ determineDataSizesDlb()

◆ domainId()

◆ estimateParameters()

◆ getLoadQuantities()

◆ getSpecifiedSolverWeights()

◆ globalNoLoadTypes()

◆ gridb() [1/2]

◆ gridb() [2/2]

◆ initDlbProperties()

◆ initTimers()

◆ isAdaptationTimeStep()

◆ isDlbTimeStep()

◆ loadBalancingCalcNewGlobalOffsets()

◆ loadBalancingCalcNoCellsToSend()

◆ loadBalancingPartition()

◆ logTimerStatistics()

◆ needLoadBalancing()

◆ noCouplers()

◆ noDomains()

◆ noSolvers()

◆ partition()

◆ printDomainStatistics()

◆ redistributeDataDlb()

◆ refineCellVec()

◆ removeCellVec()

◆ removeChildsVec()

◆ resetAllTimer()

◆ resizeGridMapVec()

◆ savePartitionFile()

◆ setDataDlb()

◆ solver() [1/2]

◆ solver() [2/2]

◆ solverLocalRootDomain()

◆ storeLoadsAndWeights()

◆ storeTimings()

◆ swapProxyVec()

◆ updateGridPartitionWorkloads()

◆ updateWeightsAndWorkloads()

◆ writeRestartFile()

Member Data Documentation

◆ m_adaptation

◆ m_adaptationInterval

◆ m_adaptationStart

◆ m_adaptationStop

◆ m_balance

◆ m_balanceAdaptationInterval

◆ m_balanceAfterAdaptation

◆ m_cellOutside

◆ m_couplers

◆ m_currentGridFileName

◆ m_debugBalance

◆ m_dlbIdleTimeAll

◆ m_dlbImbalanceThreshold

◆ m_dlbLastResetTimeStep

◆ m_dlbLastWeights

◆ m_dlbMaxWorkloadLimit

◆ m_dlbNoFinalLocalShifts

◆ m_dlbNoLocalShifts