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)
62 if( Parameter :: verbose () )
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();
128 ::Dune::Petsc::KSPCreate( comm, &ksp() );
131 Mat& A = assembledOperator_->exportMatrix();
132#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
133 ::Dune::Petsc::KSPSetOperators( ksp(), A, A, SAME_PRECONDITIONER);
135 ::Dune::Petsc::KSPSetOperators( ksp(), A, A );
139 ::Dune::Petsc::KSPSetInitialGuessNonzero( ksp(), PETSC_TRUE );
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);
147 ::Dune::Petsc::KSPSetTolerances(ksp(), tolerance, 1e-50, PETSC_DEFAULT, maxits);
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,
161 const auto& reader = parameter.parameter();
162 PetscSolver kspType = PetscSolver::gmres;
163 if( reader.exists(
"petsc.kspsolver.method") )
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";
172 kspType =
static_cast< PetscSolver
>(
173 parameter.solverMethod (
174 { SolverParameter::gmres,
175 SolverParameter::bicgstab,
177 SolverParameter::minres,
178 SolverParameter::bicg,
179 SolverParameter::preonly
180 }, {
"kspoptions" } )
183 if (kspType > PetscSolver::kspoptions)
184 solverName_ = SolverParameter::solverMethodTable(
static_cast< int >( kspType ) );
186 solverName_ =
"kspoptions";
191 case PetscSolver::cg:
192 ::Dune::Petsc::KSPSetType( ksp(), KSPCG );
194 case PetscSolver::bicgstab:
195 ::Dune::Petsc::KSPSetType( ksp(), KSPBCGS );
197 case PetscSolver::gmres:
199 ::Dune::Petsc::KSPSetType( ksp(), KSPGMRES );
200 PetscInt restart = 10;
201 if( reader.exists(
"petsc.gmresrestart") )
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";
208 restart = parameter.gmresRestart() ;
210 ::Dune::Petsc::KSPGMRESSetRestart( ksp(), restart );
213 case PetscSolver::minres:
214 ::Dune::Petsc::KSPSetType( ksp(), KSPMINRES );
216 case PetscSolver::bicg:
217 ::Dune::Petsc::KSPSetType( ksp(), KSPBICG );
219 case PetscSolver::preonly:
220 ::Dune::Petsc::KSPSetType( ksp(), KSPPREONLY );
222 case PetscSolver::kspoptions:
224 ::Dune::Petsc::KSPSetFromOptions( ksp() );
225 ::Dune::Petsc::KSPSetUp( ksp() );
228 DUNE_THROW(InvalidStateException,
"PetscInverseOperator: invalid solver choosen." );
235 int pcType = SolverParameter::none;
236 if( reader.exists(
"petsc.preconditioning.method") )
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";
247 pcType = parameter.preconditionMethod(
249 SolverParameter::none,
250 SolverParameter::oas,
251 SolverParameter::sor,
252 SolverParameter::jacobi,
253 SolverParameter::ilu,
254 SolverParameter::icc,
255 SolverParameter::superlu
267 ::Dune::Petsc::KSPGetPC( ksp(), &pc );
273 if ( kspType != PetscSolver::kspoptions )
276 ::Dune::Petsc::PCSetFromOptions( pc );
277 ::Dune::Petsc::PCSetUp( pc );
280 case SolverParameter::none:
281 ::Dune::Petsc::PCSetType( pc, PCNONE );
283 case SolverParameter::oas:
285 ::Dune::Petsc::PCSetType( pc, PCASM );
286 ::Dune::Petsc::PCSetUp( pc );
289 case SolverParameter::sor:
290 ::Dune::Petsc::PCSetType( pc, PCSOR );
292 case SolverParameter::jacobi:
293 ::Dune::Petsc::PCSetType( pc, PCJACOBI );
297#ifdef PETSC_HAVE_HYPRE
298 int hypreType = parameter.hypreMethod();
300 if ( hypreType == PetscSolverParameter::boomeramg )
302 else if ( hypreType == PetscSolverParameter::parasails )
304 else if ( hypreType == PetscSolverParameter::pilut )
307 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: invalid hypre preconditioner choosen." );
309 ::Dune::Petsc::PCSetType( pc, PCHYPRE );
310 ::Dune::Petsc::PCHYPRESetType( pc, hypre.c_str() );
311 ::Dune::Petsc::PCSetUp( pc );
313 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: petsc not build with hypre support." );
319 ::Dune::Petsc::PCSetType( pc, PCML );
321 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: petsc not build with ml support." );
324 case SolverParameter::ilu:
326 if ( MPIManager::size() > 1 )
327 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: ilu preconditioner does not work in parallel." );
331 if( reader.exists(
"petsc.preconditioning.levels") )
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";
338 pcLevel = parameter.preconditionerLevel() ;
340 ::Dune::Petsc::PCSetType( pc, PCILU );
341 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
344 ::Dune::Petsc::PCSetType( pc, PCML );
346 case SolverParameter::icc:
349 if ( MPIManager::size() > 1 )
350 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: icc preconditioner does not worl in parallel." );
354 if( reader.exists(
"petsc.preconditioning.levels") )
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";
361 pcLevel = parameter.preconditionerLevel() ;
364 ::Dune::Petsc::PCSetType( pc, PCICC );
365 ::Dune::Petsc::PCFactorSetLevels( pc, pcLevel );
367 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: petsc not build with icc support." );
372 case SolverParameter::superlu:
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());
379 ::Dune::Petsc::PCSetType( pc, PCLU );
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 );
388 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: invalid factorization package choosen." );
390 ::Dune::Petsc::PCSetUp( pc );
394 ::Dune::Petsc::PCSetType( pc, PCGAMG );
398 DUNE_THROW( InvalidStateException,
"PetscInverseOperator: invalid preconditioner choosen." );
403 if( parameter.verbose() )
405 ::Dune::Petsc::KSPView( comm, ksp() );
406 ::Dune::Petsc::KSPMonitorSet( ksp(), &monitor, PETSC_NULL, PETSC_NULL);
410 int apply(
const PetscDiscreteFunctionType& arg, PetscDiscreteFunctionType& dest )
const
413 if( dest.space().continuous() )
414 dest.dofVector().clearGhost();
417 ::Dune::Petsc::KSPSolve( *ksp_, *arg.petscVec() , *dest.petscVec() );
420 if( dest.space().continuous() )
427 ::Dune::Petsc::KSPGetIterationNumber( *ksp_, &its );
428 KSPConvergedReason reason;
429 ::Dune::Petsc::KSPGetConvergedReason( *ksp_, &reason );
430 if( parameter_->verbose() && Parameter::verbose() )
434 std::cout <<
"Converged reason:" << reason << std::endl;
437 bool converged = int(reason) >= 0 ;
439 return (converged) ? its : -its;
443 KSP & ksp () { assert( ksp_ );
return *ksp_; }
445 using BaseType :: assembledOperator_;
446 using BaseType :: parameter_;
447 using BaseType :: iterations_;
449 std::unique_ptr< KSP, KSPDeleter > ksp_;
451 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