KDDecomposition.h

#ifndef KDDECOMPOSITION2_H_
#define KDDECOMPOSITION2_H_
 
#include <mpi.h>
#include <algorithm>
#include <list>
#include <vector>
#include <limits>
 
#define KDDIM 3
 
#include "DomainDecompMPIBase.h"
#include "parallel/CollectiveCommunication.h"
#include "ParticleDataForwardDeclaration.h"
#include "particleContainer/adapter/CellProcessor.h"
#include "parallel/LoadCalc.h"
 
class KDNode;
 
 
class KDDecomposition: public DomainDecompMPIBase {
 
    friend class KDDecompositionTest;
 
 public:
    KDDecomposition(double cutoffRadius, int numParticleTypes, int updateFrequency = 100,  int fullSearchThreshold = 2);
 
    void init(Domain* domain);
 
    ~~KDDecomposition() override;
 
    virtual void readXML(XMLfileUnits& xmlconfig) override;
 
    // documentation in base class
    virtual void prepareNonBlockingStage(bool forceRebalancing,
                ParticleContainer* moleculeContainer, Domain* domain,
                unsigned int stageNumber) override;
 
    // documentation in base class
    virtual void finishNonBlockingStage(bool forceRebalancing,
            ParticleContainer* moleculeContainer, Domain* domain,
            unsigned int stageNumber) override;
 
    // documentation in base class
    bool queryBalanceAndExchangeNonBlocking(bool forceRebalancing, ParticleContainer* moleculeContainer, Domain* domain, double etime) override;
 
    void balanceAndExchange(double lastTraversalTime, bool forceRebalancing, ParticleContainer* moleculeContainer, Domain* domain) override;
 
    double getBoundingBoxMin(int dimension, Domain* domain) override;
    double getBoundingBoxMax(int dimension, Domain* domain) override;
 
    int getUpdateFrequency() const { return _frequency; }
    void setUpdateFrequency(int frequency) { _frequency = frequency; }
    std::vector<int> getNeighbourRanks() override;
    std::vector<int> getNeighbourRanksFullShell() override;
 
    std::vector<CommunicationPartner> getNeighboursFromHaloRegion(Domain* domain, const HaloRegion& haloRegion, double cutoff) override;
 
    /*
     * Create a TunerTimes object with measured values from the tuner
     *
     * Has to be an extra function since the CellProcessor is only created after the KDDecomposition was already constructed
     */
    void fillTimeVecs(CellProcessor **cellProc);
 
    void printTree(std::ostream& ostream);
 
private:
    void constructNewTree(KDNode *& newRoot, KDNode *& newOwnLeaf, ParticleContainer* moleculeContainer);
    bool migrateParticles(const KDNode& newRoot, const KDNode& newOwnLeaf, ParticleContainer* moleculeContainer, Domain* domain);
    void initCommunicationPartners(double cutoffRadius, Domain * domain, ParticleContainer* moleculeContainer);
 
 
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    //$ sonstige Methoden
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 
    void calculateCostsPar(KDNode* area, std::vector<std::vector<double> >& costsLeft, std::vector<std::vector<double> >& costsRight, MPI_Comm commGroup);
 
 
    unsigned int getGlobalIndex(int divDim, int dim1, int dim2, int index_dim, int index_dim1, int index_dim2, KDNode* localArea);
 
    int ownMod(int number, int modulo) const;
 
    void getCellBorderFromIntCoords(double * lC, double * hC, int lo[3], int hi[3]) const;
 
    void getCellIntCoordsFromRegionPeriodic(int * lo, int * hi, const double lC[3], const double hC[3], const Domain* domain) const;
 
    void completeTreeInfo(KDNode*& root, KDNode*& ownArea);
 
 
    void getOwningProcs(int low[KDDIM], int high[KDDIM], KDNode* decompTree, KDNode* testNode, std::vector<int>* procIDs, std::vector<int>* neighbHaloAreas) const;
 
 
    void calcNumParticlesPerCell(ParticleContainer* moleculeContainer);
 
    bool decompose(KDNode* fatherNode, KDNode*& ownArea, MPI_Comm commGroup);
 
    bool decompose(KDNode* fatherNode, KDNode*& ownArea, MPI_Comm commGroup, double globalMinimalDeviation);
 
    bool heteroDecompose(KDNode* fatherNode, KDNode*& ownArea, MPI_Comm commGroup);
 
    bool calculateAllPossibleSubdivisions(KDNode* node, std::list<KDNode*>& subdividedNodes, MPI_Comm commGroup);
    
    bool calculateHeteroSubdivision(KDNode* node, KDNode*& subdivededNode, MPI_Comm commGroup);
 
    void updateMeanProcessorSpeeds(std::vector<double>& processorSpeeds,
            std::vector<double>& accumulatedProcessorSpeeds, ParticleContainer* moleculeContainer);
 
    void collectMoleculesInRegion(ParticleContainer* moleculeContainer, const double startRegion[3], const double endRegion[3], std::vector<Molecule*>& mols) const;
 
    /*
     * Determines the partition rank that is needed for the "cluster" heterogeneous decomposition
     */
    int calculatePartitionRank();
    bool checkNeedRebalance(double lastTraversalTime) const;
 
    //######################################
    //###    private member variables    ###
    //######################################
    double _cellSize[KDDIM]{};
 
    int _globalCellsPerDim[KDDIM]{};
 
    int _globalNumCells{};
 
    // Pointer to the root element of the decomposition tree
    KDNode* _decompTree{nullptr};
 
    // each process owns an area in the decomposition
    KDNode* _ownArea{nullptr};
 
    std::vector<unsigned int> _numParticlesPerCell;
 
    /* TODO: This may not be equal to the number simulation steps if balanceAndExchange
     * is not called exactly once in every simulation step! */
    size_t _steps{0ul};
 
    int _frequency;
 
    double _cutoffRadius;
    int _cellsInCutoffRadius{1};
 
    /*
     * Threshold for full tree search. If a node has more than _fullSearchThreshold processors,
     * it is for each dimension divided in the middle only. Otherwise, all possible subdivisions
     * are created. // TODO hetero: change description
     */
    int _fullSearchThreshold;
 
 
    std::vector<double> _processorSpeeds;
    std::vector<double> _accumulatedProcessorSpeeds;//length nprocs+1, first element is 0.
    double _totalMeanProcessorSpeed{1.};
    double _totalProcessorSpeed{1.};
    int _processorSpeedUpdateCount{0};
    bool _heterogeneousSystems{false};
    bool _clusteredHeterogeneouseSystems{false};
    bool _splitBiggest{true};  // indicates, whether a subdomain is to be split along its biggest size
    bool _forceRatio{false};  // if you want to enable forcing the above ratio, enable this.
 
    bool _doMeasureLoadCalc {false};  // specifies if measureLoad should be used.
    int  _measureLoadInterpolationStartsAt{1};  // specifies at which number of particles per cell measureLoad should start using interpolation.
    bool _measureLoadIncreasingTimeValues{true};  // specifies if the time values should be increasing if the number of particles increases.
 
    int _splitThreshold{std::numeric_limits<int>::max()};
    int _numParticleTypes;
 
    /*
     * Stores the maximum number of particles encountered in a cell by the processor (per particle type).
     * This information is max-reduced in the destructor to print the overall maximum number of particles (per type).
     *
     * This is can be useful for the vectorization tuner to know how many values should be measured
     * and comes at the cost of only one additional conditional move per cell per load balancing.
     */
    int _maxPars{std::numeric_limits<int>::min()};
    int _maxPars2{std::numeric_limits<int>::min()};
 
    /*
     * The partition rank that is used in the "clustered" heterogeneous decomposition to determine
     * the mpi ranks associated with each cluster.
     *
     * Cluster 1: (0; partitionRank(
     * Cluster 2: (partitionRank; maxRank)
     */
    const int _partitionRank;
    LoadCalc* _loadCalc;  // stores the times (and constants) measured by the vectorization tuner
    MeasureLoad* _measureLoadCalc;  // stores the measured times of the real-world simulations
 
    /*
     * The following Variables are only used for as parameters for the Vectorization tuner constructor.
     * They are stored here since the tuner is still partly a optional OutputPlugin which doesn't
     * really fit to its use now as essential part of the load balancing of the k-d-tree.
     *
     * For documentation see public tune function of the vectorization tuner
     *
     * TODO These should probably be transferred into the tuner at some point,
     * when a decision is made whether the tuner should still be a Outputplugin or simply be a part of the
     * KDDecomposition (or both)
     */
    std::vector<int> _vecTunParticleNums;
    bool _generateNewFiles{true};
    bool _useExistingFiles{true};
    bool _vecTunerAllowMPIReduce{true};
 
 
    double _rebalanceLimit{0.}; 
 
    MPI_Op _deviationReductionOperation{MPI_SUM};
};
 
 
#endif /* KDDECOMPOSITION2_H_ */