ParallelCoarseSpaceSolver.h 12.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// ============================================================================
// ==                                                                        ==
// == 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.



21
/** \file ParallelCoarseSpaceSolver.h */
22

23
24
#ifndef AMDIS_PARALLEL_COARSE_SPACE_SOLVER_H
#define AMDIS_PARALLEL_COARSE_SPACE_SOLVER_H
25

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
#include "OEMSolver.h"
35
36

namespace AMDiS {
37
  
38
39
40
41
42
43
44
45
46
47
48
  /**
   * 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.
49
50
   * - This class also manages the creation of the corresponding non zero 
   *   structure of the matrices.
51
   */
52
  class ParallelCoarseSpaceSolver : public OEMSolver {
53
  public:
54
    /// Constructor
55
    ParallelCoarseSpaceSolver(string name);
56

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

    /// Returns the parallel DOF mapping for the interior DOFs.
    ParallelDofMapping* getDofMapping()
    {
      return interiorMap;
    }
68
    
69
70
71
72
73
74
75
76
77
78
79
80
    /** \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);

81
82
83
84
85
    /// Set mesh distributor object
    void setMeshDistributor(MeshDistributor *m);

    /// Set mesh distributor object
    void setMeshDistributor(MeshDistributor *md,
86
			    MPI::Intracomm mpiComm0,
87
			    MPI::Intracomm mpiComm1);
88
89
90
91
92
93
94
95
96
97

    /// 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);
98

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

      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
200
201
202
      return vecSol[coarseSpace + 1];
    }

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

210
211
212
213
214
215
216
217
218
219
220
221
222
    /** \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
223
    {
224
      FUNCNAME("ParallelCoarseSpaceSolver::isCoarseSpace()");
225
226
227
228
229
230
231
      
      if (coarseSpaceMap.empty())
	return false;

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

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

235

236
237
238
239
240
241
    /// Returns whether the solver has a coarse grid.
    inline bool hasCoarseSpace() 
    {
      return (!coarseSpaceMap.empty());
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
242
  protected:
243
244
245
246
247
248
249
    /// 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
250
251
    void updateSubdomainData();

252
  private:
253
254
255
256
257
258
259
260
261
262
263
    /// 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.
264
265
    vector<vector<Mat> > mat;

266
267
    /// 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
268
269
    vector<Vec> vecSol, vecRhs;

270
271
    /// Matrix of objects to control the matrix non zero structure of the 
    /// corresponding PETSc matrices stored in \ref mat.
272
273
    vector<vector<MatrixNnzStructure> > nnz;

274
275
276
277
278
279
280
281
282
283
284
285
286
    /** \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.
     */
287
    vector<int> componentIthCoarseMap;
288

289
    /// Stores a set of all coarse space DOF mapping. All entries are unique.
290
291
    vector<ParallelDofMapping*> uniqueCoarseMap;

292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
    /// 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:
307
308
309
    /// Prefix string for parameters in init file.
    string initFileStr;
    
310
    /// Pointer to a mesh distributor object.
311
312
    MeshDistributor *meshDistributor;

313
314
    /// Level of subdomain/interior discretization. Is used for multi-level
    /// methods only.
315
316
    int subdomainLevel;

317
318
319
    /// 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.
320
321
    MPI::Intracomm mpiCommLocal;
    
322
323
    /// 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.
324
325
    MPI::Intracomm mpiCommGlobal;

326
327
328
329
330
    /// 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
331
332
    int rStartInterior;

333
334
335
    /// 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
336
337
    int nGlobalOverallInterior;

338
339
    /// Parallel DOF mapping for the interior.
    ParallelDofMapping *interiorMap;
340

341
342
343
    /// Parallel DOF mapping of the (optional) coarse space. Allows to define
    /// different coarse spaces for different components.
    map<int, ParallelDofMapping*> coarseSpaceMap;   
344
  };
345
  
346
347
348
}

#endif