Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

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