dune-fem  2.8-git
istlinverseoperators.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
2 #define DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
3 
9 
12 
13 #if HAVE_DUNE_ISTL
14 #include <dune/common/version.hh>
15 
18 
19 #include <dune/istl/scalarproducts.hh>
20 #include <dune/istl/solvers.hh>
21 #include <dune/istl/preconditioner.hh>
22 
23 namespace Dune
24 {
25 
26  namespace Fem
27  {
28 
29  // wrapper for Fem Preconditioners (Operators acting as preconditioners) into ISTL preconditioners
30  template< class Preconditioner >
31  class ISTLPreconditionAdapter
32  : public Dune::Preconditioner< typename Preconditioner::RangeFunctionType::DofStorageType, typename Preconditioner::DomainFunctionType::DofStorageType >
33  {
34  typedef ISTLPreconditionAdapter< Preconditioner > ThisType;
35  typedef Dune::Preconditioner< typename Preconditioner::RangeFunctionType::DofStorageType, typename Preconditioner::DomainFunctionType::DofStorageType > BaseType;
36 
37  typedef typename Preconditioner::DomainFunctionType DomainFunctionType;
38  typedef typename Preconditioner::RangeFunctionType RangeFunctionType;
39 
40  public:
41 #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
42  enum {category=SolverCategory::sequential};
43 #endif // #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
44 
45  typedef typename BaseType::domain_type domain_type;
46  typedef typename BaseType::range_type range_type;
47  typedef typename BaseType::field_type field_type;
48 
49  typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainFunctionSpaceType;
50  typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeFunctionSpaceType;
51 
52  ISTLPreconditionAdapter ( const Preconditioner *precon, const DomainFunctionSpaceType &domainSpace, const RangeFunctionSpaceType &rangeSpace )
53  : precon_( precon ),
54  domainSpace_( domainSpace ),
55  rangeSpace_( rangeSpace )
56  {}
57 
58  // pre and post do nothing here
59  virtual void pre ( domain_type &x, range_type &y ) override {}
60  virtual void post ( domain_type &x ) override {}
61 
62  virtual void apply ( domain_type &x, const range_type &y ) override
63  {
64  // no precon
65  if( !precon_ )
66  {
67  x = y;
68  }
69  else
70  {
71  // note: ISTL switches the arguments !!!
72  // it is assumed that we have a left preconditioner
73  RangeFunctionType px( "ISTLPreconditionAdapter::apply::x", rangeSpace_, x );
74  DomainFunctionType py( "ISTLPreconditionAdapter::apply::y", domainSpace_, y );
75 
76  (*precon_)( px, py );
77  }
78  }
79 
80 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
81  SolverCategory::Category category () const override { return SolverCategory::sequential; }
82 #endif // #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
83 
84  protected:
85  const Preconditioner *precon_;
86  const DomainFunctionSpaceType &domainSpace_;
87  const RangeFunctionSpaceType &rangeSpace_;
88  };
89 
90 
91  template< class BlockVector >
92  struct ISTLSolverReduction
93  {
94  ISTLSolverReduction ( std::shared_ptr< ISTLSolverParameter > parameter )
95  : parameter_( parameter )
96  {}
97 
98  double operator() ( const Dune::LinearOperator< BlockVector, BlockVector > &op,
99  Dune::ScalarProduct< BlockVector > &scp,
100  const BlockVector &rhs, const BlockVector &x ) const
101  {
102 
103  if( parameter_->errorMeasure() == 0)
104  {
105  BlockVector residuum( rhs );
106  op.applyscaleadd( -1., x, residuum );
107  const double res = scp.norm( residuum );
108  return (res > 0 ? parameter_->tolerance() / res : 1e-3);
109  }
110  else
111  return parameter_->tolerance();
112  }
113 
114  private:
115  std::shared_ptr<ISTLSolverParameter> parameter_;
116  };
117 
118  template< int method,
119  class X,
120  class Reduction = ISTLSolverReduction< X > >
121  struct ISTLSolverAdapter
122  {
123  typedef Reduction ReductionType;
124 
125  typedef X domain_type;
126  typedef X range_type;
127 
128  ISTLSolverAdapter ( const ReductionType &reduction, std::shared_ptr<ISTLSolverParameter> parameter )
129  : reduction_( reduction ),
130  method_( method < 0 ?
131  parameter->solverMethod({ SolverParameter::gmres,
138  })
139  : method ),
140  parameter_( parameter )
141  {}
142 
143  template<class Op, class ScP, class PC >
144  void operator () ( Op& op, ScP &scp, PC &pc,
145  range_type &rhs, domain_type &x,
146  Dune::InverseOperatorResult &result ) const
147  {
148  const int verbosity = (Dune::Fem::Parameter::verbose() && parameter_->verbose()) ? 2 : 0;
149  int maxIterations = std::min( std::numeric_limits< int >::max(), parameter_->maxIterations() );
150  if( method_ == SolverParameter::cg )
151  {
152  typedef Dune::CGSolver< X > SolverType;
153  SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
154  solver.apply( x, rhs, result );
155  return ;
156  }
157  else if( method_ == SolverParameter::bicgstab )
158  {
159  typedef Dune::BiCGSTABSolver< X > SolverType;
160  SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
161  solver.apply( x, rhs, result );
162  return ;
163  }
164  else if( method_ == SolverParameter::gmres )
165  {
166  typedef Dune::RestartedGMResSolver< X > SolverType;
167  SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), parameter_->gmresRestart(), maxIterations, verbosity );
168  solver.apply( x, rhs, result );
169  return ;
170  }
171  else if( method_ == SolverParameter::minres )
172  {
173  typedef Dune::MINRESSolver< X > SolverType;
174  SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
175  solver.apply( x, rhs, result );
176  }
177  else if( method_ == SolverParameter::gradient )
178  {
179  typedef Dune::GradientSolver< X > SolverType;
180  SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
181  solver.apply( x, rhs, result );
182  return ;
183  }
184  else if( method_ == SolverParameter::loop )
185  {
186  typedef Dune::LoopSolver< X > SolverType;
187  SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
188  solver.apply( x, rhs, result );
189  return ;
190  }
191  else if( method_ == SolverParameter::superlu )
192  {
193  callSuperLU( op, rhs, x, result );
194  }
195  else
196  {
197  DUNE_THROW(NotImplemented,"ISTLSolverAdapter::operator(): wrong method solver identifier" << method_ );
198  }
199  }
200 
201  void setMaxLinearIterations( unsigned int maxIterations ) { parameter_->setMaxIterations(maxIterations); }
202  void setMaxIterations( unsigned int maxIterations ) { parameter_->setMaxIterations(maxIterations); }
203 
204  std::shared_ptr<ISTLSolverParameter> parameter () const { return parameter_; }
205 
206  protected:
207  template< class ImprovedMatrix >
208  void callSuperLU ( ISTLParallelMatrixAdapterInterface< ImprovedMatrix >& op,
209  range_type &rhs, domain_type &x,
210  Dune::InverseOperatorResult &result ) const
211  {
212 #if HAVE_SUPERLU
213  const int verbosity = (Dune::Fem::Parameter::verbose() && parameter_->verbose()) ? 2 : 0;
214  typedef typename ImprovedMatrix :: BaseType Matrix;
215  const ImprovedMatrix& matrix = op.getmat();
216  SuperLU< Matrix > solver( matrix, verbosity );
217  solver.apply( x, rhs, result );
218 #else
219  DUNE_THROW(NotImplemented,"ISTLSolverAdapter::callSuperLU: SuperLU solver selected but SuperLU not available!");
220 #endif
221  }
222 
223  template< class Op >
224  void callSuperLU ( ISTLLinearOperatorAdapter< Op >& op, range_type &rhs, domain_type &x,
225  Dune::InverseOperatorResult &result ) const
226  {
227  DUNE_THROW(NotImplemented,"ISTLSolverAdapter::callSuperLU: SuperLU only works for AssembledLinearOperators!");
228  }
229 
230  ReductionType reduction_;
231  const int method_;
232 
233  std::shared_ptr<ISTLSolverParameter> parameter_;
234  };
235 
236 
237 
238  // ISTLInverseOperator
239  // -------------------
240 
241  template< class DiscreteFunction, int method = -1,
242  class Preconditioner = const Operator< DiscreteFunction, DiscreteFunction > >
243  class ISTLInverseOperator;
244 
245  template< class DiscreteFunction, int method, class Preconditioner >
246  struct ISTLInverseOperatorTraits
247  {
248  typedef DiscreteFunction DiscreteFunctionType;
250  typedef Preconditioner PreconditionerType;
251 
252  typedef Dune::Fem::ISTLLinearOperator< DiscreteFunction, DiscreteFunction > AssembledOperatorType;
253  typedef ISTLBlockVectorDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType ;
254 
255  typedef ISTLInverseOperator< DiscreteFunction, method, Preconditioner > InverseOperatorType;
256  typedef ISTLSolverParameter SolverParameterType;
257  };
258 
259  template< class DiscreteFunction, int method, class Preconditioner >
260  class ISTLInverseOperator : public InverseOperatorInterface<
261  ISTLInverseOperatorTraits< DiscreteFunction,
262  method, Preconditioner > >
263  {
264  typedef ISTLInverseOperatorTraits< DiscreteFunction, method, Preconditioner > Traits;
265  typedef InverseOperatorInterface< Traits > BaseType;
266 
267  friend class InverseOperatorInterface< Traits >;
268  public:
269  typedef typename BaseType::DomainFunctionType DomainFunctionType;
270  typedef typename BaseType::RangeFunctionType RangeFunctionType;
271  typedef typename BaseType::OperatorType OperatorType;
272  typedef typename BaseType::PreconditionerType PreconditionerType;
273  typedef typename BaseType::AssembledOperatorType AssembledOperatorType;
274 
275  using BaseType :: bind;
276  using BaseType :: unbind;
277  using BaseType :: setMaxIterations;
278  using BaseType :: setMaxLinearIterations;
279 
280  protected:
281  typedef typename DomainFunctionType :: DiscreteFunctionSpaceType
282  DiscreteFunctionSpaceType;
283 
284  typedef ISTLLinearOperatorAdapter< OperatorType > ISTLOperatorType;
285  typedef ISTLPreconditionAdapter< PreconditionerType > ISTLPreconditionerAdapterType;
286 
287  typedef typename RangeFunctionType :: ScalarProductType ParallelScalarProductType;
288  typedef typename RangeFunctionType::DofStorageType BlockVectorType;
289 
290  typedef ISTLSolverAdapter< method, BlockVectorType > SolverAdapterType;
291  typedef typename SolverAdapterType::ReductionType ReductionType;
292  public:
293 
294  //non-deprecated constructors
295  ISTLInverseOperator ( const ISTLSolverParameter & parameter = ISTLSolverParameter() )
296  : BaseType( parameter ), solverAdapter_( ReductionType( parameter_ ), parameter_ )
297  {}
298 
299  ISTLInverseOperator ( const OperatorType &op,
300  const ISTLSolverParameter & parameter = ISTLSolverParameter() )
301  : ISTLInverseOperator ( parameter )
302  {
303  bind( op );
304  }
305 
306  ISTLInverseOperator ( const OperatorType &op, PreconditionerType &preconditioner,
307  const ISTLSolverParameter & parameter = ISTLSolverParameter() )
308  : ISTLInverseOperator( parameter )
309  {
310  bind( op, preconditioner );
311  }
312 
313 
314  protected:
315  // apply for arbitrary domain function type and matching range function type
316  template <class DomainFunction>
317  int apply( const DomainFunction& u, RangeFunctionType& w ) const
318  {
319  auto& scp = w.scalarProduct();
320  // u may not be a discrete function, therefore use w.space()
321  const DiscreteFunctionSpaceType& space = w.space();
322  ISTLPreconditionerAdapterType istlPreconditioner( preconditioner_, space, space );
323 
324  if( assembledOperator_ )
325  {
326  auto& matrix = assembledOperator_->matrixAdapter( *(solverAdapter_.parameter()) );
327  // if preconditioner_ was set use that one, otherwise the one from the matrix object
328  typedef Dune::Preconditioner< BlockVectorType, BlockVectorType > PreconditionerType;
329  PreconditionerType& matrixPre = matrix.preconditionAdapter();
330  PreconditionerType& istlPre = istlPreconditioner;
331  PreconditionerType& precon = ( preconditioner_ ) ? istlPre : matrixPre;
332  return solve( matrix, scp, precon, u, w );
333  }
334  else
335  {
336  assert( operator_ );
337  ISTLOperatorType istlOperator( *operator_, space, space );
338  return solve( istlOperator, scp, istlPreconditioner, u, w );
339  }
340  }
341 
343  template< class OperatorAdapter, class ISTLPreconditioner, class DomainFunction >
344  int solve ( OperatorAdapter &istlOperator, ParallelScalarProductType &scp,
345  ISTLPreconditioner &preconditioner,
346  const DomainFunction& u,
347  RangeFunctionType& w ) const
348  {
349  if( ! rhs_ )
350  {
351  // u may not be a discrete function, therefore use w.space()
352  rhs_.reset( new DomainFunctionType( "ISTLInvOp::rhs", w.space() ) );
353  rightHandSideCopied_ = false;
354  }
355 
356  if( ! rightHandSideCopied_ )
357  {
358  // copy right hand side since ISTL solvers seem to modify it
359  rhs_->assign( u );
360  rightHandSideCopied_ = true;
361  }
362 
363  Dune::InverseOperatorResult result;
364  solverAdapter_.setMaxIterations( parameter_->maxIterations() );
365  solverAdapter_( istlOperator, scp, preconditioner, rhs_->blockVector(), w.blockVector(), result );
366  return (result.converged) ? result.iterations : -(result.iterations);
367  }
368 
369  using BaseType :: operator_;
370  using BaseType :: assembledOperator_;
371  using BaseType :: preconditioner_;
372 
373  using BaseType :: rhs_;
374  using BaseType :: x_;
375 
376  using BaseType :: rightHandSideCopied_;
377  using BaseType :: parameter_;
378 
379  mutable SolverAdapterType solverAdapter_;
380  };
381 
382  } // namespace Fem
383 
384 } // namespace Dune
385 
386 #endif // #if HAVE_ISTL
387 
388 #endif // #ifndef DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
Definition: bindguard.hh:11
static constexpr T max(T a)
Definition: utility.hh:77
static constexpr T min(T a)
Definition: utility.hh:93
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static const int loop
Definition: solver/parameter.hh:29
static const int cg
Definition: solver/parameter.hh:24
static const int gmres
Definition: solver/parameter.hh:26
static const int minres
Definition: solver/parameter.hh:27
static const int gradient
Definition: solver/parameter.hh:28
static const int bicgstab
Definition: solver/parameter.hh:25
static const int superlu
Definition: solver/parameter.hh:30