// ============================================================================ // == == // == AMDiS - Adaptive multidimensional simulations == // == == // == http://www.amdis-fem.org == // == == // ============================================================================ // // Software License for AMDiS // // Copyright (c) 2010 Dresden University of Technology // All rights reserved. // Authors: Simon Vey, Thomas Witkowski et al. // // This file is part of AMDiS // // See also license.opensource.txt in the distribution. /** \file MeshDistributor.h */ #ifndef AMDIS_MESHDISTRIBUTOR_H #define AMDIS_MESHDISTRIBUTOR_H #include #include "parallel/DofComm.h" #include "parallel/ElementObjectDatabase.h" #include "parallel/ParallelTypes.h" #include "parallel/MeshLevelData.h" #include "parallel/MeshPartitioner.h" #include "parallel/InteriorBoundary.h" #include "parallel/ParallelDofMapping.h" #include "parallel/PeriodicMap.h" #include "parallel/StdMpi.h" #include "AMDiS_fwd.h" #include "Containers.h" #include "Global.h" #include "ProblemTimeInterface.h" #include "ProblemIterationInterface.h" #include "FiniteElemSpace.h" #include "Serializer.h" #include "BoundaryManager.h" #include "SystemVector.h" namespace AMDiS { using namespace std; struct BoundaryDofInfo { map geoDofs; }; class MeshDistributor { private: MeshDistributor(); public: ~MeshDistributor(); void initParallelization(); void exitParallelization(); /// Adds a DOFVector to the set of \ref interchangeVecs. Thus, this vector /// will be automatically interchanged between ranks when mesh is /// repartitioned. void addInterchangeVector(DOFVector *vec) { interchangeVectors.push_back(vec); } /// Adds all DOFVectors of a SystemVector to \ref interchangeVecs. void addInterchangeVector(SystemVector *vec) { for (int i = 0; i < vec->getSize(); i++) interchangeVectors.push_back(vec->getDOFVector(i)); } /** \brief * This function checks if the mesh has changed on at least on rank. In * this case, the interior boundaries are adapted on all ranks such that * they fit together on all ranks. Furthermore the function * \ref updateLocalGlobalNumbering() is called to update the DOF numberings * and mappings on all rank due to the new mesh structure. * * \param[in] tryRepartition If this parameter is true, repartitioning * may be done. This depends on several other * parameters. If the parameter is false, the * mesh is only checked and adapted but never * repartitioned. */ void checkMeshChange(bool tryRepartition = true); /// Checks if is required to repartition the mesh. If this is the case, a new /// partition will be created and the mesh will be redistributed between the /// ranks. void repartitionMesh(); void getImbalanceFactor(double &imbalance, int &minDofs, int &maxDofs, int &sumDofs); double getImbalanceFactor(); /// Calculates the imbalancing factor and prints it to screen. void printImbalanceFactor(); /// Test, if the mesh consists of macro elements only. The mesh partitioning /// of the parallelization works for macro meshes only and would fail, if the /// mesh is already refined in some way. Therefore, this function will exit /// the program if it finds a non macro element in the mesh. void testForMacroMesh(); /// Set for each element on the partitioning level the number of /// leaf elements. void setInitialElementWeights(); inline string getName() { return name; } inline Mesh* getMesh() { return mesh; } /// Returns an FE space from \ref feSpaces. inline const FiniteElemSpace* getFeSpace(unsigned int i = 0) { FUNCNAME("MeshDistributor::getFeSpace()"); TEST_EXIT_DBG(i < feSpaces.size())("Should not happen!\n"); return feSpaces[i]; } /// Returns all FE spaces, thus \ref feSpaces. inline vector& getFeSpaces() { return feSpaces; } /// Returns the DOF mapping object, \ref dofMap. inline ParallelDofMapping& getDofMap() { return dofMap; } inline ParallelDofMapping& getDofMapSd() { return dofMapSd; } /// Returns the periodic mapping handler, \ref periodicMap. inline PeriodicMap& getPeriodicMap() { return periodicMap; } DofComm& getDofComm() { return dofComm; } DofComm& getDofCommSd() { return dofCommSd; } InteriorBoundary& getIntBoundary() { return intBoundary; } InteriorBoundary& getIntBoundarySd() { return intBoundarySd; } inline long getLastMeshChangeIndex() { return lastMeshChangeIndex; } inline int getMpiRank() { return mpiRank; } inline int getMpiSize() { return mpiSize; } inline MPI::Intracomm& getMpiComm() { return mpiComm; } inline bool isInitialized() { return initialized; } // Writes all data of this object to an output stream. void serialize(ostream &out); // Reads the object data from an input stream. void deserialize(istream &in); /** \brief * This function must be used if the values of a DOFVector must be * synchronised over all ranks. That means, that each rank sends the * values of the DOFs, which are owned by the rank and lie on an interior * bounday, to all other ranks also having these DOFs. * * This function must be used, for example, after the lineary system is * solved, or after the DOFVector is set by some user defined functions, * e.g., initial solution functions. */ template void synchVector(DOFVector &vec) { StdMpi > stdMpi(mpiComm); const FiniteElemSpace *fe = vec.getFeSpace(); for (DofComm::Iterator it(dofComm.getSendDofs(), fe); !it.end(); it.nextRank()) { vector dofs; dofs.reserve(it.getDofs().size()); for (; !it.endDofIter(); it.nextDof()) dofs.push_back(vec[it.getDofIndex()]); stdMpi.send(it.getRank(), dofs); } for (DofComm::Iterator it(dofComm.getRecvDofs()); !it.end(); it.nextRank()) stdMpi.recv(it.getRank()); stdMpi.startCommunication(); for (DofComm::Iterator it(dofComm.getRecvDofs(), fe); !it.end(); it.nextRank()) for (; !it.endDofIter(); it.nextDof()) vec[it.getDofIndex()] = stdMpi.getRecvData(it.getRank())[it.getDofCounter()]; } /// Works in the same way as the function above defined for DOFVectors. Due /// to performance, this function does not call \ref synchVector for each /// DOFVector, but instead sends all values of all DOFVectors all at once. void synchVector(SystemVector &vec, int level = 0); /// Works quite similar to the function \ref synchVector, but instead the /// values of subdomain vectors are add along the boundaries. template void synchAddVector(DOFVector &vec) { StdMpi > stdMpi(mpiComm); const FiniteElemSpace *fe = vec.getFeSpace(); for (DofComm::Iterator it(dofComm.getRecvDofs(), fe); !it.end(); it.nextRank()) { vector dofs; dofs.reserve(it.getDofs().size()); for (; !it.endDofIter(); it.nextDof()) dofs.push_back(vec[it.getDofIndex()]); stdMpi.send(it.getRank(), dofs); } for (DofComm::Iterator it(dofComm.getSendDofs()); !it.end(); it.nextRank()) stdMpi.recv(it.getRank()); stdMpi.startCommunication(); for (DofComm::Iterator it(dofComm.getSendDofs(), fe); !it.end(); it.nextRank()) for (; !it.endDofIter(); it.nextDof()) vec[it.getDofIndex()] += stdMpi.getRecvData(it.getRank())[it.getDofCounter()]; synchVector(vec); } /// In 3D, a subdomain may not be a valid AMDiS mesh if it contains two /// parts which are only connected by an edge. In this case, the standard /// refinement algorithm does not work correctly, as two elements connected /// only on one edge are not neighours by definition. This functions checks /// for this situation and fix the problem. For this, the mesh is search for /// all edges connecting two elements that are otherwise not connected. void check3dValidMesh(); void setBoundaryDofRequirement(Flag flag) { createBoundaryDofFlag = flag; } BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace, int level) { FUNCNAME("MeshDistributor::getBoundaryDofInfo()"); TEST_EXIT_DBG(level < static_cast(boundaryDofInfo.size())) ("Wrong level number: %d, whereas array size is %d!\n", level, boundaryDofInfo.size()); return boundaryDofInfo[level][feSpace]; } void getAllBoundaryDofs(const FiniteElemSpace *feSpace, int level, DofContainer& dofs); const ElementObjectDatabase& getElementObjectDb() { return elObjDb; } /// Adds a stationary problem to the global mesh distributor objects. static void addProblemStatGlobal(ProblemStatSeq *probStat); MeshLevelData& getMeshLevelData() { return levelData; } void updateLocalGlobalNumbering(); /// Updates the local and global DOF numbering after the mesh has been /// changed. void updateLocalGlobalNumbering(ParallelDofMapping &dmap, DofComm &dcom, const FiniteElemSpace *feSpace); protected: void addProblemStat(ProblemStatSeq *probStat); /// Determines the interior boundaries, i.e. boundaries between ranks, and /// stores all information about them in \ref interiorBoundary. void createInteriorBoundary(bool firstCall); void createBoundaryDofs(); /// Removes all macro elements from the mesh that are not part of ranks /// partition. void removeMacroElements(); /// Calls \ref createPeriodicMap(feSpace) for all FE spaces that are /// handled by the mesh distributor. void createPeriodicMap(); /// Creates, for a specific FE space, to all DOFs in rank's partition that /// are on a periodic boundary the mapping from dof index to the other /// periodic dof indices. This information is stored in \ref periodicDofMap. void createPeriodicMap(const FiniteElemSpace *feSpace); /// This function is called only once during the initialization when the /// whole macro mesh is available on all cores. It copies the pointers of all /// macro elements to \ref allMacroElements and stores all neighbour /// information based on macro element indices (and not pointer based) in /// \ref macroElementNeighbours. These information are then used to /// reconstruct macro elements during mesh redistribution. void createMacroElementInfo(); void updateMacroElementInfo(); /** \brief * Checks for all given interior boundaries if the elements fit together on * both sides of the boundaries. If this is not the case, the mesh is * adapted. Because refinement of a certain element may forces the * refinement of other elements, it is not guaranteed that all rank's meshes * fit together after this function terminates. Hence, it must be called * until a stable mesh refinement is reached. * * \param[in] allBound Defines a map from rank to interior boundaries * which should be checked. * * \return If the mesh has been changed by this function, it returns * true. Otherwise, it returns false, i.e., the given interior * boundaries fit together on both sides. */ bool checkAndAdaptBoundary(RankToBoundMap &allBound); /// Sets \ref isRankDof to all matrices and rhs vectors in a given /// stationary problem. void setRankDofs(ProblemStatSeq *probStat); /// Sets \ref isRankDof to all matrices and rhs vectors in all /// stationary problems. void setRankDofs(); /// Removes all periodic boundary condition information from all matrices and /// vectors of all stationary problems and from the mesh itself. void removePeriodicBoundaryConditions(); /// Removes all periodic boundary condition information from all matrices and /// vector of a given stationary problem. void removePeriodicBoundaryConditions(ProblemStatSeq *probStat); // Removes all periodic boundaries from a given boundary map. void removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap); void createMeshLevelStructure(); /// Writes a vector of dof pointers to an output stream. void serialize(ostream &out, DofContainer &data); /// Writes a \ref RankToDofContainer to an output stream. void serialize(ostream &out, map > &data); /// Reads a vector of dof pointers from an input stream. void deserialize(istream &in, DofContainer &data, map &dofIndexMap); /// Reads a \ref RankToDofContainer from an input stream. void deserialize(istream &in, map > &data, map > &dofIndexMap); /// Writes a mapping from dof pointers to some values to an output stream. template void serialize(ostream &out, map &data) { FUNCNAME("ParallelDomainBase::serialize()"); int mapSize = data.size(); SerUtil::serialize(out, mapSize); for (typename map::iterator it = data.begin(); it != data.end(); ++it) { int v1 = (*(it->first)); T v2 = it->second; SerUtil::serialize(out, v1); SerUtil::serialize(out, v2); } } /// Reads a mapping from dof pointer to some values from an input stream. template void deserialize(istream &in, map &data, map &dofIndexMap) { FUNCNAME("ParallelDomainBase::deserialize()"); int mapSize = 0; SerUtil::deserialize(in, mapSize); for (int i = 0; i < mapSize; i++) { int v1 = 0; T v2; SerUtil::deserialize(in, v1); SerUtil::deserialize(in, v2); TEST_EXIT_DBG(dofIndexMap.count(v1) != 0) ("Cannot find DOF %d in map!\n", v1); data[dofIndexMap[v1]] = v2; } } protected: /// List of all stationary problems that are managed by this mesh /// distributor. vector problemStat; /// If true, the mesh distributor is already initialized; bool initialized; /// The rank of the current process. int mpiRank; /// Overall number of processes. int mpiSize; /// MPI communicator collected all processes, which should be used for /// calculation. The Debug procces is not included in this communicator. MPI::Intracomm mpiComm; /// Name of the problem (as used in the init files) string name; /// Finite element spaces of the problem. vector feSpaces; /// Mesh of the problem. Mesh *mesh; /// A refinement manager that should be used on the mesh. It is used to /// refine elements at interior boundaries in order to fit together with /// elements on the other side of the interior boundary. RefinementManager *refineManager; /// Info level. int info; /// Pointer to a mesh partitioner that is used to partition the mesh to /// the ranks. MeshPartitioner *partitioner; /// Weights for the elements, i.e., the number of leaf elements within /// this element. map elemWeights; /// Stores to every macro element index the number of the rank that owns this /// macro element. map partitionMap; /// Mapping object to map from local DOF indices to global ones. ParallelDofMapping dofMap; ParallelDofMapping dofMapSd; /// Database to store and query all sub-objects of all elements of the /// macro mesh. ElementObjectDatabase elObjDb; /// Defines the interior boundaries of the domain that result from /// partitioning the whole mesh. InteriorBoundary intBoundary; InteriorBoundary intBoundarySd; DofComm dofComm; DofComm dofCommSd; PeriodicMap periodicMap; /// This set of values must be interchanged between ranks when the mesh is /// repartitioned. vector*> interchangeVectors; /// If the problem definition has been read from a serialization file, this /// variable is true, otherwise it is false. This variable is used to stop the /// initialization function, if the problem definition has already been read /// from a serialization file. bool deserialized; /// Denotes whether there exists a filewriter for this object. bool writeSerializationFile; /// If true, it is possible to repartition the mesh during computations. bool repartitioningAllowed; /// Stores the number of mesh changes that must lie in between to /// repartitionings. int repartitionIthChange; /// Counts the number of mesh changes after the last mesh repartitioning /// was done. int nMeshChangesAfterLastRepartitioning; /// Countes the number of mesh repartitions that were done. Till now, this /// variable is used only for debug outputs. int repartitioningCounter; /// If repartitioning of the mesh fail, this variable has a positive value /// that gives the number of mesh changes the mesh distributer will wait /// before trying new mesh repartitioning. int repartitioningFailed; /// Directory name where all debug output files should be written to. string debugOutputDir; /// Stores the mesh change index. This is used to recognize changes in the /// mesh structure (e.g. through refinement or coarsening managers). long lastMeshChangeIndex; /// Stores for all macro elements of the original macro mesh the /// neighbourhood information based on element indices. Thus, each macro /// element index is mapped to a vector containing all indices of /// neighbouring macro elements. map > macroElementNeighbours; /// Store all macro elements of the overall mesh, i.e., before the /// mesh is redistributed for the first time. vector allMacroElements; Flag createBoundaryDofFlag; /// Stores on each mesh level for all FE spaces the information about /// all boundary DOFs. vector > boundaryDofInfo; MeshLevelData levelData; /// If there is no mesh adaptivity, the mesh distributor can remove some /// data structures which are only used if mesh changes or it must be /// redistributed due to some local adaptivity. By default, this variable /// is set to true, and thus no special assumption are made. bool meshAdaptivity; public: static bool sebastianMode; /// The boundary DOFs are sorted by subobject entities, i.e., first all /// face DOFs, edge DOFs and to the last vertex DOFs will be set to /// communication structure vectors, \ref sendDofs and \ref recvDofs. static const Flag BOUNDARY_SUBOBJ_SORTED; /// When boundary DOFs are created, \ref boundaryDofInfo is filled for /// all DOFs that this rank will send to other ranks (thus, rank /// owned DOFs. static const Flag BOUNDARY_FILL_INFO_SEND_DOFS; /// When boundary DOFs are created, \ref boundaryDofInfo is filled for /// all DOFs that this rank will receive from other ranks (thus, DOFs /// that are owned by another rank). static const Flag BOUNDARY_FILL_INFO_RECV_DOFS; static MeshDistributor *globalMeshDistributor; friend class ParallelDebug; }; } #endif // AMDIS_MESHDISTRIBUTOR_H