dune-fem  2.8-git
amgxsolver.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_AMGXSOLVER_HH
2 #define DUNE_FEM_AMGXSOLVER_HH
3 
4 #include <limits>
5 
9 
15 
16 #if HAVE_AMGXSOLVER
17 // AMGX solver wrapper based on Petsc data structures
18 #include <AmgXSolver.hpp>
19 #endif
20 
21 namespace Dune
22 {
23 
24  namespace Fem
25  {
26 
27  //=====================================================================
28  // Implementation of AMGX solver wrapper using PETSc matrix
29  //=====================================================================
30 
35  struct AMGXSolverParameter : public LocalParameter< SolverParameter, AMGXSolverParameter >
36  {
38 
39  public:
42 
44  : BaseType( parameter )
45  {}
46 
47 
50  {}
51 
54  {}
55 
56  virtual std::string solvermode () const
57  {
58  const std::string modes [] = { "dDDI" , "dDFI", "dFFI", "hDDI", "hDFI", "hFFI" };
59  int mode = parameter().getEnum(keyPrefix() + "amgx.mode", modes, 0 );
60  return modes[ mode ];
61  }
62 
63  virtual std::string solverconfig () const
64  {
65  return parameter().template getValue< std::string >( keyPrefix() + "amgx.config", "amgxconfig.json");
66  }
67  };
68 
69 
70 
71  // AMGXSolver
72  // --------------
73 
74  template< class DiscreteFunction >
75  class AMGXInverseOperator;
76 
77  template< class DiscreteFunction >
79  {
80  typedef DiscreteFunction DiscreteFunctionType;
82 
85 
86 #if HAVE_AMGXSOLVER
87  typedef Fem::PetscLinearOperator< DiscreteFunction, DiscreteFunction > AssembledOperatorType;
88 #else
90 #endif
91 
93 
95  };
96 
97 
98 
99 
101  template< class DF >
103  : public InverseOperatorInterface< AMGXInverseOperatorTraits< DF > >
104  {
107  friend class InverseOperatorInterface< Traits >;
108  public:
109  using BaseType :: parameter;
110 
115 
121  : BaseType( parameter )
122  {
123  }
124 
128  {
129  bind( op );
130  }
131 
135  {
136  bind( op, preconditioner );
137  }
138 
140  : AMGXInverseOperator( other.parameter() )
141  {
142  if( other.operator_ )
143  bind( *(other.operator_) );
144  }
145 
146  void bind( const OperatorType& op )
147  {
148  BaseType::bind( op );
149  init( parameter() );
150  }
151 
152  void unbind()
153  {
154 #if HAVE_AMGXSOLVER
155  amgXSolver_->finalize();
156  amgXSolver_.reset();
157 #endif
159  }
160 
161  protected:
163  {
164  if( assembledOperator_ )
165  {
166  std::string mode = parameter.solvermode();
167  std::string config = parameter.solverconfig();
168 #if HAVE_AMGXSOLVER
169  amgXSolver_.reset( new AmgXSolver() );
170  amgXSolver_->initialize(PETSC_COMM_WORLD, mode, config );
171 
172  // check that PetscMat was assembled not in block mode
173  if( assembledOperator_->blockedMode() )
174  DUNE_THROW(InvalidStateException, "AMGXInverseOperator only works with PetscLinearOperator in non-blocked mode!");
175 
176  // attach Matrix to linear solver context
177  Mat& A = const_cast<Mat &> (assembledOperator_->exportMatrix());
178 
179  // set matrix
180  amgXSolver_->setA( A );
181 #else
182  DUNE_THROW(InvalidStateException,"AMGX solver or PETSc not found during cmake config. Please reconfigure!");
183 #endif
184  }
185  }
186 
188  {
189  if( !assembledOperator_ )
190  DUNE_THROW(NotImplemented,"AMGX solver with matrix free implementations is not supported!");
191 
192 
193  int iterations = -1;
194 #if HAVE_AMGXSOLVER
195  assert( amgXSolver_ );
196 
197  // need to have a 'distributed' destination vector for continuous spaces
198  if( dest.space().continuous() )
199  dest.dofVector().clearGhost();
200 
201  // call PETSc solvers, dest = x, arg = rhs
202  Vec& x = *dest.petscVec();
203  Vec& rhs = *(const_cast< SolverDiscreteFunctionType& > (arg).petscVec());
204  amgXSolver_->solve( x, rhs );
205 
206  // a continuous solution is 'distributed' so need a communication here
207  if( dest.space().continuous() )
208  {
209  dest.communicate();
210  }
211 
212  // get number of iterations
213  amgXSolver_->getIters( iterations );
214 #else
215  DUNE_THROW(InvalidStateException,"AMGX solver or PETSc not found during cmake config. Please reconfigure!");
216 #endif
217  return iterations;
218  }
219 
220  protected:
221 #if HAVE_AMGXSOLVER
222  mutable std::unique_ptr< AmgXSolver > amgXSolver_;
223 #endif
226  };
227 
229 
230  } // namespace Fem
231 
232 } // namespace Dune
233 
234 #endif // #ifndef DUNE_FEM_PETSCSOLVER_HH
Definition: bindguard.hh:11
int getEnum(const std::string &key, const std::string(&values)[n]) const
Definition: reader.hh:225
static ParameterContainer & container()
Definition: io/parameter.hh:193
Definition: io/parameter.hh:551
Definition: amgxsolver.hh:36
AMGXSolverParameter(const SolverParameter &sp)
Definition: amgxsolver.hh:52
virtual std::string solvermode() const
Definition: amgxsolver.hh:56
virtual std::string solverconfig() const
Definition: amgxsolver.hh:63
LocalParameter< SolverParameter, AMGXSolverParameter > BaseType
Definition: amgxsolver.hh:37
AMGXSolverParameter(const std::string &keyPrefix, const ParameterReader &parameter=Parameter::container())
Definition: amgxsolver.hh:48
AMGXSolverParameter(const ParameterReader &parameter=Parameter::container())
Definition: amgxsolver.hh:43
AMGX solver context for PETSc Mat and PETSc Vec.
Definition: amgxsolver.hh:104
BaseType::SolverDiscreteFunctionType SolverDiscreteFunctionType
Definition: amgxsolver.hh:111
void init(const AMGXSolverParameter &parameter)
Definition: amgxsolver.hh:162
AMGXInverseOperator(const AMGXSolverParameter &parameter=AMGXSolverParameter())
constructor
Definition: amgxsolver.hh:120
int apply(const SolverDiscreteFunctionType &arg, SolverDiscreteFunctionType &dest) const
Definition: amgxsolver.hh:187
void bind(const OperatorType &op)
Definition: amgxsolver.hh:146
AMGXInverseOperator(const AMGXInverseOperator &other)
Definition: amgxsolver.hh:139
void unbind()
Definition: amgxsolver.hh:152
BaseType::PreconditionerType PreconditionerType
Definition: amgxsolver.hh:113
AMGXInverseOperator(const OperatorType &op, const AMGXSolverParameter &parameter=AMGXSolverParameter())
Definition: amgxsolver.hh:125
BaseType::AssembledOperatorType AssembledOperatorType
Definition: amgxsolver.hh:114
AMGXInverseOperator(const OperatorType &op, PreconditionerType &preconditioner, const AMGXSolverParameter &parameter=AMGXSolverParameter())
Definition: amgxsolver.hh:132
BaseType::OperatorType OperatorType
Definition: amgxsolver.hh:112
Definition: amgxsolver.hh:79
OperatorType PreconditionerType
Definition: amgxsolver.hh:84
DiscreteFunction DiscreteFunctionType
Definition: amgxsolver.hh:80
PetscDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType
Definition: amgxsolver.hh:81
Dune::Fem::Operator< DiscreteFunction, DiscreteFunction > OperatorType
Definition: amgxsolver.hh:83
OperatorType AssembledOperatorType
Definition: amgxsolver.hh:89
AMGXSolverParameter SolverParameterType
Definition: amgxsolver.hh:94
AMGXInverseOperator< DiscreteFunction > InverseOperatorType
Definition: amgxsolver.hh:92
Definition: inverseoperatorinterface.hh:17
Traits ::PreconditionerType PreconditionerType
Definition: inverseoperatorinterface.hh:32
const OperatorType * operator_
Definition: inverseoperatorinterface.hh:197
void unbind()
reset all pointers and internal temporary memory
Definition: inverseoperatorinterface.hh:100
Traits ::OperatorType OperatorType
Definition: inverseoperatorinterface.hh:30
Traits ::SolverDiscreteFunctionType SolverDiscreteFunctionType
Definition: inverseoperatorinterface.hh:29
std::shared_ptr< SolverParameterType > parameter_
Definition: inverseoperatorinterface.hh:195
SolverParameterType & parameter() const
Definition: inverseoperatorinterface.hh:127
const AssembledOperatorType * assembledOperator_
Definition: inverseoperatorinterface.hh:198
Traits ::AssembledOperatorType AssembledOperatorType
Definition: inverseoperatorinterface.hh:31
void bind(const OperatorType &op)
store pointer to linear operator
Definition: inverseoperatorinterface.hh:81
int iterations() const
return number of iterations used in previous call of application operator
Definition: inverseoperatorinterface.hh:103
Definition: solver/parameter.hh:15
const std::string & keyPrefix() const
Definition: solver/parameter.hh:66
const ParameterReader & parameter() const
Definition: solver/parameter.hh:68
Definition: cachedcommmanager.hh:47