dune-pdelab  2.7-git
iterativeblockjacobipreconditioner.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 #ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ITERATIVEBLOCKJACOBIPRECONDITIONER_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ITERATIVEBLOCKJACOBIPRECONDITIONER_HH
5 
6 #include<dune/grid/common/scsgmapper.hh>
7 #include<dune/istl/solvers.hh>
8 #include<dune/istl/operators.hh>
9 
12 
15 
16 namespace Dune {
17  namespace PDELab {
18 
19  namespace impl{
20 
21  /* \brief Point Jacobi preconditioner for local diagonal block systems
22  *
23  * This class will do the Jacobi preconditioning for a diagonal
24  * block. The values of the inverse of the diagonal need to be provided
25  * during construction.
26  */
27  template<typename X>
28  class LocalPointJacobiPreconditioner : public Dune::Preconditioner<X,X>
29  {
30  public :
31  typedef X domain_type;
32  typedef X range_type;
33  typedef typename X::BaseContainer InvDiagonal;
34  typedef typename X::value_type value_type;
35 
36  Dune::SolverCategory::Category category() const override
37  {
38  return Dune::SolverCategory::sequential;
39  }
40 
48  const value_type diagonalWeight,
49  const bool precondition=true)
50  : _precondition(precondition)
51  , _invDiagonal(invDiagonal)
52  , _diagonalWeight(diagonalWeight)
53  {}
54 
55  void pre (domain_type& x, range_type& b) override {}
56 
57  void apply (domain_type& v, const range_type& d) override
58  {
59  if(_precondition){
60  std::transform(d.data(),
61  d.data()+d.size(),
62  _invDiagonal.begin(),
63  v.data(),
64  [=](const value_type& a, const value_type& b){
65  return _diagonalWeight*a*b;
66  });
67  }
68  else{
69  std::copy(d.data(),d.data()+d.size(),v.data());
70  }
71  }
72 
73  void post (domain_type& x) override {}
74 
75  private :
76  const bool _precondition;
77  const InvDiagonal& _invDiagonal;
78  const value_type _diagonalWeight;
79  };
80 
81 
91  template<typename BlockDiagonalLocalOperator, typename W, typename XView, typename EG, typename LFSU, typename LFSV>
93  : public Dune::LinearOperator<W,W>
94  {
95  using Base = Dune::LinearOperator<W,W>;
96  using field_type = typename Base::field_type;
97  using weight_type = typename W::WeightedAccumulationView::weight_type;
98  static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
99 
100  Dune::SolverCategory::Category category() const override
101  {
102  return Dune::SolverCategory::sequential;
103  }
104 
105  BlockDiagonalOperator(const BlockDiagonalLocalOperator& blockDiagonalLocalOperator,
106  const EG& eg,
107  const LFSU& lfsu,
108  const LFSV& lfsv)
109  : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
110  , _eg(eg)
111  , _lfsu(lfsu)
112  , _lfsv(lfsv)
113  , _tmp(lfsu.size(),0.0)
114  , _u(nullptr)
115  {}
116 
117  void apply(const W& z, W& y) const override
118  {
119  y = 0.0;
120  typename W::WeightedAccumulationView y_view(y, _weight);
121  if (isLinear)
122  applyLocalDiagonalBlock(_blockDiagonalLocalOperator, _eg, _lfsu, z, z, _lfsv, y_view);
123  else
124  applyLocalDiagonalBlock(_blockDiagonalLocalOperator, _eg, _lfsu, *_u, z, _lfsv, y_view);
125  }
126 
127  void applyscaleadd(field_type alpha, const W& z, W& y) const override
128  {
129  apply(z, _tmp);
130  y.axpy(alpha, _tmp);
131  }
132 
133  void setLinearizationPoint(const XView& u)
134  {
135  _u = &u;
136  }
137 
138  void setWeight(const weight_type& weight)
139  {
140  _weight = weight;
141  }
142 
143  private :
144  const BlockDiagonalLocalOperator& _blockDiagonalLocalOperator;
145  const EG& _eg;
146  const LFSU& _lfsu;
147  const LFSV& _lfsv;
148  mutable W _tmp;
149  const XView* _u;
150  weight_type _weight;
151  };
152 
153  template<typename C>
155  {
156  public :
157  using Container = C;
158 
159  using value_type = typename Container::value_type;
160  using size_type = typename Container::size_type;
162 
164  {
165  return _weight;
166  }
167 
168  auto data()
169  {
170  return _container.data();
171  }
172 
173  const auto data() const
174  {
175  return _container.data();
176  }
177 
178  template<typename LFSV>
179  void accumulate(const LFSV& lfsv, size_type i, value_type value)
180  {
181  _container.data()[lfsv.localIndex(i)] += _weight * value;
182  }
183 
185  : _container(container)
186  , _weight(weight)
187  {}
188 
189  private :
190  Container& _container;
191  weight_type _weight;
192  };
193 
194  } // namespace impl
195 
196 
202  {
211  BlockSolverOptions(const double resreduction=1e-5,
212  const size_t maxiter=100,
213  const bool precondition=true,
214  const double weight=1.0,
215  const int verbose=0)
216  : _resreduction(resreduction)
217  , _maxiter(maxiter)
218  , _precondition(precondition)
219  , _weight(weight)
220  , _verbose(verbose)
221  {}
225  size_t _maxiter;
229  double _weight;
231  int _verbose;
232  }; // end struct BlockSolverOptions
233 
234 
261  template<typename BlockDiagonalLocalOperator,
262  typename PointDiagonalLocalOperator,
263  typename GridFunctionSpace,
264  typename DomainField,
265  template<typename> class IterativeSolver>
269  {
270  // Extract some types
271  using GridView = typename GridFunctionSpace::Traits::GridViewType;
272  using Grid = typename GridView::Traits::Grid;
273  using EntityType = typename Grid::template Codim<0>::Entity;
274  using value_type = DomainField;
276  using InvDiagonal = typename LocalVector::BaseContainer;
277 
278  public:
279  // Since this class implements something like D^{-1}v for a diagonal
280  // block matrix D we will only have volume methods. The underlying local
281  // operator that describes the block diagonal might of course have
282  // skeleton or boundary parts.
283  static constexpr bool doPatternVolume = true;
284  static constexpr bool doAlphaVolume = true;
285  static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
286 
301  IterativeBlockJacobiPreconditionerLocalOperator(const BlockDiagonalLocalOperator& blockDiagonalLocalOperator,
302  const PointDiagonalLocalOperator& pointDiagonalLocalOperator,
303  const GridFunctionSpace& gridFunctionSpace,
304  SolverStatistics<int>& solverStatiscits,
305  BlockSolverOptions solveroptions,
306  const bool verbose=0)
307  : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
308  , _pointDiagonalLocalOperator(pointDiagonalLocalOperator)
309  , _gridFunctionSpace(gridFunctionSpace)
310  , _solverStatistics(solverStatiscits)
311  , _mapper(gridFunctionSpace.gridView())
312  , _invDiagonalCache(_mapper.size())
313  , _solveroptions(solveroptions)
314  , _verbose(verbose)
315  , _requireSetup(true)
316  {}
317 
318  // Before we can call the jacobian_apply methods we need to assemble the
319  // point diagonal of the diagonal block. This will be used as a preconditioner
320  // for the iterative matrix free solver on the diagonal block.
322  {
323  return _requireSetup;
324  }
325  void setRequireSetup(bool v)
326  {
327  _requireSetup = v;
328  }
329 
336  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
337  void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
338  {
339  const int size = lfsu.size();
340  assert(lfsv.size() == size);
341 
342  // Assemble point diagonal
343  std::size_t cache_idx = _mapper.index(eg.entity());
344  _invDiagonalCache[cache_idx].resize(lfsu.size());
346  TmpView view(_invDiagonalCache[cache_idx], y.weight());
347  assembleLocalPointDiagonal(_pointDiagonalLocalOperator, eg, lfsu, x, lfsv, view);
348 
349  // Invert this and store it (will later be used as a preconditioner)
350  for(auto& val : _invDiagonalCache[cache_idx]){
351  val = 1. / val;
352  }
353  }
354 
355 
357  template<typename EG, typename LFSU, typename Z, typename LFSV, typename Y>
358  void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const Z& z, const LFSV& lfsv, Y& y) const
359  {
360  assert(not _requireSetup);
361 
362  // Get correct point diagonal from the cache
363  std::size_t cache_idx = _mapper.index(eg.entity());
364  const auto& invDiagonal = _invDiagonalCache[cache_idx];
365 
366  // Create iterative solver with point Jacobi preconditioner for solving
367  // the local block diagonal system
369  _solveroptions._weight,
370  _solveroptions._precondition);
372  op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
373  op.setWeight(y.weight());
374  IterativeSolver<LocalVector> solver(op,
375  pointJacobi,
376  _solveroptions._resreduction,
377  _solveroptions._maxiter,
378  _solveroptions._verbose);
379  Dune::InverseOperatorResult stat;
380 
381  // Set up local right hand side
382  LocalVector z_tmp(lfsu.size(), 0.0);
383  std::copy(z.data(), z.data()+z.size(), z_tmp.data());
384 
385  // Solve local blocks iteratively
386  LocalVector y_tmp(lfsv.size(), 0.0);
387  solver.apply(y_tmp, z_tmp, stat);
388  _solverStatistics.append(stat.iterations);
389  std::transform(y.data(),
390  y.data()+y.size(),
391  y_tmp.data(),
392  y.data(),
393  std::plus<value_type>{});
394  } // end jacobian_apply_volume
395 
396 
402  template<typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
403  void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const
404  {
405  assert(not _requireSetup);
406 
407  // Get correct point diagonal from the cache
408  std::size_t cache_idx = _mapper.index(eg.entity());
409  const auto& invDiagonal = _invDiagonalCache[cache_idx];
410 
411  // Create iterative solver with point Jacobi preconditioner for solving
412  // the local block diagonal system
414  _solveroptions._weight,
415  _solveroptions._precondition);
417  op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
418  op.setLinearizationPoint(x);
419  op.setWeight(y.weight());
420  IterativeSolver<LocalVector> solver(op,
421  pointJacobi,
422  _solveroptions._resreduction,
423  _solveroptions._maxiter,
424  _solveroptions._verbose);
425  Dune::InverseOperatorResult stat;
426 
427  // Set up local right hand side
428  LocalVector z_tmp(lfsu.size(), 0.0);
429  std::copy(z.data(), z.data()+z.size(), z_tmp.data());
430 
431  // Solve local blocks iteratively
432  LocalVector y_tmp(lfsv.size(), 0.0);
433  solver.apply(y_tmp, z_tmp, stat);
434  _solverStatistics.append(stat.iterations);
435  std::transform(y.data(),
436  y.data()+y.size(),
437  y_tmp.data(),
438  y.data(),
439  std::plus<value_type>{});
440  } // end jacobian_apply_volume
441 
442  private :
443  BlockDiagonalLocalOperator _blockDiagonalLocalOperator;
444  PointDiagonalLocalOperator _pointDiagonalLocalOperator;
445  const GridFunctionSpace& _gridFunctionSpace;
446  SolverStatistics<int>& _solverStatistics;
447  typename Dune::SingleCodimSingleGeomTypeMapper<GridView, 0> _mapper;
448  mutable std::vector<InvDiagonal> _invDiagonalCache;
449  mutable BlockSolverOptions _solveroptions;
450  const int _verbose;
451  bool _requireSetup;
452  }; // end class IterativeBlockJacobiPreconditionerLocalOperator
453 
454  } // namespace PDELab
455 } // namespace Dune
456 #endif
Provides a class for collecting statistics on the number of block-solves.
const Entity & e
Definition: localfunctionspace.hh:123
GridView GridViewType
Definition: gridfunctionspace.hh:135
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
void applyLocalDiagonalBlock(const LocalOperator &localOperator, const EG &eg, const LFSU &lfsu, const X &x, const Z &z, const LFSV &lfsv, Y &y)
A function for applying a single diagonal block.
Definition: blockdiagonalwrapper.hh:261
void assembleLocalPointDiagonal(const LocalOperator &localOperator, const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y)
A function for assembling the point diagonal of a single block.
Definition: pointdiagonalwrapper.hh:139
Definition: iterativeblockjacobipreconditioner.hh:29
void apply(domain_type &v, const range_type &d) override
Definition: iterativeblockjacobipreconditioner.hh:57
X domain_type
Definition: iterativeblockjacobipreconditioner.hh:31
LocalPointJacobiPreconditioner(const InvDiagonal &invDiagonal, const value_type diagonalWeight, const bool precondition=true)
Constructor.
Definition: iterativeblockjacobipreconditioner.hh:47
Dune::SolverCategory::Category category() const override
Definition: iterativeblockjacobipreconditioner.hh:36
void pre(domain_type &x, range_type &b) override
Definition: iterativeblockjacobipreconditioner.hh:55
void post(domain_type &x) override
Definition: iterativeblockjacobipreconditioner.hh:73
X range_type
Definition: iterativeblockjacobipreconditioner.hh:32
X::BaseContainer InvDiagonal
Definition: iterativeblockjacobipreconditioner.hh:33
X::value_type value_type
Definition: iterativeblockjacobipreconditioner.hh:34
Create ISTL operator that applies a local block diagonal.
Definition: iterativeblockjacobipreconditioner.hh:94
Dune::SolverCategory::Category category() const override
Definition: iterativeblockjacobipreconditioner.hh:100
typename W::WeightedAccumulationView::weight_type weight_type
Definition: iterativeblockjacobipreconditioner.hh:97
typename Base::field_type field_type
Definition: iterativeblockjacobipreconditioner.hh:96
void applyscaleadd(field_type alpha, const W &z, W &y) const override
Definition: iterativeblockjacobipreconditioner.hh:127
static constexpr bool isLinear
Definition: iterativeblockjacobipreconditioner.hh:98
void setLinearizationPoint(const XView &u)
Definition: iterativeblockjacobipreconditioner.hh:133
Dune::LinearOperator< W, W > Base
Definition: iterativeblockjacobipreconditioner.hh:95
BlockDiagonalOperator(const BlockDiagonalLocalOperator &blockDiagonalLocalOperator, const EG &eg, const LFSU &lfsu, const LFSV &lfsv)
Definition: iterativeblockjacobipreconditioner.hh:105
void apply(const W &z, W &y) const override
Definition: iterativeblockjacobipreconditioner.hh:117
void setWeight(const weight_type &weight)
Definition: iterativeblockjacobipreconditioner.hh:138
Definition: iterativeblockjacobipreconditioner.hh:155
auto data()
Definition: iterativeblockjacobipreconditioner.hh:168
typename Container::value_type value_type
Definition: iterativeblockjacobipreconditioner.hh:159
C Container
Definition: iterativeblockjacobipreconditioner.hh:157
const auto data() const
Definition: iterativeblockjacobipreconditioner.hh:173
WeightedPointDiagonalAccumulationView(Container &container, weight_type weight)
Definition: iterativeblockjacobipreconditioner.hh:184
const weight_type & weight()
Definition: iterativeblockjacobipreconditioner.hh:163
value_type weight_type
Definition: iterativeblockjacobipreconditioner.hh:161
typename Container::size_type size_type
Definition: iterativeblockjacobipreconditioner.hh:160
void accumulate(const LFSV &lfsv, size_type i, value_type value)
Definition: iterativeblockjacobipreconditioner.hh:179
Options for IterativeBlockJacobiPreconditionerLocalOperator.
Definition: iterativeblockjacobipreconditioner.hh:202
double _weight
Weight for point-jacobi.
Definition: iterativeblockjacobipreconditioner.hh:229
size_t _maxiter
Maximal number of iterations.
Definition: iterativeblockjacobipreconditioner.hh:225
BlockSolverOptions(const double resreduction=1e-5, const size_t maxiter=100, const bool precondition=true, const double weight=1.0, const int verbose=0)
Constructor.
Definition: iterativeblockjacobipreconditioner.hh:211
double _resreduction
Residual reduction, i.e. solver accuracy.
Definition: iterativeblockjacobipreconditioner.hh:223
int _verbose
verbosity level
Definition: iterativeblockjacobipreconditioner.hh:231
bool _precondition
Precondition with point-Jacobi?
Definition: iterativeblockjacobipreconditioner.hh:227
Local operator that can be used to create a fully matrix-free Jacobi preconditioner.
Definition: iterativeblockjacobipreconditioner.hh:269
static constexpr bool doPatternVolume
Definition: iterativeblockjacobipreconditioner.hh:283
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const Z &z, const LFSV &lfsv, Y &y) const
Apply fully matrix-free preconditioner - linear case.
Definition: iterativeblockjacobipreconditioner.hh:358
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Prepare fully matrix-free preconditioner.
Definition: iterativeblockjacobipreconditioner.hh:337
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const Z &z, const LFSV &lfsv, Y &y) const
Apply fully matrix-free preconditioner - nonlinear case.
Definition: iterativeblockjacobipreconditioner.hh:403
IterativeBlockJacobiPreconditionerLocalOperator(const BlockDiagonalLocalOperator &blockDiagonalLocalOperator, const PointDiagonalLocalOperator &pointDiagonalLocalOperator, const GridFunctionSpace &gridFunctionSpace, SolverStatistics< int > &solverStatiscits, BlockSolverOptions solveroptions, const bool verbose=0)
Constructor.
Definition: iterativeblockjacobipreconditioner.hh:301
void setRequireSetup(bool v)
Definition: iterativeblockjacobipreconditioner.hh:325
static constexpr bool doAlphaVolume
Definition: iterativeblockjacobipreconditioner.hh:284
static constexpr bool isLinear
Definition: iterativeblockjacobipreconditioner.hh:285
bool requireSetup()
Definition: iterativeblockjacobipreconditioner.hh:321
void append(const T x)
Add new data point.
Definition: solverstatistics.hh:52
A grid function space.
Definition: gridfunctionspace.hh:191
std::vector< value_type > BaseContainer
The type of the underlying storage container.
Definition: localvector.hh:188
auto data()
Access underlying container.
Definition: localvector.hh:244
Default flags for all local operators.
Definition: flags.hh:19
sparsity pattern generator
Definition: pattern.hh:14
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139