PETScLUSolver.hpp 9.95 KB
Newer Older
1
/*
2
 * Copyright (c) 2017-2022 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Nils Kohl.
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * 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
22
23
#pragma once

#include <memory>

Dominik Thoennes's avatar
Dominik Thoennes committed
24
#include "hyteg/solvers/Solver.hpp"
25

26
27
28
#include "PETScSparseMatrix.hpp"
#include "PETScVector.hpp"

Dominik Thoennes's avatar
Dominik Thoennes committed
29
#ifdef HYTEG_BUILD_WITH_PETSC
30

31
#if ( PETSC_VERSION_MAJOR > 3 ) || ( PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 9 )
Dominik Thoennes's avatar
Dominik Thoennes committed
32
#define HYTEG_PCFactorSetMatSolverType PCFactorSetMatSolverType
Daniel Drzisga's avatar
Daniel Drzisga committed
33
#else
Dominik Thoennes's avatar
Dominik Thoennes committed
34
#define HYTEG_PCFactorSetMatSolverType PCFactorSetMatSolverPackage
Daniel Drzisga's avatar
Daniel Drzisga committed
35
36
#endif

37
namespace hyteg {
38

39
40
41
42
43
44
enum class PETScDirectSolverType
{
   MUMPS,
   SUPER_LU
};

45
template < class OperatorType >
46
47
class PETScLUSolver : public Solver< OperatorType >
{
48
49
50
51
 public:
   typedef typename OperatorType::srcType FunctionType;

   PETScLUSolver( const std::shared_ptr< PrimitiveStorage >& storage, const uint_t& level )
52
53
   : storage_( storage )
   , allocatedLevel_( level )
54
   , petscCommunicator_( storage->getSplitCommunicatorByPrimitiveDistribution() )
55
   , num( "numerator", storage, level, level )
56
57
58
59
60
   , Amat( "Amat", petscCommunicator_ )
   , AmatUnsymmetric( "AmatUnsymmetric", petscCommunicator_ )
   , AmatTmp( "AmatTmp", petscCommunicator_ )
   , xVec( "xVec", petscCommunicator_ )
   , bVec( "bVec", petscCommunicator_ )
61
#if 0
62
  , inKernel( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ) )
63
#endif
64
65
   , flag_( hyteg::All )
   , verbose_( false )
66
   , manualAssemblyAndFactorization_( false )
67
   , reassembleMatrix_( false )
68
   , assumeSymmetry_( true )
69
   , solverType_( PETScDirectSolverType::MUMPS )
70
71
   {
      num.enumerate( level );
72
      KSPCreate( petscCommunicator_, &ksp );
73
74
75
      KSPSetType( ksp, KSPPREONLY );
      KSPSetFromOptions( ksp );
   }
76

77
   ~PETScLUSolver() { KSPDestroy( &ksp ); }
78

79
80
81
82
83
84
85
86
87
88
89
#if 0
  void setNullSpace( FunctionType & inKernel, const uint_t & level )
  {
    inKernel.createVectorFromFunction( inKernel, *num, level, All );
    VecNormalize(inKernel.get(), NULL);
    MatNullSpace nullspace;
    MatNullSpaceCreate( walberla::MPIManager::instance()->comm(), PETSC_FALSE, 1, &(inKernel.get()), &nullspace );
    MatSetNullSpace( Amat.get(), nullspace );
  }
#endif

90
91
   void setDirectSolverType( PETScDirectSolverType solverType ) { solverType_ = solverType; }

92
93
94
   void setConstantNullSpace()
   {
      MatNullSpace nullspace;
95
      MatNullSpaceCreate( petscCommunicator_, PETSC_TRUE, 0, NULL, &nullspace );
96
97
98
99
100
      MatSetNullSpace( Amat.get(), nullspace );
   }

   void setVerbose( bool verbose ) { verbose_ = verbose; }

101
102
103
   /// \brief If set to true, the symmetry of the operator is exploited by the solver.
   void assumeSymmetry( bool assumeSymmetry ) { assumeSymmetry_ = assumeSymmetry; }

104
105
106
107
   /// \brief If set to true, no assembly and no factorization will happen during the solve() call.
   ///        For successful solution of the system, assembleAndFactorize() has to be called before
   ///        the first solve() and after each modification of the operator.
   void setManualAssemblyAndFactorization( bool manualAssemblyAndFactorization )
108
   {
109
110
      manualAssemblyAndFactorization_ = manualAssemblyAndFactorization;
   }
111

112
113
114
   /// \brief If set to true, the operator is reassembled for every solve / manual assembly call.
   ///        Default is false.
   void reassembleMatrix( bool reassembleMatrix ) { reassembleMatrix_ = reassembleMatrix; }
115

116
   void setMUMPSIcntrl( uint_t key, int value ) { mumpsIcntrl_[key] = value; }
Nils Kohl's avatar
Nils Kohl committed
117
118
   void setMUMPSCntrl( uint_t key, real_t value ) { mumpsCntrl_[key] = value; }

119
120
121
   void assembleAndFactorize( const OperatorType& A )
   {
      storage_->getTimingTree()->start( "Matrix assembly" );
122

123
124
125
      bool matrixAssembledForTheFirstTime;
      if ( reassembleMatrix_ )
      {
126
127
         AmatUnsymmetric.zeroEntries();
         AmatUnsymmetric.createMatrixFromOperator( A, allocatedLevel_, num, All );
128
129
130
131
         matrixAssembledForTheFirstTime = true;
      }
      else
      {
132
         matrixAssembledForTheFirstTime = AmatUnsymmetric.createMatrixFromOperatorOnce( A, allocatedLevel_, num, All );
133
      }
134

135
      storage_->getTimingTree()->stop( "Matrix assembly" );
136
137
138

      if ( matrixAssembledForTheFirstTime )
      {
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
         Amat.zeroEntries();
         Amat.createMatrixFromOperator( A, allocatedLevel_, num, All );
         AmatTmp.zeroEntries();
         AmatTmp.createMatrixFromOperator( A, allocatedLevel_, num, All );

         MatCopy( AmatUnsymmetric.get(), Amat.get(), DIFFERENT_NONZERO_PATTERN );

         if ( assumeSymmetry_ )
         {
            Amat.applyDirichletBCSymmetrically( num, allocatedLevel_ );
         }
         else
         {
            Amat.applyDirichletBC( num, allocatedLevel_ );
         }
154

155
         KSPSetOperators( ksp, Amat.get(), Amat.get() );
156
         KSPGetPC( ksp, &pc );
157
158
159
160
161
162
163
164
165
166

         if ( assumeSymmetry_ )
         {
            PCSetType( pc, PCCHOLESKY );
         }
         else
         {
            PCSetType( pc, PCLU );
         }

167
168
169
170
         MatSolverType petscSolverType;
         switch ( solverType_ )
         {
         case PETScDirectSolverType::MUMPS:
171
#ifdef PETSC_HAVE_MUMPS
172
173
            petscSolverType = MATSOLVERMUMPS;
            break;
Dominik Thoennes's avatar
Dominik Thoennes committed
174
#else
175
176
            WALBERLA_ABORT( "PETSc is not build with MUMPS support." )
#endif
177
178
179
180
181
182
183
         case PETScDirectSolverType::SUPER_LU:
            petscSolverType = MATSOLVERSUPERLU_DIST;
            break;
         default:
            WALBERLA_ABORT( "Invalid PETSc solver type." )
         }
         HYTEG_PCFactorSetMatSolverType( pc, petscSolverType );
Nils Kohl's avatar
Nils Kohl committed
184
185
186
187
188

         if ( solverType_ == PETScDirectSolverType::MUMPS )
         {
            PCFactorSetUpMatSolverType( pc );
            PCFactorGetMatrix( pc, &F );
189
#ifdef PETSC_HAVE_MUMPS
Nils Kohl's avatar
Nils Kohl committed
190
191
192
193
194
195
196
197
            for ( auto it : mumpsIcntrl_ )
            {
               MatMumpsSetIcntl( F, it.first, it.second );
            }
            for ( auto it : mumpsCntrl_ )
            {
               MatMumpsSetCntl( F, it.first, it.second );
            }
198
#endif
Nils Kohl's avatar
Nils Kohl committed
199
         }
200
201
202
         storage_->getTimingTree()->start( "Factorization" );
         PCSetUp( pc );
         storage_->getTimingTree()->stop( "Factorization" );
203
      }
204
205
206
207
208
209
210
211
212
213
214
   }

   void solve( const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level )
   {
      WALBERLA_CHECK_EQUAL( level, allocatedLevel_ );

      walberla::WcTimer timer;

      storage_->getTimingTree()->start( "PETSc LU Solver" );
      storage_->getTimingTree()->start( "Setup" );

215
      timer.start();
216
217
218
219
      if ( !manualAssemblyAndFactorization_ )
      {
         assembleAndFactorize( A );
      }
220
221
      timer.end();
      const double matrixAssemblyAndFactorizationTime = timer.last();
222

223
224
      storage_->getTimingTree()->start( "RHS vector setup" );

225
      b.assign( { 1.0 }, { x }, level, DirichletBoundary );
226
      bVec.createVectorFromFunction( b, num, level, All );
Nils Kohl's avatar
Nils Kohl committed
227
      xVec.createVectorFromFunction( x, num, level, All );
228
229
230
231
232
233
234
235
236
237

      if ( assumeSymmetry_ )
      {
         AmatTmp.zeroEntries();
         MatCopy( AmatUnsymmetric.get(), AmatTmp.get(), DIFFERENT_NONZERO_PATTERN );
         AmatTmp.applyDirichletBCSymmetrically( x, num, bVec, allocatedLevel_ );
      }

      storage_->getTimingTree()->stop( "RHS vector setup" );

238
      storage_->getTimingTree()->stop( "Setup" );
239

240
      storage_->getTimingTree()->start( "Solver" );
241
242
243
244
      timer.start();
      KSPSolve( ksp, bVec.get(), xVec.get() );
      timer.end();
      const double petscKSPTimer = timer.last();
245
      storage_->getTimingTree()->stop( "Solver" );
246
247
248
249
250
251

      xVec.createFunctionFromVector( x, num, level, flag_ );

      if ( verbose_ )
      {
         WALBERLA_LOG_INFO_ON_ROOT( "[PETScLUSolver] "
252
253
                                    << "PETSc KSPSolver time: " << petscKSPTimer
                                    << ", assembly and fact time: " << matrixAssemblyAndFactorizationTime );
254
255
      }

256
      storage_->getTimingTree()->stop( "PETSc LU Solver" );
257
258
259
   }

 private:
260
   std::shared_ptr< PrimitiveStorage >                                                           storage_;
261
   uint_t                                                                                        allocatedLevel_;
Nils Kohl's avatar
Nils Kohl committed
262
   MPI_Comm                                                                                      petscCommunicator_;
263
   typename OperatorType::srcType::template FunctionType< idx_t >                                num;
264
265
266
   PETScSparseMatrix< OperatorType >                                                             Amat;
   PETScSparseMatrix< OperatorType >                                                             AmatUnsymmetric;
   PETScSparseMatrix< OperatorType >                                                             AmatTmp;
267
   PETScVector< typename FunctionType::valueType, OperatorType::srcType::template FunctionType > xVec;
268
   PETScVector< typename FunctionType::valueType, OperatorType::dstType::template FunctionType > bVec;
269
270
271
272
#if 0
  PETScVector<typename FunctionType::valueType, OperatorType::srcType::template FunctionType> inKernel;
#endif

Nils Kohl's avatar
Nils Kohl committed
273
274
275
276
277
278
279
280
281
   KSP                        ksp;
   PC                         pc;
   hyteg::DoFType             flag_;
   bool                       verbose_;
   Mat                        F; //factored Matrix
   bool                       manualAssemblyAndFactorization_;
   bool                       reassembleMatrix_;
   bool                       assumeSymmetry_;
   PETScDirectSolverType      solverType_;
282
   std::map< uint_t, int >    mumpsIcntrl_;
Nils Kohl's avatar
Nils Kohl committed
283
   std::map< uint_t, real_t > mumpsCntrl_;
284
285
};

286
} // namespace hyteg
287

288
#endif