dune-fem 2.8.0
Loading...
Searching...
No Matches
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 {
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
128 ::Dune::Petsc::KSPCreate( comm, &ksp() );
129
130 // attach Matrix to linear solver context
131 Mat& A = assembledOperator_->exportMatrix();
132#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
133 ::Dune::Petsc::KSPSetOperators( ksp(), A, A, SAME_PRECONDITIONER);
134#else
135 ::Dune::Petsc::KSPSetOperators( ksp(), A, A );
136#endif
137
138 // allow for non-zero initial guess
139 ::Dune::Petsc::KSPSetInitialGuessNonzero( ksp(), PETSC_TRUE );
140
141 // set prescribed tolerances
142 PetscInt maxits = parameter_->maxIterations();
143 PetscReal tolerance = parameter_->tolerance();
144 if (parameter_->errorMeasure() == 0)
145 ::Dune::Petsc::KSPSetTolerances(ksp(), 1e-50, tolerance, PETSC_DEFAULT, maxits);
146 else
147 ::Dune::Petsc::KSPSetTolerances(ksp(), tolerance, 1e-50, PETSC_DEFAULT, maxits);
148
149 enum class PetscSolver {
150 cg = SolverParameter::cg,
151 bicgstab = SolverParameter::bicgstab,
152 gmres = SolverParameter::gmres,
153 minres = SolverParameter::minres,
154 bicg = SolverParameter::bicg,
155 preonly = SolverParameter::preonly,
156 kspoptions = 0
157 };
158
159 // if special petsc solver parameter exists use that one, otherwise
160 // use solverMethod from SolverParameter
161 const auto& reader = parameter.parameter();
162 PetscSolver kspType = PetscSolver::gmres;
163 if( reader.exists("petsc.kspsolver.method") )
164 {
165 // see PETSc docu for more types
166 const std::string kspNames[] = { "default", "cg", "bicgstab", "gmres", "minres", "gradient", "loop", "superlu", "bicg", "preonly" };
167 kspType = static_cast< PetscSolver >( reader.getEnum("petsc.kspsolver.method", kspNames, int(PetscSolver::gmres) ) );
168 std::cout << "WARNING: using deprecated parameter 'petsc.kpsolver.method' use "
169 << parameter.keyPrefix() << "method instead\n";
170 }
171 else
172 kspType = static_cast< PetscSolver >(
173 parameter.solverMethod (
174 { SolverParameter::gmres,
175 SolverParameter::bicgstab,
176 SolverParameter::cg,
177 SolverParameter::minres,
178 SolverParameter::bicg,
179 SolverParameter::preonly
180 }, { "kspoptions" } )
181 );
182
183 if (kspType > PetscSolver::kspoptions)
184 solverName_ = SolverParameter::solverMethodTable( static_cast< int >( kspType ) );
185 else
186 solverName_ = "kspoptions";
187
188 // select linear solver
189 switch( kspType )
190 {
191 case PetscSolver::cg:
192 ::Dune::Petsc::KSPSetType( ksp(), KSPCG );
193 break;
194 case PetscSolver::bicgstab:
195 ::Dune::Petsc::KSPSetType( ksp(), KSPBCGS );
196 break;
197 case PetscSolver::gmres:
198 {
199 ::Dune::Petsc::KSPSetType( ksp(), KSPGMRES );
200 PetscInt restart = 10;
201 if( reader.exists("petsc.gmresrestart") )
202 {
203 restart = reader.getValue<int>("petsc.gmresrestart", restart );
204 std::cout << "WARNING: using deprecated parameter 'petsc.gmresrestart' use "
205 << parameter.keyPrefix() << "gmres.restart instead\n";
206 }
207 else
208 restart = parameter.gmresRestart() ;
209
210 ::Dune::Petsc::KSPGMRESSetRestart( ksp(), restart );
211 break;
212 }
213 case PetscSolver::minres:
214 ::Dune::Petsc::KSPSetType( ksp(), KSPMINRES );
215 break;
216 case PetscSolver::bicg:
217 ::Dune::Petsc::KSPSetType( ksp(), KSPBICG );
218 break;
219 case PetscSolver::preonly:
220 ::Dune::Petsc::KSPSetType( ksp(), KSPPREONLY );
221 break;
222 case PetscSolver::kspoptions:
223 // setup solver context from database/cmdline options
224 ::Dune::Petsc::KSPSetFromOptions( ksp() );
225 ::Dune::Petsc::KSPSetUp( ksp() );
226 break;
227 default:
228 DUNE_THROW(InvalidStateException,"PetscInverseOperator: invalid solver choosen." );
229 }
230
232 // preconditioning
234
235 int pcType = SolverParameter::none;
236 if( reader.exists("petsc.preconditioning.method") )
237 {
238 const std::string pcNames[] = { "default", "none", "asm", "sor", "jacobi", "ilu", "icc", "superlu", "hypre", "ml", "lu" };
239 pcType = reader.getEnum("petsc.preconditioning.method", pcNames, 0 );
240 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.method' use "
241 << parameter.keyPrefix() << "preconditioning.method instead\n";
242 if (pcType >= 8)
243 pcType = 7-pcType; // hypre=-1, ml=-2, lu=-3
244 }
245 else
246 {
247 pcType = parameter.preconditionMethod(
248 {
249 SolverParameter::none, // no preconditioning
250 SolverParameter::oas, // Overlapping Additive Schwarz
251 SolverParameter::sor, // SOR and SSOR
252 SolverParameter::jacobi, // Jacobi preconditioning
253 SolverParameter::ilu, // ILU preconditioning
254 SolverParameter::icc, // Incomplete Cholesky factorization
255 SolverParameter::superlu // SuperLU direct factorization
256 },
257 {"kspoptions", // = 0, // use command line options -ksp...
258 "hypre", // = -1, // Hypre preconditioning
259 "ml", // = -2, // ML preconditioner (from Trilinos)
260 "lu", // = -3, // LU factorization
261 "pcgamg", // = -4 // Petsc internal AMG
262 });
263 }
264
265 // setup preconditioning context
266 PC pc;
267 ::Dune::Petsc::KSPGetPC( ksp(), &pc );
268
269 switch( pcType )
270 {
271 case 0:
272 // don't setup the pc context twice
273 if ( kspType != PetscSolver::kspoptions )
274 {
275 // setup pc context from database/cmdline options
276 ::Dune::Petsc::PCSetFromOptions( pc );
277 ::Dune::Petsc::PCSetUp( pc );
278 }
279 break;
280 case SolverParameter::none:
281 ::Dune::Petsc::PCSetType( pc, PCNONE );
282 break;
283 case SolverParameter::oas:
284 {
285 ::Dune::Petsc::PCSetType( pc, PCASM );
286 ::Dune::Petsc::PCSetUp( pc );
287 break;
288 }
289 case SolverParameter::sor:
290 ::Dune::Petsc::PCSetType( pc, PCSOR );
291 break;
292 case SolverParameter::jacobi:
293 ::Dune::Petsc::PCSetType( pc, PCJACOBI );
294 break;
295 case -1: // PetscPrec::hypre:
296 {
297#ifdef PETSC_HAVE_HYPRE
298 int hypreType = parameter.hypreMethod();
299 std::string hypre;
300 if ( hypreType == PetscSolverParameter::boomeramg )
301 hypre = "boomeramg";
302 else if ( hypreType == PetscSolverParameter::parasails )
303 hypre = "parasails";
304 else if ( hypreType == PetscSolverParameter::pilut )
305 hypre = "pilut";
306 else
307 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid hypre preconditioner choosen." );
308
309 ::Dune::Petsc::PCSetType( pc, PCHYPRE );
310 ::Dune::Petsc::PCHYPRESetType( pc, hypre.c_str() );
311 ::Dune::Petsc::PCSetUp( pc );
312#else // PETSC_HAVE_HYPRE
313 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with hypre support." );
314#endif // PETSC_HAVE_HYPRE
315 break;
316 }
317 case -2: // PetscPrec::ml:
318#ifdef PETSC_HAVE_ML
319 ::Dune::Petsc::PCSetType( pc, PCML );
320#else // PETSC_HAVE_ML
321 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with ml support." );
322#endif // PETSC_HAVE_ML
323 break;
324 case SolverParameter::ilu:
325 {
326 if ( MPIManager::size() > 1 )
327 DUNE_THROW( InvalidStateException, "PetscInverseOperator: ilu preconditioner does not work in parallel." );
328
329 // get fill-in level
330 PetscInt pcLevel;
331 if( reader.exists("petsc.preconditioning.levels") )
332 {
333 pcLevel = reader.getValue<int>("petsc.preconditioning.levels", 0 );
334 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.levels' use "
335 << parameter.keyPrefix() << "preconditioning.level instead\n";
336 }
337 else
338 pcLevel = parameter.preconditionerLevel() ;
339
340 ::Dune::Petsc::PCSetType( pc, PCILU );
341 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
342 break;
343 }
344 ::Dune::Petsc::PCSetType( pc, PCML );
345 break;
346 case SolverParameter::icc:
347 {
348#ifdef PETSC_HAVE_ICC
349 if ( MPIManager::size() > 1 )
350 DUNE_THROW( InvalidStateException, "PetscInverseOperator: icc preconditioner does not worl in parallel." );
351
352 // get fill-in level
353 PetscInt pcLevel;
354 if( reader.exists("petsc.preconditioning.levels") )
355 {
356 pcLevel = reader.getValue<int>("petsc.preconditioning.levels", 0 );
357 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.levels' use "
358 << parameter.keyPrefix() << "preconditioning.level instead\n";
359 }
360 else
361 pcLevel = parameter.preconditionerLevel() ;
362
363
364 ::Dune::Petsc::PCSetType( pc, PCICC );
365 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
366#else // PETSC_HAVE_ICC
367 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with icc support." );
368#endif // PETSC_HAVE_ICC
369 break;
370 }
371 case -3: // PetscPrec::lu:
372 case SolverParameter::superlu:
373 {
374 enum class Factorization { petsc = 0, superlu = 1, mumps = 2 };
375 Factorization factorType = Factorization::superlu;
376 if (pcType != SolverParameter::superlu)
377 factorType = static_cast<Factorization>(parameter.superluMethod());
378
379 ::Dune::Petsc::PCSetType( pc, PCLU );
380
381 if ( factorType == Factorization::petsc )
382 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERPETSC );
383 else if ( factorType == Factorization::superlu )
384 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERSUPERLU_DIST );
385 else if ( factorType == Factorization::mumps )
386 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERMUMPS );
387 else
388 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid factorization package choosen." );
389
390 ::Dune::Petsc::PCSetUp( pc );
391 break;
392 }
393 case -4: // PetscPrec::pcgamg:
394 ::Dune::Petsc::PCSetType( pc, PCGAMG );
395 break;
396
397 default:
398 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid preconditioner choosen." );
399 }
400
401 // set monitor in verbose mode for all cores
402 // (and then check Parameter::verbose locally inside monitor)
403 if( parameter.verbose() )
404 {
405 ::Dune::Petsc::KSPView( comm, ksp() );
406 ::Dune::Petsc::KSPMonitorSet( ksp(), &monitor, PETSC_NULL, PETSC_NULL);
407 }
408 }
409
410 int apply( const PetscDiscreteFunctionType& arg, PetscDiscreteFunctionType& dest ) const
411 {
412 // need to have a 'distributed' destination vector for continuous spaces
413 if( dest.space().continuous() )
414 dest.dofVector().clearGhost();
415
416 // call PETSc solvers
417 ::Dune::Petsc::KSPSolve( *ksp_, *arg.petscVec() , *dest.petscVec() );
418
419 // a continuous solution is 'distributed' so need a communication here
420 if( dest.space().continuous() )
421 {
422 dest.communicate();
423 }
424
425 // get number of iterations
426 PetscInt its ;
427 ::Dune::Petsc::KSPGetIterationNumber( *ksp_, &its );
428 KSPConvergedReason reason;
429 ::Dune::Petsc::KSPGetConvergedReason( *ksp_, &reason );
430 if( parameter_->verbose() && Parameter::verbose() )
431 {
432 // list of reasons:
433 // https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPConvergedReason.html
434 std::cout << "Converged reason:" << reason << std::endl;
435 }
436
437 bool converged = int(reason) >= 0 ;
438
439 return (converged) ? its : -its;
440 }
441
442 protected:
443 KSP & ksp () { assert( ksp_ ); return *ksp_; }
444
445 using BaseType :: assembledOperator_;
446 using BaseType :: parameter_;
447 using BaseType :: iterations_;
448
449 std::unique_ptr< KSP, KSPDeleter > ksp_; // PETSc Krylov Space solver context
450
451 std::string solverName_;
452 };
453
455
456 } // namespace Fem
457
458} // namespace Dune
459
460#endif // #if HAVE_PETSC
461
462#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