ParallelCoarseSpaceMatVec.h 11.6 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
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 <vector>
#include <map>
#include <petsc.h>
#include "AMDiS_fwd.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
30
#include "parallel/MatrixNnzStructure.h"
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.
47
48
49
   */
  class ParallelCoarseSpaceMatVec {
  public:
50
    ParallelCoarseSpaceMatVec();
51

52
53
54
55
56
    /// Set parallel DOF mapping for the interior DOFs.
    void setDofMapping(ParallelDofMapping *interiorDofs)
    {
      interiorMap = interiorDofs;
    }
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<DOFMatrix*>& seqMat);
89

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
90
91
    /// Run PETSc's matrix assembly routines.
    void matAssembly();
92

93
    /// Run PETSc's vector assembly routines on rhs vectors.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
94
95
    void vecRhsAssembly();

96
    /// Run PETSc's vector assembly routines on solution vectors.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
97
98
    void vecSolAssembly();

99
    /// Destroys PETSc matrix objects.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
100
101
    void matDestroy();

102
    /// Destroys PETSc vector objects.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
103
    void vecDestroy();
104

105
106
    /// Get interior matrix.
    inline Mat& getMatInterior()
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)
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)
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)
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's avatar
Thomas Witkowski committed
135
    inline Mat& getMatCoarseByComponent(int rowComp, int colComp = -1)
136
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
139
      int rowMatIndex = componentIthCoarseMap[rowComp] + 1;
      int colMatIndex = componentIthCoarseMap[(colComp == -1 ? rowComp : colComp)] + 1;
      return mat[rowMatIndex][colMatIndex];
140
141
    }

142
143
144
    /// Get coupling matrix of the interior and the coarse space of a 
    /// system component.
    inline Mat& getMatInteriorCoarseByComponent(int comp)
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)
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's avatar
Blub    
Thomas Witkowski committed
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's avatar
Blub    
Thomas Witkowski committed
166
167
168
169
    {
      return vecRhs[coarseSpace + 1];
    }

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's avatar
Blub    
Thomas Witkowski committed
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's avatar
Blub    
Thomas Witkowski committed
185
186
187
188
    {
      return vecSol[coarseSpace + 1];
    }

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's avatar
Blub    
Thomas Witkowski committed
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's avatar
Blub    
Thomas Witkowski committed
221
222
    }

223

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
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's avatar
Blub    
Thomas Witkowski committed
232
233
    void updateSubdomainData();

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.
246
247
    vector<vector<Mat> > 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's avatar
Blub    
Thomas Witkowski committed
250
251
    vector<Vec> vecSol, vecRhs;

252
253
    /// Matrix of objects to control the matrix non zero structure of the 
    /// corresponding PETSc matrices stored in \ref mat.
254
255
    vector<vector<MatrixNnzStructure> > 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.
     */
269
    vector<int> componentIthCoarseMap;
270

271
    /// Stores a set of all coarse space DOF mapping. All entries are unique.
272
273
    vector<ParallelDofMapping*> 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.
290
291
    MeshDistributor *meshDistributor;

292
293
    /// Level of subdomain/interior discretization. Is used for multi-level
    /// methods only.
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.
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.
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's avatar
Blub    
Thomas Witkowski committed
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's avatar
Blub    
Thomas Witkowski committed
315
316
    int nGlobalOverallInterior;

317
318
    /// Parallel DOF mapping for the interior.
    ParallelDofMapping *interiorMap;
319

320
321
322
    /// Parallel DOF mapping of the (optional) coarse space. Allows to define
    /// different coarse spaces for different components.
    map<int, ParallelDofMapping*> coarseSpaceMap;   
323
324
325
326
  };
}

#endif