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 
17 namespace 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  {
62  if( Parameter :: verbose () )
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 {
152  minres = SolverParameter::minres,
153  bicg = SolverParameter::bicg,
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;
196  case PetscSolver::gmres:
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 bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static ParameterContainer & container()
Definition: io/parameter.hh:193
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