DUNE-FEM (unstable)

petscinverseoperators.hh
1#ifndef DUNE_FEM_PETSCINVERSEOPERATORS_HH
2#define DUNE_FEM_PETSCINVERSEOPERATORS_HH
3
4#include <limits>
5
6#include <dune/fem/function/common/scalarproducts.hh>
7#include <dune/fem/operator/common/operator.hh>
8#include <dune/fem/io/parameter.hh>
9#include <dune/fem/solver/inverseoperatorinterface.hh>
10
11
12#if HAVE_PETSC
13#include <dune/fem/operator/linear/petscoperator.hh>
14#include <dune/fem/misc/petsc/petsccommon.hh>
15#include <dune/fem/function/petscdiscretefunction.hh>
16#include <dune/fem/solver/parameter.hh>
17#include <dune/fem/solver/petscavailable.hh>
18
19namespace Dune
20{
21
22 namespace Fem
23 {
24
25 //=====================================================================
26 // Implementation for PETSc matrix based Krylov solvers
27 //=====================================================================
28
33 // PETScSolver
34 // --------------
35 template< class DF, class Op = Dune::Fem::Operator< DF, DF > >
36 class PetscInverseOperator;
37
38 template <class DF, class Op >
39 struct PetscInverseOperatorTraits
40 {
41 private:
42 typedef typename DF :: DiscreteFunctionSpaceType SpaceType ;
43 public:
44 typedef DF DiscreteFunctionType;
45 typedef Op OperatorType;
46 typedef OperatorType PreconditionerType;
47 typedef PetscDiscreteFunction< SpaceType > SolverDiscreteFunctionType;
48 typedef PetscLinearOperator< DF, DF > AssembledOperatorType;
49 typedef PetscInverseOperator< DF, Op > InverseOperatorType;
50 typedef PetscSolverParameter SolverParameterType;
51 };
52
53
55 template< class DF, class Op >
56 class PetscInverseOperator
57 : public InverseOperatorInterface< PetscInverseOperatorTraits< DF, Op > >,
58 public PetscInverseOperatorAvailable
59 {
60 protected:
61 // monitor function for PETSc solvers
62 static PetscErrorCode
63 monitor (KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
64 {
65 if( Parameter :: verbose ( Parameter::solverStatistics ) )
66 {
67 std::cout << "PETSc::KSP: it = "
68 << std::setw(3) << std::left << it
69 << " res = " << rnorm << std::endl;
70 }
71 return PetscErrorCode(0);
72 }
73
74 // destroy solver context
75 struct KSPDeleter
76 {
77 void operator() ( KSP* p ) const
78 {
79 if( !p )
80 return;
81
82 ::Dune::Petsc::KSPDestroy( p );
83 delete p;
84 }
85 };
86
87 typedef PetscInverseOperatorTraits< DF, Op > Traits;
88 typedef InverseOperatorInterface< Traits > BaseType;
89 friend class InverseOperatorInterface< Traits >;
90
91 using PetscInverseOperatorAvailable::PetscSolver;
92
93 public:
94 using PetscInverseOperatorAvailable::supportedSolverMethods;
95 using PetscInverseOperatorAvailable::supportedPreconditionMethods;
96 using PetscInverseOperatorAvailable::extraPreconditionMethods;
97
101 static const bool preconditioningAvailable = false;
102
103 typedef typename BaseType :: SolverDiscreteFunctionType PetscDiscreteFunctionType;
104 typedef typename BaseType :: OperatorType OperatorType;
105 typedef typename BaseType :: PreconditionerType PreconditionerType;
106
107 PetscInverseOperator ( const PetscSolverParameter &parameter = PetscSolverParameter(Parameter::container()) )
108 : BaseType( parameter )
109 {}
110
111 PetscInverseOperator ( const OperatorType &op, const PetscSolverParameter &parameter = PetscSolverParameter(Parameter::container()) )
112 : BaseType( parameter )
113 {
114 bind( op );
115 }
116
117 void bind ( const OperatorType &op )
118 {
119 BaseType :: bind( op );
120 initialize( *parameter_ );
121 }
122
123 void bind ( const OperatorType &op, const PreconditionerType& p )
124 {
125 DUNE_THROW(NotImplemented,"PetscInverseOperator::bind: preconditioners cannot be set from outside!");
126 }
127
128 void unbind ()
129 {
130 BaseType :: unbind();
131 ksp_.reset();
132 }
133
134 void printTexInfo(std::ostream& out) const
135 {
136 out << "Solver: " << solverName_ << " eps = " << parameter_->tolerance() ;
137 out << "\\\\ \n";
138 }
139
140 protected:
141 template <class Reader>
142 bool parseAndSetOptions(const Reader& reader, const std::string& optionsKey) const
143 {
144 if( reader.exists(optionsKey) )
145 {
146 std::string opts = reader.template getValue< std::string > (optionsKey);
147
148 std::stringstream options( opts );
149
150 std::string key;
151 std::string value;
152 std::string nextKey;
153 bool readKey = true;
154
155 while (! options.eof() )
156 {
157 if( readKey )
158 options >> key;
159 else
160 {
161 key = nextKey;
162 readKey = true;
163 }
164
165 if( key.length() == 0 )
166 break;
167
168 options >> value;
169 // if a key was instead simply store for next iteration
170 if ( value.length() > 0 && value[0] == '-')
171 {
172 // store key for next iteration
173 nextKey = value;
174 value = ""; // empty value
175 readKey = false;
176 }
177
178 // add key,value to global options data base
179 ::Dune::Petsc::PetscOptionsSetValue( key, value);
180 }
181 return true;
182 }
183 else
184 {
185 std::cout << "WARNING: 'kspoptions' selected but no options delivered with '" << optionsKey << "'!" << std::endl;
186 return false;
187 }
188 }
189
190
191 void initialize ( const PetscSolverParameter& parameter )
192 {
193 if( !assembledOperator_ )
194 DUNE_THROW(NotImplemented,"Petsc solver with matrix free implementations not yet supported!");
195
196 // Create linear solver context
197 ksp_.reset( new KSP() );
198 const auto& comm = assembledOperator_->domainSpace().gridPart().comm();
199
200 ::Dune::Petsc::KSPCreate( comm, &ksp() );
201
202 // attach Matrix to linear solver context
203 Mat& A = assembledOperator_->exportMatrix();
204#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
205 ::Dune::Petsc::KSPSetOperators( ksp(), A, A, SAME_PRECONDITIONER);
206#else
207 ::Dune::Petsc::KSPSetOperators( ksp(), A, A );
208#endif
209
210 // allow for non-zero initial guess
211 ::Dune::Petsc::KSPSetInitialGuessNonzero( ksp(), PETSC_TRUE );
212
213 // set prescribed tolerances
214 PetscInt maxits = parameter_->maxIterations();
215 PetscReal tolerance = parameter_->tolerance();
216 PetscReal omega = parameter_->relaxation();
217 if (parameter_->errorMeasure() == 0)
218 ::Dune::Petsc::KSPSetTolerances(ksp(), 1e-50, tolerance, PETSC_DEFAULT, maxits);
219 else
220 ::Dune::Petsc::KSPSetTolerances(ksp(), tolerance, 1e-50, PETSC_DEFAULT, maxits);
221
222 // if special petsc solver parameter exists use that one, otherwise
223 // use solverMethod from SolverParameter
224 const auto& reader = parameter.parameter();
225 PetscSolver kspType = PetscSolver::gmres;
226 if( reader.exists("petsc.kspsolver.method") )
227 {
228 // see PETSc docu for more types
229 const std::string kspNames[] = { "default", "cg", "bicgstab", "gmres", "minres", "gradient", "loop", "superlu", "bicg", "preonly" };
230 kspType = static_cast< PetscSolver >( reader.getEnum("petsc.kspsolver.method", kspNames, int(PetscSolver::gmres) ) );
231 std::cout << "WARNING: using deprecated parameter 'petsc.kpsolver.method' use "
232 << parameter.keyPrefix() << "method instead\n";
233 }
234 else
235 kspType = static_cast< PetscSolver >(
236 parameter.solverMethod (
237 supportedSolverMethods(), { "kspoptions" } )
238 );
239
240 if (kspType > PetscSolver::kspoptions)
241 solverName_ = SolverParameter::solverMethodTable( static_cast< int >( kspType ) );
242 else
243 solverName_ = "kspoptions";
244
245 // select linear solver
246 switch( kspType )
247 {
248 case PetscSolver::cg:
249 ::Dune::Petsc::KSPSetType( ksp(), KSPCG );
250 break;
251 case PetscSolver::bicgstab:
252 ::Dune::Petsc::KSPSetType( ksp(), KSPBCGS );
253 break;
254 case PetscSolver::gmres:
255 {
256 ::Dune::Petsc::KSPSetType( ksp(), KSPGMRES );
257 PetscInt restart = 10;
258 if( reader.exists("petsc.gmresrestart") )
259 {
260 restart = reader.getValue<int>("petsc.gmresrestart", restart );
261 std::cout << "WARNING: using deprecated parameter 'petsc.gmresrestart' use "
262 << parameter.keyPrefix() << "gmres.restart instead\n";
263 }
264 else
265 restart = parameter.gmresRestart() ;
266
267 ::Dune::Petsc::KSPGMRESSetRestart( ksp(), restart );
268 break;
269 }
270 case PetscSolver::minres:
271 ::Dune::Petsc::KSPSetType( ksp(), KSPMINRES );
272 break;
273 case PetscSolver::bicg:
274 ::Dune::Petsc::KSPSetType( ksp(), KSPBICG );
275 break;
276 case PetscSolver::preonly:
277 ::Dune::Petsc::KSPSetType( ksp(), KSPPREONLY );
278 break;
279 case PetscSolver::kspoptions:
280 {
281 const std::string key = parameter.keyPrefix() + "kspoptions";
282 const bool foundOptions = parseAndSetOptions( reader, key );
283
284 // setup solver context from database/cmdline options
285 ::Dune::Petsc::KSPSetFromOptions( ksp() );
286 ::Dune::Petsc::KSPSetUp( ksp() );
287 if( foundOptions )
288 {
289 // clear global options data base
290 ::Dune::Petsc::PetscOptionsClear();
291 }
292
293 break;
294 }
295 default:
296 DUNE_THROW(InvalidStateException,"PetscInverseOperator: invalid solver choosen." );
297 }
298 if ( parameter_->knollTrick() )
299 {
300 ::Dune::Petsc::KSPSetInitialGuessKnoll( ksp(), PETSC_TRUE );
301 }
302
304 // preconditioning
306
307 int pcType = SolverParameter::none;
308 if( reader.exists("petsc.preconditioning.method") )
309 {
310 const std::string pcNames[] = { "default", "none", "asm", "sor", "jacobi", "ilu", "icc", "superlu", "hypre", "ml", "lu" };
311 pcType = reader.getEnum("petsc.preconditioning.method", pcNames, 0 );
312 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.method' use "
313 << parameter.keyPrefix() << ".preconditioning.method instead\n";
314 if (pcType >= 8)
315 pcType = 7-pcType; // hypre=-1, ml=-2, lu=-3
316 }
317 else
318 {
319 pcType = parameter.preconditionMethod(
320 supportedPreconditionMethods(),
321 extraPreconditionMethods() );
322 }
323
324 // setup preconditioning context
325 PC pc;
326 ::Dune::Petsc::KSPGetPC( ksp(), &pc );
327
328 switch( pcType )
329 {
330 case 0:
331 // don't setup the pc context twice
332 if ( kspType != PetscSolver::kspoptions )
333 {
334 const std::string key = parameter.keyPrefix() + "kspoptions";
335 const bool foundOptions = parseAndSetOptions( reader, key );
336
337 // setup pc context from database/cmdline options
338 ::Dune::Petsc::PCSetFromOptions( pc );
339 ::Dune::Petsc::PCSetUp( pc );
340 if( foundOptions )
341 {
342 // clear global options data base
343 ::Dune::Petsc::PetscOptionsClear();
344 }
345 }
346 break;
347 case SolverParameter::none:
348 ::Dune::Petsc::PCSetType( pc, PCNONE );
349 break;
350 case SolverParameter::oas:
351 {
352 ::Dune::Petsc::PCSetType( pc, PCASM );
353 ::Dune::Petsc::PCSetUp( pc );
354 break;
355 }
356 case SolverParameter::gauss_seidel:
357 ::Dune::Petsc::PCSetType( pc, PCSOR );
358 ::Dune::Petsc::PCSORSetOmega( pc, 1.0 );
359 break;
360 case SolverParameter::sor:
361 ::Dune::Petsc::PCSetType( pc, PCSOR );
362 ::Dune::Petsc::PCSORSetOmega( pc, omega );
363 break;
364 case SolverParameter::ssor:
365 ::Dune::Petsc::PCSetType( pc, PCSOR );
366 // set symmetric version
367 ::Dune::Petsc::PCSORSetSymmetric( pc, SOR_LOCAL_SYMMETRIC_SWEEP );
368 ::Dune::Petsc::PCSORSetOmega( pc, omega );
369 break;
370 case SolverParameter::jacobi:
371 ::Dune::Petsc::PCSetType( pc, PCJACOBI );
372 break;
373 case -1: // PetscPrec::hypre:
374 {
375#ifdef PETSC_HAVE_HYPRE
376 // set with parameter ...petsc.preconditioning.hypre.method
377 int hypreType = parameter.hypreMethod();
378 std::string hypre;
379 if ( hypreType == PetscSolverParameter::boomeramg )
380 hypre = "boomeramg";
381 else if ( hypreType == PetscSolverParameter::parasails )
382 hypre = "parasails";
383 else if ( hypreType == PetscSolverParameter::pilut )
384 hypre = "pilut";
385 else
386 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid hypre preconditioner choosen." );
387
388 ::Dune::Petsc::PCSetType( pc, PCHYPRE );
389 ::Dune::Petsc::PCHYPRESetType( pc, hypre.c_str() );
390 ::Dune::Petsc::PCSetUp( pc );
391#else // PETSC_HAVE_HYPRE
392 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with hypre support." );
393#endif // PETSC_HAVE_HYPRE
394 break;
395 }
396 case -2: // PetscPrec::ml:
397#ifdef PETSC_HAVE_ML
398 ::Dune::Petsc::PCSetType( pc, PCML );
399#else // PETSC_HAVE_ML
400 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with ml support." );
401#endif // PETSC_HAVE_ML
402 break;
403 case SolverParameter::ilu:
404 {
405 if ( MPIManager::size() > 1 )
406 DUNE_THROW( InvalidStateException, "PetscInverseOperator: ilu preconditioner does not work in parallel." );
407
408 // get fill-in level
409 PetscInt pcLevel;
410 if( reader.exists("petsc.preconditioning.levels") )
411 {
412 pcLevel = reader.getValue<int>("petsc.preconditioning.levels", 0 );
413 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.levels' use "
414 << parameter.keyPrefix() << "preconditioning.level instead\n";
415 }
416 else
417 pcLevel = parameter.preconditionerLevel() ;
418
419 ::Dune::Petsc::PCSetType( pc, PCILU );
420 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
421 break;
422 }
423 ::Dune::Petsc::PCSetType( pc, PCML );
424 break;
425 case SolverParameter::icc:
426 {
427#ifdef PETSC_HAVE_ICC
428 if ( MPIManager::size() > 1 )
429 DUNE_THROW( InvalidStateException, "PetscInverseOperator: icc preconditioner does not worl in parallel." );
430
431 // get fill-in level
432 PetscInt pcLevel;
433 if( reader.exists("petsc.preconditioning.levels") )
434 {
435 pcLevel = reader.getValue<int>("petsc.preconditioning.levels", 0 );
436 std::cout << "WARNING: using deprecated parameter 'petsc.preconditioning.levels' use "
437 << parameter.keyPrefix() << "preconditioning.level instead\n";
438 }
439 else
440 pcLevel = parameter.preconditionerLevel() ;
441
442
443 ::Dune::Petsc::PCSetType( pc, PCICC );
444 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
445#else // PETSC_HAVE_ICC
446 DUNE_THROW( InvalidStateException, "PetscInverseOperator: petsc not build with icc support." );
447#endif // PETSC_HAVE_ICC
448 break;
449 }
450 case -3: // PetscPrec::lu:
451 case SolverParameter::superlu:
452 {
453 enum class Factorization { petsc = 0, superlu = 1, mumps = 2 };
454 Factorization factorType = Factorization::superlu;
455 if (pcType != SolverParameter::superlu)
456 factorType = static_cast<Factorization>(parameter.superluMethod());
457
458 ::Dune::Petsc::PCSetType( pc, PCLU );
459
460 if ( factorType == Factorization::petsc )
461 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERPETSC );
462 else if ( factorType == Factorization::superlu )
463 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERSUPERLU_DIST );
464 else if ( factorType == Factorization::mumps )
465 ::Dune::Petsc::PCFactorSetMatSolverPackage( pc, MATSOLVERMUMPS );
466 else
467 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid factorization package choosen." );
468
469 ::Dune::Petsc::PCSetUp( pc );
470 break;
471 }
472 case -4: // PetscPrec::pcgamg:
473 // requires MATRIX_AIJ, i.e. not blocking of entries
474 if( parameter.blockedMode() )
475 {
476 DUNE_THROW(NotImplemented,"PetscInverseOperator: 'pcgamg' requires 'aij' matrix. Set 'petsc.blockedmode' to false!");
477 }
478 ::Dune::Petsc::PCSetType( pc, PCGAMG );
479 break;
480
481 default:
482 DUNE_THROW( InvalidStateException, "PetscInverseOperator: invalid preconditioner choosen." );
483 }
484
485 // set monitor in verbose mode for all cores
486 // (and then check Parameter::verbose locally inside monitor)
487 if( parameter.verbose() && Parameter::verbose( Parameter::solverStatistics ) )
488 {
489 // only print information about solver type and pc type in extended mode
490 if( Parameter::verbose( Parameter::extendedStatistics ) )
491 ::Dune::Petsc::KSPView( comm, ksp() );
492
493 ::Dune::Petsc::KSPMonitorSet( ksp(), &monitor, PETSC_NULLPTR, PETSC_NULLPTR);
494 }
495 }
496
497 int apply( const PetscDiscreteFunctionType& arg, PetscDiscreteFunctionType& dest ) const
498 {
499 // need to have a 'distributed' destination vector for continuous spaces
500 if( dest.space().continuous() )
501 dest.dofVector().clearGhost();
502
503 // call PETSc solvers
504 ::Dune::Petsc::KSPSolve( *ksp_, *arg.petscVec() , *dest.petscVec() );
505
506 // a continuous solution is 'distributed' so need a communication here
507 if( dest.space().continuous() )
508 {
509 dest.communicate();
510 }
511
512 // get number of iterations
513 PetscInt its ;
514 ::Dune::Petsc::KSPGetIterationNumber( *ksp_, &its );
515 KSPConvergedReason reason;
516 ::Dune::Petsc::KSPGetConvergedReason( *ksp_, &reason );
517
518 bool converged = int(reason) >= 0 ;
519
520 if( parameter_->verbose() && Parameter::verbose( 1 ) )
521 {
522 // list of reasons:
523 // https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPConvergedReason.html
524 if( converged )
525 std::cout << "Converged reason: ";
526 else
527 std::cout << "**Diverged** reason: ";
528 std::cout << reason << " linear iterations: " << its << std::endl;
529 }
530
531 return (converged) ? its : -its;
532 }
533
534 protected:
535 KSP & ksp () { assert( ksp_ ); return *ksp_; }
536
537 using BaseType :: assembledOperator_;
538 using BaseType :: parameter_;
539 using BaseType :: iterations_;
540
541 std::unique_ptr< KSP, KSPDeleter > ksp_; // PETSc Krylov Space solver context
542
543 std::string solverName_;
544 };
545
547
548 } // namespace Fem
549
550} // namespace Dune
551
552#endif // #if HAVE_PETSC
553
554#endif // #ifndef DUNE_FEM_PETSCINVERSEOPERATORS_HH
Dune::Fem::CGInverseOperator
Inverse operator base on CG method. Uses a runtime parameter fem.preconditioning which enables diagon...
Definition: cginverseoperator.hh:408
Dune::Fem::FunctionSpace
A vector valued function space.
Definition: functionspace.hh:60
Dune::Fem::Parameter::verbose
static bool verbose()
obtain the cached value for fem.verbose with default verbosity level 2
Definition: parameter.hh:466
Dune::Fem::TupleDiscreteFunction
forward declaration
Definition: discretefunction.hh:51
DUNE_THROW
#define DUNE_THROW(E,...)
Definition: exceptions.hh:314
Dune::GeometryTypes::none
constexpr GeometryType none(unsigned int dim)
Returns a GeometryType representing a singular of dimension dim.
Definition: type.hh:471
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Nov 3, 23:36, 2025)