dune-fem  2.8-git
krylovinverseoperators.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_INVERSEOPERATORS_HH
2 #define DUNE_FEM_SOLVER_INVERSEOPERATORS_HH
3 
8 
10 
15 
17 
18 namespace Dune
19 {
20 
21  namespace Fem
22  {
23 
24  // KrylovInverseOperator
25  // ---------------------
26 
27  template< class DiscreteFunction, int method = -1 >
28  class KrylovInverseOperator;
29 
30  template< class DiscreteFunction, int method >
32  {
33  typedef DiscreteFunction DiscreteFunctionType;
36 
38  typedef DiscreteFunction SolverDiscreteFunctionType ;
39 
42  };
43 
44 
45  template< class DiscreteFunction, int method >
47  : public InverseOperatorInterface< KrylovInverseOperatorTraits< DiscreteFunction, method > >
48  {
51 
52  friend class InverseOperatorInterface< Traits >;
53  public:
59 
60  using BaseType :: bind;
61  using BaseType :: unbind;
64 
65  public:
68  {
69  bind( op );
70  }
71 
72  KrylovInverseOperator ( const OperatorType &op, const PreconditionerType& preconditioner,
74  : KrylovInverseOperator( op, &preconditioner, parameter )
75  {}
76 
79  : BaseType( parameter ),
80  precondObj_(),
81  useDiagonalPreconder_(false),
83  method_( method < 0 ?
84  parameter.solverMethod({ SolverParameter::gmres,
87  : method )
88  {
89  assert( parameter_->errorMeasure() == 0 );
90  if( ! precondObj_ )
91  {
92  bool preconditioning = parameter.parameter().template getValue< bool >( "fem.preconditioning", false );
93  if (!preconditioning)
94  preconditioning = parameter.preconditionMethod(
96  useDiagonalPreconder_ = preconditioning;
97  }
98  }
99 
100  template <class Operator>
101  void bind ( const Operator &op )
102  {
104  {
105  // create diagonal preconditioner
107  preconditioner_ = precondObj_.operator->();
108  }
109  if( precondObj_ )
110  BaseType::bind( op, *precondObj_ );
111  else
112  BaseType::bind( op );
113  }
114 
115  protected:
116  int apply( const DomainFunctionType &u, RangeFunctionType &w ) const
117  {
118  std::ostream* os = nullptr;
119  // only set output when general verbose mode is enabled
120  // (basically to avoid output on every rank)
121  if( verbose_ && Parameter :: verbose() )
122  {
123  os = &std::cout;
124  }
125 
126  int numIter = 0;
127 
129  {
130  if( v_.empty() )
131  {
132  v_.reserve( parameter_->gmresRestart()+1 );
133  for( int i=0; i<=parameter_->gmresRestart(); ++i )
134  {
135  v_.emplace_back( DiscreteFunction( "GMRes::v", u.space() ) );
136  }
137  if( preconditioner_ )
138  v_.emplace_back( DiscreteFunction( "GMRes::z", u.space() ) );
139  }
140 
141  // if solver convergence failed numIter will be negative
143  v_, w, u, parameter_->gmresRestart(),
144  parameter_->tolerance(), parameter_->maxIterations(),
145  parameter_->errorMeasure(), os );
146  }
147  else if( method_ == SolverParameter::bicgstab )
148  {
149  if( v_.empty() )
150  {
151  v_.emplace_back( DomainFunctionType( "BiCGStab::r", u.space() ) );
152  v_.emplace_back( DomainFunctionType( "BiCGStab::r*", u.space() ) );
153  v_.emplace_back( DomainFunctionType( "BiCGStab::p", u.space() ) );
154  v_.emplace_back( DomainFunctionType( "BiCGStab::s", u.space() ) );
155  v_.emplace_back( DomainFunctionType( "BiCGStab::tmp", u.space() ) );
156  if( preconditioner_ )
157  v_.emplace_back( DomainFunctionType( "BiCGStab::z", u.space() ) );
158  }
159 
160  // if solver convergence failed numIter will be negative
162  v_, w, u,
163  parameter_->tolerance(), parameter_->maxIterations(),
164  parameter_->errorMeasure(), os );
165  }
166  else if( method_ == SolverParameter::cg )
167  {
168  if( v_.empty() )
169  {
170  v_.emplace_back( DomainFunctionType( "CG::h", u.space() ) );
171  v_.emplace_back( DomainFunctionType( "CG::r", u.space() ) );
172  v_.emplace_back( DomainFunctionType( "CG::p", u.space() ) );
173 
174  if( preconditioner_ )
175  {
176  v_.emplace_back( DomainFunctionType( "CG::s", u.space() ) );
177  v_.emplace_back( DomainFunctionType( "CG::q", u.space() ) );
178  }
179  }
180 
181  // if solver convergence failed numIter will be negative
183  v_, w, u,
184  parameter_->tolerance(), parameter_->maxIterations(),
185  parameter_->errorMeasure(), os );
186  }
187  else
188  {
189  DUNE_THROW(InvalidStateException,"KrylovInverseOperator: invalid method " << method_ );
190  }
191 
192  return numIter;
193  }
194 
195  template <class LinearOperator>
197  const PreconditionerType *preconditioner,
200  {
201  bind(op);
202  if( ! preconditioner_ )
203  {
204  bool preconditioning = parameter.parameter().template getValue< bool >( "fem.preconditioning", false );
205  if (!preconditioning)
206  preconditioning = parameter.preconditionMethod(
208  if( preconditioning && std::is_base_of< AssembledOperator< DomainFunctionType, DomainFunctionType >, LinearOperator > :: value )
209  {
210  // create diagonal preconditioner
212  preconditioner_ = precondObj_.operator->();
213  }
214  }
215  }
216 
217  protected:
218  using BaseType :: operator_;
220 
223 
224  std::unique_ptr< PreconditionerType > precondObj_;
226 
227  mutable std::vector< DomainFunctionType > v_;
228 
229  const bool verbose_;
230 
231  const int method_;
232  };
233 
234 
235  // CgInverseOperator
236  // -----------------
237 
238  template< class DiscreteFunction >
240 
241 
242  // BicgstabInverseOperator
243  // -----------------------
244 
245  template< class DiscreteFunction >
247 
248 
249  // GmresInverseOperator
250  // --------------------
251 
252  template< class DiscreteFunction >
254 
255 
256  // ParDGGeneralizedMinResInverseOperator
257  // -------------------------------------
258 
259  template< class DiscreteFunction >
261 
262  // ParDGBiCGStabInverseOperator
263  // ----------------------------
264 
265  template< class DiscreteFunction >
267 
268  template <class DiscreteFunctionType, class OpType >
270 
271  template <class DiscreteFunctionType, class OpType >
273 
274  template <class DiscreteFunctionType, class OpType >
276 
277  template <class DiscreteFunctionType, class OpType >
279 
280  template <class DiscreteFunctionType, class OpType >
282 
283  } // namespace Fem
284 
285 } // namespace Dune
286 
287 #endif // DUNE_FEM_SOLVER_INVERSEOPERATORS_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
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static ParameterContainer & container()
Definition: io/parameter.hh:193
abstract affine-linear operator
Definition: operator.hh:87
abstract matrix operator
Definition: operator.hh:124
Precondtioner, multiplies with inverse of the diagonal works with.
Definition: diagonalpreconditioner.hh:152
Definition: inverseoperatorinterface.hh:17
virtual void setMaxIterations(const int iter)
Definition: inverseoperatorinterface.hh:113
int iterations_
Definition: inverseoperatorinterface.hh:205
const PreconditionerType * preconditioner_
Definition: inverseoperatorinterface.hh:199
Traits ::PreconditionerType PreconditionerType
Definition: inverseoperatorinterface.hh:32
const OperatorType * operator_
Definition: inverseoperatorinterface.hh:197
virtual void setMaxLinearIterations(const int iter)
set number of max linear iterations to be used before an exception is thrown
Definition: inverseoperatorinterface.hh:108
void unbind()
reset all pointers and internal temporary memory
Definition: inverseoperatorinterface.hh:100
Traits ::OperatorType OperatorType
Definition: inverseoperatorinterface.hh:30
BaseType ::DomainFunctionType DomainFunctionType
Definition: inverseoperatorinterface.hh:26
BaseType ::RangeFunctionType RangeFunctionType
Definition: inverseoperatorinterface.hh:27
std::shared_ptr< SolverParameterType > parameter_
Definition: inverseoperatorinterface.hh:195
SolverParameterType & parameter() const
Definition: inverseoperatorinterface.hh:127
Traits ::AssembledOperatorType AssembledOperatorType
Definition: inverseoperatorinterface.hh:31
void bind(const OperatorType &op)
store pointer to linear operator
Definition: inverseoperatorinterface.hh:81
Definition: krylovinverseoperators.hh:48
int apply(const DomainFunctionType &u, RangeFunctionType &w) const
Definition: krylovinverseoperators.hh:116
BaseType::AssembledOperatorType AssembledOperatorType
Definition: krylovinverseoperators.hh:58
BaseType::RangeFunctionType RangeFunctionType
Definition: krylovinverseoperators.hh:55
void bind(const Operator &op)
Definition: krylovinverseoperators.hh:101
BaseType::DomainFunctionType DomainFunctionType
Definition: krylovinverseoperators.hh:54
BaseType::PreconditionerType PreconditionerType
Definition: krylovinverseoperators.hh:57
BaseType::OperatorType OperatorType
Definition: krylovinverseoperators.hh:56
const int method_
Definition: krylovinverseoperators.hh:231
bool useDiagonalPreconder_
Definition: krylovinverseoperators.hh:225
KrylovInverseOperator(const SolverParameter &parameter=SolverParameter(Parameter::container()))
main constructor
Definition: krylovinverseoperators.hh:78
KrylovInverseOperator(const LinearOperator &op, const PreconditionerType *preconditioner, const SolverParameter &parameter=SolverParameter(Parameter::container()))
Definition: krylovinverseoperators.hh:196
std::unique_ptr< PreconditionerType > precondObj_
Definition: krylovinverseoperators.hh:224
KrylovInverseOperator(const OperatorType &op, const SolverParameter &parameter=SolverParameter(Parameter::container()))
Definition: krylovinverseoperators.hh:66
std::vector< DomainFunctionType > v_
Definition: krylovinverseoperators.hh:227
const bool verbose_
Definition: krylovinverseoperators.hh:229
KrylovInverseOperator(const OperatorType &op, const PreconditionerType &preconditioner, const SolverParameter &parameter=SolverParameter(Parameter::container()))
Definition: krylovinverseoperators.hh:72
Definition: krylovinverseoperators.hh:32
KrylovInverseOperator< DiscreteFunction, method > InverseOperatorType
Definition: krylovinverseoperators.hh:40
DiscreteFunction DiscreteFunctionType
Definition: krylovinverseoperators.hh:33
OperatorType PreconditionerType
Definition: krylovinverseoperators.hh:35
SolverParameter SolverParameterType
Definition: krylovinverseoperators.hh:41
OperatorType AssembledOperatorType
Definition: krylovinverseoperators.hh:37
Operator< DiscreteFunction, DiscreteFunction > OperatorType
Definition: krylovinverseoperators.hh:34
DiscreteFunction SolverDiscreteFunctionType
Definition: krylovinverseoperators.hh:38
Definition: solver/parameter.hh:15
static const int none
Definition: solver/parameter.hh:40
static const int cg
Definition: solver/parameter.hh:24
static const int gmres
Definition: solver/parameter.hh:26
static const int bicgstab
Definition: solver/parameter.hh:25
static const int jacobi
Definition: solver/parameter.hh:45