dune-fem 2.8-git
petscinverseoperators.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_PETSCINVERSEOPERATORS_HH
2#define DUNE_FEM_PETSCINVERSEOPERATORS_HH
3
4#include <limits>
5
10
11#if HAVE_PETSC
16
17namespace Dune
18{
19
20 namespace Fem
21 {
22
23 //=====================================================================
24 // Implementation for PETSc matrix based Krylov solvers
25 //=====================================================================
26
31 // PETScSolver
32 // --------------
33 template< class DF, class Op = Dune::Fem::Operator< DF, DF > >
34 class PetscInverseOperator;
35
36 template <class DF, class Op >
37 struct PetscInverseOperatorTraits
38 {
39 private:
40 typedef typename DF :: DiscreteFunctionSpaceType SpaceType ;
41 public:
42 typedef DF DiscreteFunctionType;
43 typedef Op OperatorType;
44 typedef OperatorType PreconditionerType;
45 typedef PetscDiscreteFunction< SpaceType > SolverDiscreteFunctionType;
46 typedef PetscLinearOperator< DF, DF > AssembledOperatorType;
47 typedef PetscInverseOperator< DF, Op > InverseOperatorType;
48 typedef PetscSolverParameter SolverParameterType;
49 };
50
51
53 template< class DF, class Op >
54 class PetscInverseOperator
55 : public InverseOperatorInterface< PetscInverseOperatorTraits< DF, Op > >
56 {
57 protected:
58 // monitor function for PETSc solvers
59 static PetscErrorCode
60 monitor (KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
61 {
63 {
64 std::cout << "PETSc::KSP: it = " << it << " res = " << rnorm << std::endl;
65 }
66 return PetscErrorCode(0);
67 }
68
69 // destroy solver context
70 struct KSPDeleter
71 {
72 void operator() ( KSP* p ) const
73 {
74 if( !p )
75 return;
76
77 ::Dune::Petsc::KSPDestroy( p );
78 delete p;
79 }
80 };
81
82 typedef PetscInverseOperatorTraits< DF, Op > Traits;
83 typedef InverseOperatorInterface< Traits > BaseType;
84 friend class InverseOperatorInterface< Traits >;
85 public:
86
87 typedef typename BaseType :: SolverDiscreteFunctionType PetscDiscreteFunctionType;
88 typedef typename BaseType :: OperatorType OperatorType;
89
90 PetscInverseOperator ( const PetscSolverParameter &parameter = PetscSolverParameter(Parameter::container()) )
91 : BaseType( parameter )
92 {}
93
94 PetscInverseOperator ( const OperatorType &op, const PetscSolverParameter &parameter = PetscSolverParameter(Parameter::container()) )
95 : BaseType( parameter )
96 {
97 bind( op );
98 }
99
100 void bind ( const OperatorType &op )
101 {
102 BaseType :: bind( op );
103 initialize( *parameter_ );
104 }
105
106 void unbind ()
107 {
108 BaseType :: unbind();
109 ksp_.reset();
110 }
111
112 void printTexInfo(std::ostream& out) const
113 {
114 out << "Solver: " << solverName_ << " eps = " << parameter_->tolerance() ;
115 out << "\\\\ \n";
116 }
117
118 protected:
119 void initialize ( const PetscSolverParameter& parameter )
120 {
121 if( !assembledOperator_ )
122 DUNE_THROW(NotImplemented,"Petsc solver with matrix free implementations not yet supported!");
123
124 // Create linear solver context
125 ksp_.reset( new KSP() );
126 const auto& comm = assembledOperator_->domainSpace().gridPart().comm();
127 ::Dune::Petsc::KSPCreate( comm, &ksp() );
128
129 // attach Matrix to linear solver context
130 Mat& A = assembledOperator_->exportMatrix();
131#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
132 ::Dune::Petsc::KSPSetOperators( ksp(), A, A, SAME_PRECONDITIONER);
133#else
134 ::Dune::Petsc::KSPSetOperators( ksp(), A, A );
135#endif
136
137 // allow for non-zero initial guess
138 ::Dune::Petsc::KSPSetInitialGuessNonzero( ksp(), PETSC_TRUE );
139
140 // set prescribed tolerances
141 PetscInt maxits = parameter_->maxIterations();
142 PetscReal tolerance = parameter_->tolerance();
143 if (parameter_->errorMeasure() == 0)
144 ::Dune::Petsc::KSPSetTolerances(ksp(), 1e-50, tolerance, PETSC_DEFAULT, maxits);
145 else
146 ::Dune::Petsc::KSPSetTolerances(ksp(), tolerance, 1e-50, PETSC_DEFAULT, maxits);
147
148 enum class PetscSolver {
154 preonly = SolverParameter::preonly,
155 kspoptions = 0
156 };
157
158 // if special petsc solver parameter exists use that one, otherwise
159 // use solverMethod from SolverParameter
160 const auto& reader = parameter.parameter();
161 PetscSolver kspType = PetscSolver::gmres;
162 if( reader.exists("petsc.kspsolver.method") )
163 {
164 // see PETSc docu for more types
165 const std::string kspNames[] = { "default", "cg", "bicgstab", "gmres", "minres", "gradient", "loop", "superlu", "bicg", "preonly" };
166 kspType = static_cast< PetscSolver >( reader.getEnum("petsc.kspsolver.method", kspNames, int(PetscSolver::gmres) ) );
167 std::cout << "WARNING: using deprecated parameter 'petsc.kpsolver.method' use "
168 << parameter.keyPrefix() << "method instead\n";
169 }
170 else
171 kspType = static_cast< PetscSolver >(
172 parameter.solverMethod (
179 }, { "kspoptions" } )
180 );
181
182 if (kspType > PetscSolver::kspoptions)
183 solverName_ = SolverParameter::solverMethodTable( static_cast< int >( kspType ) );
184 else
185 solverName_ = "kspoptions";
186
187 // select linear solver
188 switch( kspType )
189 {
190 case PetscSolver::cg:
191 ::Dune::Petsc::KSPSetType( ksp(), KSPCG );
192 break;
194 ::Dune::Petsc::KSPSetType( ksp(), KSPBCGS );
195 break;
197 {
198 ::Dune::Petsc::KSPSetType( ksp(), KSPGMRES );
199 PetscInt restart = 10;
200 if( reader.exists("petsc.gmresrestart") )
201 {
202 restart = reader.getValue<int>("petsc.gmresrestart", restart );
203 std::cout << "WARNING: using deprecated parameter 'petsc.gmresrestart' use "
204 << parameter.keyPrefix() << "gmres.restart instead\n";
205 }
206 else
207 restart = parameter.gmresRestart() ;
208
209 ::Dune::Petsc::KSPGMRESSetRestart( ksp(), restart );
210 break;
211 }
212 case PetscSolver::minres:
213 ::Dune::Petsc::KSPSetType( ksp(), KSPMINRES );
214 break;
215 case PetscSolver::bicg:
216 ::Dune::Petsc::KSPSetType( ksp(), KSPBICG );
217 break;
218 case PetscSolver::preonly:
219 ::Dune::Petsc::KSPSetType( ksp(), KSPPREONLY );
220 break;
221 case PetscSolver::kspoptions:
222 // setup solver context from database/cmdline options
223 ::Dune::Petsc::KSPSetFromOptions( ksp() );
224 ::Dune::Petsc::KSPSetUp( ksp() );
225 break;
226 default:
227 DUNE_THROW(InvalidStateException,"PetscInverseOperator: invalid solver choosen." );
228 }
229
231 // preconditioning
233
234 int pcType = SolverParameter::none;
235 if( reader.exists("petsc.preconditioning.method") )
236 {
237 const std::string pcNames[] = { "default", "none", "asm", "sor", "jacobi", "ilu", "icc", "superlu", "hypre", "ml", "lu" };
238 pcType = reader.getEnum("petsc.preconditioning.method", pcNames, 0 );
239 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.method' use "
240 << parameter.keyPrefix() << "preconditioning.method instead\n";
241 if (pcType >= 8)
242 pcType = 7-pcType; // hypre=-1, ml=-2, lu=-3
243 }
244 else
245 {
246 pcType = parameter.preconditionMethod(
247 {
248 SolverParameter::none, // no preconditioning
249 SolverParameter::oas, // Overlapping Additive Schwarz
250 SolverParameter::sor, // SOR and SSOR
251 SolverParameter::jacobi, // Jacobi preconditioning
252 SolverParameter::ilu, // ILU preconditioning
253 SolverParameter::icc, // Incomplete Cholesky factorization
254 SolverParameter::superlu // SuperLU direct factorization
255 },
256 {"kspoptions", // = 0, // use command line options -ksp...
257 "hypre", // = -1, // Hypre preconditioning
258 "ml", // = -2, // ML preconditioner (from Trilinos)
259 "lu", // = -3, // LU factorization
260 "pcgamg", // = -4 // Petsc internal AMG
261 });
262 }
263
264 // setup preconditioning context
265 PC pc;
266 ::Dune::Petsc::KSPGetPC( ksp(), &pc );
267
268 switch( pcType )
269 {
270 case 0:
271 // don't setup the pc context twice
272 if ( kspType != PetscSolver::kspoptions )
273 {
274 // setup pc context from database/cmdline options
275 ::Dune::Petsc::PCSetFromOptions( pc );
276 ::Dune::Petsc::PCSetUp( pc );
277 }
278 break;
280 ::Dune::Petsc::PCSetType( pc, PCNONE );
281 break;
283 {
284 ::Dune::Petsc::PCSetType( pc, PCASM );
285 ::Dune::Petsc::PCSetUp( pc );
286 break;
287 }
289 ::Dune::Petsc::PCSetType( pc, PCSOR );
290 break;
292 ::Dune::Petsc::PCSetType( pc, PCJACOBI );
293 break;
294 case -1: // PetscPrec::hypre:
295 {
296#ifdef PETSC_HAVE_HYPRE
297 int hypreType = parameter.hypreMethod();
298 std::string hypre;
299 if ( hypreType == PetscSolverParameter::boomeramg )
300 hypre = "boomeramg";
301 else if ( hypreType == PetscSolverParameter::parasails )
302 hypre = "parasails";
303 else if ( hypreType == PetscSolverParameter::pilut )
304 hypre = "pilut";
305 else
306 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid hypre preconditioner choosen." );
307
308 ::Dune::Petsc::PCSetType( pc, PCHYPRE );
309 ::Dune::Petsc::PCHYPRESetType( pc, hypre.c_str() );
310 ::Dune::Petsc::PCSetUp( pc );
311#else // PETSC_HAVE_HYPRE
312 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with hypre support." );
313#endif // PETSC_HAVE_HYPRE
314 break;
315 }
316 case -2: // PetscPrec::ml:
317#ifdef PETSC_HAVE_ML
318 ::Dune::Petsc::PCSetType( pc, PCML );
319#else // PETSC_HAVE_ML
320 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with ml support." );
321#endif // PETSC_HAVE_ML
322 break;
324 {
325 if ( MPIManager::size() > 1 )
326 DUNE_THROW( InvalidStateException, "PetscInverseOperator: ilu preconditioner does not work in parallel." );
327
328 // get fill-in level
329 PetscInt pcLevel;
330 if( reader.exists("petsc.preconditioning.levels") )
331 {
332 pcLevel = reader.getValue<int>("petsc.preconditioning.levels", 0 );
333 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.levels' use "
334 << parameter.keyPrefix() << "preconditioning.level instead\n";
335 }
336 else
337 pcLevel = parameter.preconditionerLevel() ;
338
339 ::Dune::Petsc::PCSetType( pc, PCILU );
340 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
341 break;
342 }
343 ::Dune::Petsc::PCSetType( pc, PCML );
344 break;
346 {
347#ifdef PETSC_HAVE_ICC
348 if ( MPIManager::size() > 1 )
349 DUNE_THROW( InvalidStateException, "PetscInverseOperator: icc preconditioner does not worl in parallel." );
350
351 // get fill-in level
352 PetscInt pcLevel;
353 if( reader.exists("petsc.preconditioning.levels") )
354 {
355 pcLevel = reader.getValue<int>("petsc.preconditioning.levels", 0 );
356 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.levels' use "
357 << parameter.keyPrefix() << "preconditioning.level instead\n";
358 }
359 else
360 pcLevel = parameter.preconditionerLevel() ;
361
362
363 ::Dune::Petsc::PCSetType( pc, PCICC );
364 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
365#else // PETSC_HAVE_ICC
366 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with icc support." );
367#endif // PETSC_HAVE_ICC
368 break;
369 }
370 case -3: // PetscPrec::lu:
372 {
373 enum class Factorization { petsc = 0, superlu = 1, mumps = 2 };
374 Factorization factorType = Factorization::superlu;
375 if (pcType != SolverParameter::superlu)
376 factorType = static_cast<Factorization>(parameter.superluMethod());
377
378 ::Dune::Petsc::PCSetType( pc, PCLU );
379
380 if ( factorType == Factorization::petsc )
381 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERPETSC );
382 else if ( factorType == Factorization::superlu )
383 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERSUPERLU_DIST );
384 else if ( factorType == Factorization::mumps )
385 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERMUMPS );
386 else
387 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid factorization package choosen." );
388
389 ::Dune::Petsc::PCSetUp( pc );
390 break;
391 }
392 case -4: // PetscPrec::pcgamg:
393 ::Dune::Petsc::PCSetType( pc, PCGAMG );
394 break;
395
396 default:
397 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid preconditioner choosen." );
398 }
399
400 // set monitor in verbose mode
401 if( parameter.verbose() && Parameter::verbose() )
402 {
403 ::Dune::Petsc::KSPView( comm, ksp() );
404 ::Dune::Petsc::KSPMonitorSet( ksp(), &monitor, PETSC_NULL, PETSC_NULL);
405 }
406 }
407
408 int apply( const PetscDiscreteFunctionType& arg, PetscDiscreteFunctionType& dest ) const
409 {
410 // need to have a 'distributed' destination vector for continuous spaces
411 if( dest.space().continuous() )
412 dest.dofVector().clearGhost();
413
414 // call PETSc solvers
415 ::Dune::Petsc::KSPSolve( *ksp_, *arg.petscVec() , *dest.petscVec() );
416
417 // a continuous solution is 'distributed' so need a communication here
418 if( dest.space().continuous() )
419 {
420 dest.communicate();
421 }
422
423 // get number of iterations
424 PetscInt its ;
425 ::Dune::Petsc::KSPGetIterationNumber( *ksp_, &its );
426 KSPConvergedReason reason;
427 ::Dune::Petsc::KSPGetConvergedReason( *ksp_, &reason );
428 if( parameter_->verbose() && Parameter::verbose() )
429 {
430 // list of reasons:
431 // https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPConvergedReason.html
432 std::cout << "Converged reason:" << reason << std::endl;
433 }
434
435 bool converged = int(reason) >= 0 ;
436
437 return (converged) ? its : -its;
438 }
439
440 protected:
441 KSP & ksp () { assert( ksp_ ); return *ksp_; }
442
443 using BaseType :: assembledOperator_;
444 using BaseType :: parameter_;
445 using BaseType :: iterations_;
446
447 std::unique_ptr< KSP, KSPDeleter > ksp_; // PETSc Krylov Space solver context
448
449 std::string solverName_;
450 };
451
453
454 } // namespace Fem
455
456} // namespace Dune
457
458#endif // #if HAVE_PETSC
459
460#endif // #ifndef DUNE_FEM_PETSCINVERSEOPERATORS_HH
Definition: bindguard.hh:11
int gmres(Operator &op, Preconditioner *preconditioner, std::vector< DiscreteFunction > &v, DiscreteFunction &u, const DiscreteFunction &b, const int m, const double tolerance, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: gmres.hh:120
int cg(Operator &op, Precoditioner *preconditioner, std::vector< DiscreteFunction > &tempMem, DiscreteFunction &x, const DiscreteFunction &b, const double epsilon, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: cg.hh:24
int bicgstab(Operator &op, Precoditioner *preconditioner, std::vector< DiscreteFunction > &tempMem, DiscreteFunction &x, const DiscreteFunction &b, const double tolerance, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: bicgstab.hh:64
static ParameterContainer & container()
Definition: io/parameter.hh:193
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static int size()
Definition: mpimanager.hh:160
static const std::string solverMethodTable(int i)
Definition: solver/parameter.hh:33
static const int none
Definition: solver/parameter.hh:40
static const int cg
Definition: solver/parameter.hh:24
static const int icc
Definition: solver/parameter.hh:50
static const int gmres
Definition: solver/parameter.hh:26
static const int bicg
Definition: solver/parameter.hh:31
static const int sor
Definition: solver/parameter.hh:42
static const int minres
Definition: solver/parameter.hh:27
static const int ilu
Definition: solver/parameter.hh:43
static const int bicgstab
Definition: solver/parameter.hh:25
static const int superlu
Definition: solver/parameter.hh:30
static const int jacobi
Definition: solver/parameter.hh:45
static const int preonly
Definition: solver/parameter.hh:32
static const int oas
Definition: solver/parameter.hh:49