ParallelCoarseSpaceMatVec.h 11.9 KB
Newer Older
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
// ============================================================================
// ==                                                                        ==
// == 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

26
#include <mpi.h>
27 28 29 30
#include <vector>
#include <map>
#include <petsc.h>
#include "AMDiS_fwd.h"
31 32
#include "parallel/ParallelDofMapping.h"
#include "parallel/MeshDistributor.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
33
#include "parallel/MatrixNnzStructure.h"
34 35 36 37 38 39 40 41 42 43 44 45 46 47

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.
48 49
   * - This class also manages the creation of the corresponding non zero 
   *   structure of the matrices.
50 51 52
   */
  class ParallelCoarseSpaceMatVec {
  public:
53
    /// Constructor
54
    ParallelCoarseSpaceMatVec();
55

56 57 58 59 60
    /// Set parallel DOF mapping for the interior DOFs.
    void setDofMapping(ParallelDofMapping *interiorDofs)
    {
      interiorMap = interiorDofs;
    }
61 62 63 64 65 66

    /// Returns the parallel DOF mapping for the interior DOFs.
    ParallelDofMapping* getDofMapping()
    {
      return interiorMap;
    }
67
    
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
    /** \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<DOFMatrix*>& seqMat);
99

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
100 101
    /// Run PETSc's matrix assembly routines.
    void matAssembly();
102

103
    /// Run PETSc's vector assembly routines on rhs vectors.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
104 105
    void vecRhsAssembly();

106
    /// Run PETSc's vector assembly routines on solution vectors.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
107 108
    void vecSolAssembly();

109
    /// Destroys PETSc matrix objects.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
110 111
    void matDestroy();

112
    /// Destroys PETSc vector objects.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
113
    void vecDestroy();
114

115 116
    /// Get interior matrix.
    inline Mat& getMatInterior()
117 118 119 120 121
    {
      TEST_EXIT_DBG(mat.size() > 0)("No matrix data!\n");
      return mat[0][0];
    }

122 123
    /// Get coarse space matrix.
    inline Mat& getMatCoarse(int coarseSpace0 = 0, int coarseSpace1 = 0)
124 125 126 127 128 129
    {
      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];
    }

130 131
    /// Get coupling matrix of the interior and some coarse space.
    inline Mat& getMatInteriorCoarse(int coarseSpace = 0)
132 133 134 135 136
    {
      TEST_EXIT_DBG(mat.size() > coarseSpace + 1)("No matrix data!\n");
      return mat[0][coarseSpace + 1];
    }

137 138
    /// Get coupling of some coarse space matrix and the interior.
    inline Mat& getMatCoarseInterior(int coarseSpace = 0)
139 140 141 142 143
    {
      TEST_EXIT_DBG(mat.size() > coarseSpace + 1)("No matrix data!\n");
      return mat[coarseSpace + 1][0];
    }

144
    /// Get the coarse space matrix of some system component.
Thomas Witkowski's avatar
Thomas Witkowski committed
145
    inline Mat& getMatCoarseByComponent(int rowComp, int colComp = -1)
146
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
147 148 149
      int rowMatIndex = componentIthCoarseMap[rowComp] + 1;
      int colMatIndex = componentIthCoarseMap[(colComp == -1 ? rowComp : colComp)] + 1;
      return mat[rowMatIndex][colMatIndex];
150 151
    }

152 153 154
    /// Get coupling matrix of the interior and the coarse space of a 
    /// system component.
    inline Mat& getMatInteriorCoarseByComponent(int comp)
155 156 157 158 159
    {
      int matIndex = componentIthCoarseMap[comp] + 1;
      return mat[0][matIndex];
    }

160 161 162
    /// Get coupling matrix of the coarse space of a system component and the
    /// interior matrix.
    inline Mat& getMatCoarseInteriorByComponent(int comp)
163 164 165 166 167
    {
      int matIndex = componentIthCoarseMap[comp] + 1;
      return mat[matIndex][0];
    }

168 169
    /// Get the RHS vector of the interior.
    inline Vec& getVecRhsInterior()
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
170 171 172 173
    {
      return vecRhs[0];
    }

174 175
    /// Get the RHS vector of some coarse space.
    inline Vec& getVecRhsCoarse(int coarseSpace = 0)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
176 177 178 179
    {
      return vecRhs[coarseSpace + 1];
    }

180 181 182 183 184 185 186
    /// 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];
    }

187 188
    /// Get the solution vector of the interior.
    inline Vec& getVecSolInterior()
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
189 190 191 192
    {
      return vecSol[0];
    }

193 194
    /// Get the solution vector of some coarse space.
    inline Vec& getVecSolCoarse(int coarseSpace = 0)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
195
    {
196 197 198 199 200
      FUNCNAME("ParallelCoarseSpaceMatVec::getVecSolCoarse()");

      TEST_EXIT_DBG(coarseSpace + 1 < vecSol.size())
	("Wrong component %d, vecSol has only %d!\n", coarseSpace + 1, vecSol.size());

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
201 202 203
      return vecSol[coarseSpace + 1];
    }

204 205 206 207 208 209 210
    /// 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];
    }

211 212 213 214 215 216 217 218 219 220 221 222 223
    /** \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]  dof        DOF index
     *
     * \return     True, if the dof is a coarse space DOF in the component. 
     *             False otherwise.
     */
    inline bool isCoarseSpace(int component,
			      DegreeOfFreedom dof)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
224
    {
225 226 227 228 229 230 231 232
      FUNCNAME("ParallelCoarseSpaceMatVec::isCoarseSpace()");
      
      if (coarseSpaceMap.empty())
	return false;

      TEST_EXIT_DBG(coarseSpaceMap.count(component))
	("Component %d has no coarse space defined!\n", component);

233
      return (*(coarseSpaceMap[component]))[component].isSet(dof);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
234 235
    }

236

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
237
  protected:
238 239 240 241 242 243 244
    /// 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's avatar
Blub  
Thomas Witkowski committed
245 246
    void updateSubdomainData();

247
  private:
248 249 250 251 252 253 254 255 256 257 258
    /// 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.
259 260
    vector<vector<Mat> > mat;

261 262
    /// Solution and RHS vectors. vec[0] is the interior vector, vec[1] the 
    /// first coarse space vector and so on.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
263 264
    vector<Vec> vecSol, vecRhs;

265 266
    /// Matrix of objects to control the matrix non zero structure of the 
    /// corresponding PETSc matrices stored in \ref mat.
267 268
    vector<vector<MatrixNnzStructure> > nnz;

269 270 271 272 273 274 275 276 277 278 279 280 281
    /** \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.
     */
282
    vector<int> componentIthCoarseMap;
283

284
    /// Stores a set of all coarse space DOF mapping. All entries are unique.
285 286
    vector<ParallelDofMapping*> uniqueCoarseMap;

287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302
    /// 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.
303 304
    MeshDistributor *meshDistributor;

305 306
    /// Level of subdomain/interior discretization. Is used for multi-level
    /// methods only.
307 308
    int subdomainLevel;

309 310 311
    /// 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.
312 313
    MPI::Intracomm mpiCommLocal;
    
314 315
    /// 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.
316 317
    MPI::Intracomm mpiCommGlobal;

318 319 320 321 322
    /// 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's avatar
Blub  
Thomas Witkowski committed
323 324
    int rStartInterior;

325 326 327
    /// 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's avatar
Blub  
Thomas Witkowski committed
328 329
    int nGlobalOverallInterior;

330 331
    /// Parallel DOF mapping for the interior.
    ParallelDofMapping *interiorMap;
332

333 334 335
    /// Parallel DOF mapping of the (optional) coarse space. Allows to define
    /// different coarse spaces for different components.
    map<int, ParallelDofMapping*> coarseSpaceMap;   
336 337 338 339
  };
}

#endif