1#ifndef DUNE_FEM_PETSCINVERSEOPERATORS_HH
2#define DUNE_FEM_PETSCINVERSEOPERATORS_HH
33 template<
class DF,
class Op = Dune::Fem::Operator< DF, DF > >
34 class PetscInverseOperator;
36 template <
class DF,
class Op >
37 struct PetscInverseOperatorTraits
40 typedef typename DF :: DiscreteFunctionSpaceType SpaceType ;
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;
53 template<
class DF,
class Op >
54 class PetscInverseOperator
55 :
public InverseOperatorInterface< PetscInverseOperatorTraits< DF, Op > >
60 monitor (KSP ksp, PetscInt it, PetscReal rnorm,
void *mctx)
64 std::cout <<
"PETSc::KSP: it = " << it <<
" res = " << rnorm << std::endl;
66 return PetscErrorCode(0);
72 void operator() ( KSP* p )
const
77 ::Dune::Petsc::KSPDestroy( p );
82 typedef PetscInverseOperatorTraits< DF, Op > Traits;
83 typedef InverseOperatorInterface< Traits > BaseType;
84 friend class InverseOperatorInterface< Traits >;
87 typedef typename BaseType :: SolverDiscreteFunctionType PetscDiscreteFunctionType;
88 typedef typename BaseType :: OperatorType OperatorType;
90 PetscInverseOperator (
const PetscSolverParameter ¶meter = PetscSolverParameter(
Parameter::container()) )
91 : BaseType( parameter )
94 PetscInverseOperator (
const OperatorType &op,
const PetscSolverParameter ¶meter = PetscSolverParameter(
Parameter::container()) )
95 : BaseType( parameter )
100 void bind (
const OperatorType &op )
102 BaseType :: bind( op );
103 initialize( *parameter_ );
108 BaseType :: unbind();
112 void printTexInfo(std::ostream& out)
const
114 out <<
"Solver: " << solverName_ <<
" eps = " << parameter_->tolerance() ;
119 void initialize (
const PetscSolverParameter& parameter )
121 if( !assembledOperator_ )
122 DUNE_THROW(NotImplemented,
"Petsc solver with matrix free implementations not yet supported!");
125 ksp_.reset(
new KSP() );
126 const auto& comm = assembledOperator_->domainSpace().gridPart().comm();
127 ::Dune::Petsc::KSPCreate( comm, &ksp() );
130 Mat& A = assembledOperator_->exportMatrix();
131#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
132 ::Dune::Petsc::KSPSetOperators( ksp(), A, A, SAME_PRECONDITIONER);
134 ::Dune::Petsc::KSPSetOperators( ksp(), A, A );
138 ::Dune::Petsc::KSPSetInitialGuessNonzero( ksp(), PETSC_TRUE );
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);
146 ::Dune::Petsc::KSPSetTolerances(ksp(), tolerance, 1e-50, PETSC_DEFAULT, maxits);
148 enum class PetscSolver {
160 const auto& reader = parameter.parameter();
162 if( reader.exists(
"petsc.kspsolver.method") )
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";
171 kspType =
static_cast< PetscSolver
>(
172 parameter.solverMethod (
179 }, {
"kspoptions" } )
182 if (kspType > PetscSolver::kspoptions)
185 solverName_ =
"kspoptions";
191 ::Dune::Petsc::KSPSetType( ksp(), KSPCG );
194 ::Dune::Petsc::KSPSetType( ksp(), KSPBCGS );
198 ::Dune::Petsc::KSPSetType( ksp(), KSPGMRES );
199 PetscInt restart = 10;
200 if( reader.exists(
"petsc.gmresrestart") )
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";
207 restart = parameter.gmresRestart() ;
209 ::Dune::Petsc::KSPGMRESSetRestart( ksp(), restart );
212 case PetscSolver::minres:
213 ::Dune::Petsc::KSPSetType( ksp(), KSPMINRES );
215 case PetscSolver::bicg:
216 ::Dune::Petsc::KSPSetType( ksp(), KSPBICG );
218 case PetscSolver::preonly:
219 ::Dune::Petsc::KSPSetType( ksp(), KSPPREONLY );
221 case PetscSolver::kspoptions:
223 ::Dune::Petsc::KSPSetFromOptions( ksp() );
224 ::Dune::Petsc::KSPSetUp( ksp() );
227 DUNE_THROW(InvalidStateException,
"PetscInverseOperator: invalid solver choosen." );
235 if( reader.exists(
"petsc.preconditioning.method") )
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";
246 pcType = parameter.preconditionMethod(
266 ::Dune::Petsc::KSPGetPC( ksp(), &pc );
272 if ( kspType != PetscSolver::kspoptions )
275 ::Dune::Petsc::PCSetFromOptions( pc );
276 ::Dune::Petsc::PCSetUp( pc );
280 ::Dune::Petsc::PCSetType( pc, PCNONE );
284 ::Dune::Petsc::PCSetType( pc, PCASM );
285 ::Dune::Petsc::PCSetUp( pc );
289 ::Dune::Petsc::PCSetType( pc, PCSOR );
292 ::Dune::Petsc::PCSetType( pc, PCJACOBI );
296#ifdef PETSC_HAVE_HYPRE
297 int hypreType = parameter.hypreMethod();
299 if ( hypreType == PetscSolverParameter::boomeramg )
301 else if ( hypreType == PetscSolverParameter::parasails )
303 else if ( hypreType == PetscSolverParameter::pilut )
306 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: invalid hypre preconditioner choosen." );
308 ::Dune::Petsc::PCSetType( pc, PCHYPRE );
309 ::Dune::Petsc::PCHYPRESetType( pc, hypre.c_str() );
310 ::Dune::Petsc::PCSetUp( pc );
312 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: petsc not build with hypre support." );
318 ::Dune::Petsc::PCSetType( pc, PCML );
320 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: petsc not build with ml support." );
326 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: ilu preconditioner does not work in parallel." );
330 if( reader.exists(
"petsc.preconditioning.levels") )
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";
337 pcLevel = parameter.preconditionerLevel() ;
339 ::Dune::Petsc::PCSetType( pc, PCILU );
340 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
343 ::Dune::Petsc::PCSetType( pc, PCML );
349 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: icc preconditioner does not worl in parallel." );
353 if( reader.exists(
"petsc.preconditioning.levels") )
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";
360 pcLevel = parameter.preconditionerLevel() ;
363 ::Dune::Petsc::PCSetType( pc, PCICC );
364 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
366 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: petsc not build with icc support." );
373 enum class Factorization { petsc = 0, superlu = 1, mumps = 2 };
374 Factorization factorType = Factorization::superlu;
376 factorType =
static_cast<Factorization
>(parameter.superluMethod());
378 ::Dune::Petsc::PCSetType( pc, PCLU );
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 );
387 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: invalid factorization package choosen." );
389 ::Dune::Petsc::PCSetUp( pc );
393 ::Dune::Petsc::PCSetType( pc, PCGAMG );
397 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: invalid preconditioner choosen." );
403 ::Dune::Petsc::KSPView( comm, ksp() );
404 ::Dune::Petsc::KSPMonitorSet( ksp(), &monitor, PETSC_NULL, PETSC_NULL);
408 int apply(
const PetscDiscreteFunctionType& arg, PetscDiscreteFunctionType& dest )
const
411 if( dest.space().continuous() )
412 dest.dofVector().clearGhost();
415 ::Dune::Petsc::KSPSolve( *ksp_, *arg.petscVec() , *dest.petscVec() );
418 if( dest.space().continuous() )
425 ::Dune::Petsc::KSPGetIterationNumber( *ksp_, &its );
426 KSPConvergedReason reason;
427 ::Dune::Petsc::KSPGetConvergedReason( *ksp_, &reason );
432 std::cout <<
"Converged reason:" << reason << std::endl;
435 bool converged = int(reason) >= 0 ;
437 return (converged) ? its : -its;
441 KSP & ksp () { assert( ksp_ );
return *ksp_; }
443 using BaseType :: assembledOperator_;
444 using BaseType :: parameter_;
445 using BaseType :: iterations_;
447 std::unique_ptr< KSP, KSPDeleter > ksp_;
449 std::string solverName_;
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