ParallelCoarseSpaceMatVec.h 11.6 KB
 Thomas Witkowski committed Aug 20, 2012 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 // ============================================================================ // == == // == 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 ParallelCoarseSpaceMatVec.h */ #ifndef AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H #define AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H #include #include #include #include "AMDiS_fwd.h"  Thomas Witkowski committed Aug 21, 2012 30 #include "parallel/MatrixNnzStructure.h"  Thomas Witkowski committed Aug 20, 2012 31 32 33 34 35 36 37 38 39 40 41 42 43 44  namespace AMDiS { /** * This class implements a block structured PETSc matrix/vec which seperates * the discretization of the interior of subdomains and the discretization * of the coarse space. Thus, we have one matrix block for the interior and * one matrix block for the coarse space plus the coupling blocks. Some notes: * - For a single level domain decomposition method (e.g. the standad * FETI-DP method), the interior matrix is local to the current rank and the * coarse space matrix is a globally distributed matrix. * - There are different coarse spaces for different components possible. In * this case, there are as many blocks as there are different coarse spaces * plus one block for the interior matrix.  45 46  * - This class also manages the creation of the corresponding non zero * structure of the matrices.  Thomas Witkowski committed Aug 20, 2012 47 48 49  */ class ParallelCoarseSpaceMatVec { public:  Thomas Witkowski committed Aug 20, 2012 50  ParallelCoarseSpaceMatVec();  Thomas Witkowski committed Aug 20, 2012 51   52 53 54 55 56  /// Set parallel DOF mapping for the interior DOFs. void setDofMapping(ParallelDofMapping *interiorDofs) { interiorMap = interiorDofs; }  Thomas Witkowski committed Aug 20, 2012 57   58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88  /** \brief * Sets the coarse space for all or a specific component. * * \param[in] coarseDofs Coarse space DOF mapping. * \param[in] component If the standard value -1 is used, the coarse * space DOF mapping is set for all components * of the equation. Otherwise, the coarse space * DOF mapping is set only for the given one. */ void setCoarseSpaceDofMapping(ParallelDofMapping *coarseDofs, int component = -1); /// Set mesh distributor object and MPI communicators. void setMeshDistributor(MeshDistributor *m, MPI::Intracomm mpiComm0, MPI::Intracomm mpiComm1) { meshDistributor = m; mpiCommGlobal = mpiComm0; mpiCommLocal = mpiComm1; } /// Set level of the interior discretization. Is only used for /// multi-level methods. void setLevel(int l) { subdomainLevel = l; } /// Creates matrices and vectors with respect to the coarse space. void createMatVec(Matrix& seqMat);  Thomas Witkowski committed Aug 20, 2012 89   Thomas Witkowski committed Aug 21, 2012 90 91  /// Run PETSc's matrix assembly routines. void matAssembly();  Thomas Witkowski committed Aug 20, 2012 92   93  /// Run PETSc's vector assembly routines on rhs vectors.  Thomas Witkowski committed Aug 21, 2012 94 95  void vecRhsAssembly();  96  /// Run PETSc's vector assembly routines on solution vectors.  Thomas Witkowski committed Aug 21, 2012 97 98  void vecSolAssembly();  99  /// Destroys PETSc matrix objects.  Thomas Witkowski committed Aug 21, 2012 100 101  void matDestroy();  102  /// Destroys PETSc vector objects.  Thomas Witkowski committed Aug 21, 2012 103  void vecDestroy();  Thomas Witkowski committed Aug 20, 2012 104   105 106  /// Get interior matrix. inline Mat& getMatInterior()  Thomas Witkowski committed Aug 20, 2012 107 108 109 110 111  { TEST_EXIT_DBG(mat.size() > 0)("No matrix data!\n"); return mat[0][0]; }  112 113  /// Get coarse space matrix. inline Mat& getMatCoarse(int coarseSpace0 = 0, int coarseSpace1 = 0)  Thomas Witkowski committed Aug 20, 2012 114 115 116 117 118 119  { TEST_EXIT_DBG(mat.size() > coarseSpace0 + 1)("No matrix data!\n"); TEST_EXIT_DBG(mat.size() > coarseSpace1 + 1)("No matrix data!\n"); return mat[coarseSpace0 + 1][coarseSpace1 + 1]; }  120 121  /// Get coupling matrix of the interior and some coarse space. inline Mat& getMatInteriorCoarse(int coarseSpace = 0)  Thomas Witkowski committed Aug 20, 2012 122 123 124 125 126  { TEST_EXIT_DBG(mat.size() > coarseSpace + 1)("No matrix data!\n"); return mat[0][coarseSpace + 1]; }  127 128  /// Get coupling of some coarse space matrix and the interior. inline Mat& getMatCoarseInterior(int coarseSpace = 0)  Thomas Witkowski committed Aug 20, 2012 129 130 131 132 133  { TEST_EXIT_DBG(mat.size() > coarseSpace + 1)("No matrix data!\n"); return mat[coarseSpace + 1][0]; }  134  /// Get the coarse space matrix of some system component.  Thomas Witkowski committed Aug 22, 2012 135  inline Mat& getMatCoarseByComponent(int rowComp, int colComp = -1)  Thomas Witkowski committed Aug 20, 2012 136  {  Thomas Witkowski committed Aug 22, 2012 137 138 139  int rowMatIndex = componentIthCoarseMap[rowComp] + 1; int colMatIndex = componentIthCoarseMap[(colComp == -1 ? rowComp : colComp)] + 1; return mat[rowMatIndex][colMatIndex];  Thomas Witkowski committed Aug 20, 2012 140 141  }  142 143 144  /// Get coupling matrix of the interior and the coarse space of a /// system component. inline Mat& getMatInteriorCoarseByComponent(int comp)  Thomas Witkowski committed Aug 20, 2012 145 146 147 148 149  { int matIndex = componentIthCoarseMap[comp] + 1; return mat[0][matIndex]; }  150 151 152  /// Get coupling matrix of the coarse space of a system component and the /// interior matrix. inline Mat& getMatCoarseInteriorByComponent(int comp)  Thomas Witkowski committed Aug 20, 2012 153 154 155 156 157  { int matIndex = componentIthCoarseMap[comp] + 1; return mat[matIndex][0]; }  158 159  /// Get the RHS vector of the interior. inline Vec& getVecRhsInterior()  Thomas Witkowski committed Aug 21, 2012 160 161 162 163  { return vecRhs[0]; }  164 165  /// Get the RHS vector of some coarse space. inline Vec& getVecRhsCoarse(int coarseSpace = 0)  Thomas Witkowski committed Aug 21, 2012 166 167 168 169  { return vecRhs[coarseSpace + 1]; }  Thomas Witkowski committed Aug 22, 2012 170 171 172 173 174 175 176  /// Get the RHS vector of the coarse space of a system component. inline Vec& getVecRhsCoarseByComponent(int comp) { int vecIndex = componentIthCoarseMap[comp] + 1; return vecRhs[vecIndex]; }  177 178  /// Get the solution vector of the interior. inline Vec& getVecSolInterior()  Thomas Witkowski committed Aug 21, 2012 179 180 181 182  { return vecSol[0]; }  183 184  /// Get the solution vector of some coarse space. inline Vec& getVecSolCoarse(int coarseSpace = 0)  Thomas Witkowski committed Aug 21, 2012 185 186 187 188  { return vecSol[coarseSpace + 1]; }  Thomas Witkowski committed Aug 22, 2012 189 190 191 192 193 194 195  /// Get the solution vector of the coarse space of a system component. inline Vec& getVecSolCoarseByComponent(int comp) { int vecIndex = componentIthCoarseMap[comp] + 1; return vecSol[vecIndex]; }  196 197 198 199 200 201 202 203 204 205 206 207 208 209 210  /** \brief * Checks whether a given DOF index in some component is a coarse space DOF. * Note (TODO): The specification of both, the component number and FE * space is not really necessary. Rewrite this! * * \param[in] component Component number of the system. * \param[in] feSpace Finite element space of the component. * \param[in] dof DOF index * * \return True, if the dof is a coarse space DOF in the component. * False otherwise. */ inline bool isCoarseSpace(int component, const FiniteElemSpace *feSpace, DegreeOfFreedom dof)  Thomas Witkowski committed Aug 21, 2012 211  {  212 213 214 215 216 217 218 219 220  FUNCNAME("ParallelCoarseSpaceMatVec::isCoarseSpace()"); if (coarseSpaceMap.empty()) return false; TEST_EXIT_DBG(coarseSpaceMap.count(component)) ("Component %d has no coarse space defined!\n", component); return (*(coarseSpaceMap[component]))[feSpace].isSet(dof);  Thomas Witkowski committed Aug 21, 2012 221 222  }  223   Thomas Witkowski committed Aug 21, 2012 224  protected:  225 226 227 228 229 230 231  /// Prepare internal data structures. First, it create \ref uniqueCoarseMap /// and \ref componentIthCoarseMap . Both are used to create the correct /// number of matrix and vectors in \ref mat and \ref vec. void prepare(); /// Computes the values of \ref rStartInterior and /// \ref nGlobalOverallInterior.  Thomas Witkowski committed Aug 21, 2012 232 233  void updateSubdomainData();  Thomas Witkowski committed Aug 20, 2012 234  private:  235 236 237 238 239 240 241 242 243 244 245  /// Checks for mesh changes. Returns true if the mesh has been changed /// until the last matrix creation. Is used to rebuild matrix non /// zero stricture. bool checkMeshChange(); private: /// Matrix of PETSc matrices. mat[0][0] is the interior discretization /// matrix, mat[1][1] corresponds to the first coarse space and so on. /// mat[i][j], with i not equal to j, are the coupling between the interior /// and the coarse space matrices, and between the different coarse spaces /// respectively.  Thomas Witkowski committed Aug 20, 2012 246 247  vector > mat;  248 249  /// Solution and RHS vectors. vec[0] is the interior vector, vec[1] the /// first coarse space vector and so on.  Thomas Witkowski committed Aug 21, 2012 250 251  vector vecSol, vecRhs;  252 253  /// Matrix of objects to control the matrix non zero structure of the /// corresponding PETSc matrices stored in \ref mat.  Thomas Witkowski committed Aug 20, 2012 254 255  vector > nnz;  256 257 258 259 260 261 262 263 264 265 266 267 268  /** \brief * Stores for each system component (i.e. each PDE variable) the coarse * space that is used for its discretization. * * Example: We solve the Stokes equation in 2D with a different coarse * space for the velocity unknowns (component 0 and 1) and the pressure * (component 2). Than: * componentIthCoarseMap[0] = 0 * componentIthCoarseMap[1] = 0 * componentIthCoarseMap[2] = 1 * The indices can directly be used to access the correspondig parallel * DOF mapping in \ref uniqueCoarseMap. */  Thomas Witkowski committed Aug 20, 2012 269  vector componentIthCoarseMap;  Thomas Witkowski committed Aug 20, 2012 270   271  /// Stores a set of all coarse space DOF mapping. All entries are unique.  Thomas Witkowski committed Aug 20, 2012 272 273  vector uniqueCoarseMap;  274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289  /// Stores the mesh change index of the mesh the nnz structure was created /// for. Therefore, if the mesh change index is higher than this value, we /// have to create a new nnz structure for PETSc matrices, because the mesh /// has been changed and therefore also the assembled matrix structure. int lastMeshNnz; /// If this variable is set to true, the non-zero matrix structure is /// created each time from scratch by calling \ref createPetscNnzStrcuture. /// This can be necessary if the number of non-zeros in the matrix varies /// though the mesh does not change. This may happen if there are many /// operators using DOFVectors from old timestep containing many zeros due to /// some phase fields. bool alwaysCreateNnzStructure; protected: /// Pointer to a mesh distributor object.  Thomas Witkowski committed Aug 20, 2012 290 291  MeshDistributor *meshDistributor;  292 293  /// Level of subdomain/interior discretization. Is used for multi-level /// methods only.  Thomas Witkowski committed Aug 20, 2012 294 295  int subdomainLevel;  296 297 298  /// MPI communicator on the subdomain level. If no multi-level is used /// this is alway MPI_COMM_SELF. In the case of a multi-level method, this /// is a subset of MPI_COMM_WORLD.  Thomas Witkowski committed Aug 20, 2012 299 300  MPI::Intracomm mpiCommLocal;  301 302  /// MPI communicator on the coarse space level. If no multi-level method /// is used, this is always MPI_COMM_WORLD, otherwise a subset of it.  Thomas Witkowski committed Aug 20, 2012 303 304  MPI::Intracomm mpiCommGlobal;  305 306 307 308 309  /// Offset for the interior DOFs of the local interior with respect to the /// subdomain. In the case of a one-level method, each local interior /// is exactly one subdomain. In the case of a multi-level method, one /// subdomain may consists of several rank domains. This value defines than /// the offset ot rank's interior rows to the subdomain's interior rows.  Thomas Witkowski committed Aug 21, 2012 310 311  int rStartInterior;  312 313 314  /// Number of overall rows in subdomain's interior. For one-level methods, /// this value is equal to the number of rows in rank's interior. See also /// explenation for \ref rStarInterior.  Thomas Witkowski committed Aug 21, 2012 315 316  int nGlobalOverallInterior;  317 318  /// Parallel DOF mapping for the interior. ParallelDofMapping *interiorMap;  Thomas Witkowski committed Aug 20, 2012 319   320 321 322  /// Parallel DOF mapping of the (optional) coarse space. Allows to define /// different coarse spaces for different components. map coarseSpaceMap;  Thomas Witkowski committed Aug 20, 2012 323 324 325 326  }; } #endif