PETScSparseMatrix.hpp 11.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*
 * Copyright (c) 2017-2019 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Marcus Mohr, Nils Kohl.
 *
 * This file is part of HyTeG
 * (see https://i10git.cs.fau.de/hyteg/hyteg).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
20
21
#pragma once

22
#include "hyteg/types/types.hpp"
23

Nils Kohl's avatar
Nils Kohl committed
24
#include "PETScWrapper.hpp"
25

Nils Kohl's avatar
Nils Kohl committed
26
#ifdef HYTEG_BUILD_WITH_PETSC
27

Nils Kohl's avatar
Nils Kohl committed
28
#include "hyteg/composites/UnsteadyDiffusion.hpp"
Dominik Thoennes's avatar
Dominik Thoennes committed
29
#include "hyteg/composites/petsc/P1StokesPetsc.hpp"
30
#include "hyteg/composites/StrongFreeSlipWrapper.hpp"
Dominik Thoennes's avatar
Dominik Thoennes committed
31
#include "hyteg/composites/petsc/P2P1TaylorHoodPetsc.hpp"
32
#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
33
#include "hyteg/elementwiseoperators/DiagonalNonConstantOperator.hpp"
34
#include "hyteg/elementwiseoperators/ElementwiseOperatorPetsc.hpp"
Nils Kohl's avatar
Nils Kohl committed
35
36
#include "hyteg/p1functionspace/P1Petsc.hpp"
#include "hyteg/p2functionspace/P2Petsc.hpp"
37
#include "hyteg/petsc/PETScSparseMatrixInfo.hpp"
38
#include "hyteg/petsc/PETScSparseMatrixProxy.hpp"
39
#include "hyteg/petsc/PETScVector.hpp"
40

41
namespace hyteg {
42

43
/// Wrapper class for PETSc sparse matrix usage
Nils Kohl's avatar
Nils Kohl committed
44
45
46
47
48
49
template < class OperatorType, template < class > class FunctionType >
class PETScSparseMatrix
{
 public:
   PETScSparseMatrix() = delete;

50
51
52
53
   PETScSparseMatrix( uint_t          localSize,
                      uint_t          globalSize,
                      const char      name[]            = "Mat",
                      const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
54
   : petscCommunicator_( petscCommunicator ), assembled_( false )
Nils Kohl's avatar
Nils Kohl committed
55
   {
56
      MatCreate( petscCommunicator, &mat );
Nils Kohl's avatar
Nils Kohl committed
57
      MatSetType( mat, MATMPIAIJ );
58
59
60
61
62
      MatSetSizes( mat,
                   static_cast< PetscInt >( localSize ),
                   static_cast< PetscInt >( localSize ),
                   static_cast< PetscInt >( globalSize ),
                   static_cast< PetscInt >( globalSize ) );
Nils Kohl's avatar
Nils Kohl committed
63
64
65
66
67
68
      // Roughly overestimate number of non-zero entries for faster assembly of matrix
      MatMPIAIJSetPreallocation( mat, 500, NULL, 500, NULL );
      setName( name );
      reset();
   }

69
70
71
72
73
74
75
76
77
78
   PETScSparseMatrix( const std::shared_ptr< PrimitiveStorage >& storage,
                      const uint_t&                              level,
                      const char                                 name[] = "Mat",
                      const MPI_Comm& petscCommunicator                 = walberla::mpi::MPIManager::instance()->comm() )
   : PETScSparseMatrix( numberOfLocalDoFs< typename OperatorType::dstType::Tag >( *storage, level ),
                        numberOfGlobalDoFs< typename OperatorType::dstType::Tag >( *storage, level, petscCommunicator ),
                        name,
                        petscCommunicator )
   {}

Nils Kohl's avatar
Nils Kohl committed
79
80
   virtual ~PETScSparseMatrix() { MatDestroy( &mat ); }

81
82
83
84
   inline void createMatrixFromOperator( const OperatorType&          op,
                                         uint_t                       level,
                                         const FunctionType< idx_t >& numerator,
                                         DoFType                      flag = All )
Nils Kohl's avatar
Nils Kohl committed
85
   {
86
87
      auto proxy = std::make_shared< PETScSparseMatrixProxy >( mat );
      hyteg::petsc::createMatrix< OperatorType >( op, numerator, numerator, proxy, level, flag );
Nils Kohl's avatar
Nils Kohl committed
88
89
90

      MatAssemblyBegin( mat, MAT_FINAL_ASSEMBLY );
      MatAssemblyEnd( mat, MAT_FINAL_ASSEMBLY );
91
      assembled_ = true;
Nils Kohl's avatar
Nils Kohl committed
92
93
   }

94
95
96
97
   inline bool createMatrixFromOperatorOnce( const OperatorType&          op,
                                             uint_t                       level,
                                             const FunctionType< idx_t >& numerator,
                                             DoFType                      flag = All )
Nils Kohl's avatar
Nils Kohl committed
98
   {
99
      if ( assembled_ )
Nils Kohl's avatar
Nils Kohl committed
100
101
102
103
104
         return false;
      createMatrixFromOperator( op, level, numerator, flag );
      return true;
   }

105
   inline void print( const std::string& name, bool binary = false, PetscViewerFormat format = PETSC_VIEWER_ASCII_MATRIXMARKET )
Nils Kohl's avatar
Nils Kohl committed
106
107
   {
      PetscViewer viewer;
108
109
110
111
112
113
114
115
116
      if ( binary )
      {
         PetscViewerBinaryOpen( petscCommunicator_, name.c_str(), FILE_MODE_WRITE, &viewer );
      }
      else
      {
         PetscViewerASCIIOpen( petscCommunicator_, name.c_str(), &viewer );
         PetscViewerPushFormat( viewer, format );
      }
Nils Kohl's avatar
Nils Kohl committed
117
118
119
120
      MatView( mat, viewer );
      PetscViewerDestroy( &viewer );
   }

121
122
   template < typename PETSCINT >
   std::vector< PETSCINT > convertToPetscVector( std::vector< idx_t > idx_vector )
Nils Kohl's avatar
Nils Kohl committed
123
   {
124
      if constexpr ( std::is_same_v< idx_t, PETSCINT > )
125
      {
126
         return idx_vector;
127
128
129
      }
      else
      {
130
         return std::vector< PETSCINT >( idx_vector.begin(), idx_vector.end() );
131
      }
132
133
134
135
136
137
138
   }

   void applyDirichletBC( const FunctionType< idx_t >& numerator, uint_t level )
   {
      std::vector< idx_t > bcIndices;
      hyteg::petsc::applyDirichletBC( numerator, bcIndices, level );
      std::vector< PetscInt > PetscIntBcIndices = convertToPetscVector< PetscInt >( bcIndices );
Nils Kohl's avatar
Nils Kohl committed
139

140
141
142
143
144
      // This is required as the implementation of MatZeroRows() checks (for performance reasons?!)
      // if there are zero diagonals in the matrix. If there are, the function halts.
      // To disable that check, we need to allow setting MAT_NEW_NONZERO_LOCATIONS to true.
      MatSetOption( mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE );

145
      MatZeroRows( mat, static_cast< PetscInt >( PetscIntBcIndices.size() ), PetscIntBcIndices.data(), 1.0, nullptr, nullptr );
Nils Kohl's avatar
Nils Kohl committed
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

      MatAssemblyBegin( mat, MAT_FINAL_ASSEMBLY );
      MatAssemblyEnd( mat, MAT_FINAL_ASSEMBLY );
   }

   /// \brief Applies Dirichlet BCs to a linear system without losing symmetry.
   ///
   /// Uses the PETSc function MatZeroRowsColumns() which does that automatically.
   /// Still, we need to think how we can easily integrate this to use more efficient
   /// solvers in HyTeG, because the RHS is modified depending on the original system.
   ///
   /// So far I do not know any solution to this without re-assembling the system every time
   /// we solve it since we need to also rebuild the RHS.
   /// It should be possible to store a copy of the original system and circumvent re-assembling by
   /// copying it and applying only MatZeroRowsColumns() (without re-assembly) before calling the solver.
   /// If PETSc is only used as a coarse grid solver this might be a good solution.
   ///
   /// \param dirichletSolution a function that has the respective values interpolated on the Dirichlet boundary
   /// \param numerator an enumerated function
   /// \param rhsVec RHS of the system as PETSc vector - NOTE THAT THIS IS MODIFIED IN PLACE
   /// \param level the refinement level
   ///
   void applyDirichletBCSymmetrically( const FunctionType< real_t >&        dirichletSolution,
169
                                       const FunctionType< idx_t >&         numerator,
Nils Kohl's avatar
Nils Kohl committed
170
171
172
                                       PETScVector< real_t, FunctionType >& rhsVec,
                                       const uint_t&                        level )
   {
173
174
175
      std::vector< idx_t > bcIndices;
      hyteg::petsc::applyDirichletBC( numerator, bcIndices, level );
      std::vector< PetscInt > PetscIntBcIndices = convertToPetscVector< PetscInt >( bcIndices );
Nils Kohl's avatar
Nils Kohl committed
176

177
178
      PETScVector< real_t, FunctionType > dirichletSolutionVec(
          dirichletSolution, numerator, level, All, "dirichletSolutionVec", rhsVec.getCommunicator() );
Nils Kohl's avatar
Nils Kohl committed
179
180
181
182
183
184

      // This is required as the implementation of MatZeroRowsColumns() checks (for performance reasons?!)
      // if there are zero diagonals in the matrix. If there are, the function halts.
      // To disable that check, we need to allow setting MAT_NEW_NONZERO_LOCATIONS to true.
      MatSetOption( mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE );

185
      MatZeroRowsColumns(
186
          mat, PetscIntBcIndices.size(), PetscIntBcIndices.data(), 1.0, dirichletSolutionVec.get(), rhsVec.get() );
Nils Kohl's avatar
Nils Kohl committed
187
188
189
190
191
   }

   /// \brief Variant of applyDirichletBCSymmetrically() that only modifies the matrix itself
   ///
   /// \return Vector with global indices of the Dirichlet DoFs
192
   std::vector< idx_t > applyDirichletBCSymmetrically( const FunctionType< idx_t >& numerator, const uint_t& level )
Nils Kohl's avatar
Nils Kohl committed
193
   {
194
195
196
      std::vector< idx_t > bcIndices;
      hyteg::petsc::applyDirichletBC( numerator, bcIndices, level );
      std::vector< PetscInt > PetscIntBcIndices = convertToPetscVector< PetscInt >( bcIndices );
Nils Kohl's avatar
Nils Kohl committed
197
198
199
200
201
202

      // This is required as the implementation of MatZeroRowsColumns() checks (for performance reasons?!)
      // if there are zero diagonals in the matrix. If there are, the function halts.
      // To disable that check, we need to allow setting MAT_NEW_NONZERO_LOCATIONS to true.
      MatSetOption( mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE );

203
      MatZeroRowsColumns(
204
          mat, static_cast< PetscInt >( PetscIntBcIndices.size() ), PetscIntBcIndices.data(), 1.0, nullptr, nullptr );
Nils Kohl's avatar
Nils Kohl committed
205

206
      return bcIndices;
Nils Kohl's avatar
Nils Kohl committed
207
208
   }

209
   inline void reset() { assembled_ = false; }
Nils Kohl's avatar
Nils Kohl committed
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234

   /// \brief Sets all entries of the matrix to zero.
   inline void zeroEntries() { MatZeroEntries( mat ); }

   inline void setName( const char name[] ) { PetscObjectSetName( (PetscObject) mat, name ); }

   inline Mat& get() { return mat; }

   bool isSymmetric( real_t tol = real_c( 1e-13 ) )
   {
      Mat       B;
      PetscReal norm;
      MatTranspose( mat, MAT_INITIAL_MATRIX, &B );
      MatAYPX( B, -1.0, mat, DIFFERENT_NONZERO_PATTERN );
      MatNorm( B, NORM_INFINITY, &norm );
      // WALBERLA_LOG_DEVEL_ON_ROOT( "PETSC_NORM = " << norm );
      MatDestroy( &B );
      return norm < tol;
   }

   bool isDiagonal( real_t tol = real_c( 1e-13 ) )
   {
      Mat       B;
      PetscReal norm;
      Vec       diag;
235
      MatCreate( petscCommunicator_, &B );
Nils Kohl's avatar
Nils Kohl committed
236
237
238
239
240
241
242
      MatSetType( B, MATMPIAIJ );
      PetscInt localSize, globalSize;
      MatGetSize( mat, &localSize, &globalSize );
      MatSetSizes( B, localSize, localSize, globalSize, globalSize );
      MatSetUp( B );
      MatAssemblyBegin( B, MAT_FINAL_ASSEMBLY );
      MatAssemblyEnd( B, MAT_FINAL_ASSEMBLY );
243
      VecCreate( petscCommunicator_, &diag );
Nils Kohl's avatar
Nils Kohl committed
244
245
246
247
248
249
250
251
252
253
254
255
256
      VecSetType( diag, VECMPI );
      VecSetSizes( diag, localSize, globalSize );
      VecSetUp( diag );
      MatCopy( mat, B, DIFFERENT_NONZERO_PATTERN );
      MatGetDiagonal( B, diag );
      VecScale( diag, -1.0 );
      MatDiagonalSet( B, diag, ADD_VALUES );
      MatNorm( B, NORM_INFINITY, &norm );
      // WALBERLA_LOG_DEVEL_ON_ROOT( "PETSC_NORM = " << norm );
      MatDestroy( &B );
      VecDestroy( &diag );
      return norm < tol;
   }
257

258
259
260
261
262
263
264
265
266
267
268
269
   PETScSparseMatrixInfo getInfo()
   {
      if ( !assembled_ )
      {
         WALBERLA_ABORT( "Matrix assembly must be complete before calling getInfo()!" );
      }
      MatInfo info;
      MatGetInfo( mat, MAT_GLOBAL_SUM, &info );
      PETScSparseMatrixInfo matInfo( info );
      return matInfo;
   };

270
271
272
 protected:
   MPI_Comm petscCommunicator_;
   Mat      mat;
273
   bool     assembled_;
274
275
};

Nils Kohl's avatar
Nils Kohl committed
276
} // namespace hyteg
Daniel Drzisga's avatar
Daniel Drzisga committed
277
278

#endif