dune-pdelab  2.7-git
blocksorpreconditioner.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_BLOCKSORPRECONDITIONER_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_BLOCKSORPRECONDITIONER_HH
5 
6 namespace Dune {
7  namespace PDELab {
8 
37  template<typename JacobianLOP, typename BlockOffDiagonalLOP, typename GridFunctionSpace>
41  {
42  using value_type = typename GridFunctionSpace::Traits::GridView::Grid::ctype;
44  using W = typename LocalVector::WeightedAccumulationView;
45 
46  public :
47  // A SOR preconditioner includes non-diagonal blocks so we need volume
48  // and skeleton methods. We do two-sided assembly!
49  static constexpr bool doPatternVolume = true;
50  static constexpr bool doPatternSkeleton = true;
51  static constexpr bool doPatternVolumePostSkeleton = true;
52  static constexpr bool doAlphaVolume = true;
53  static constexpr bool doAlphaSkeleton = true;
54  static constexpr bool doAlphaVolumePostSkeleton = true;
55  static constexpr bool doSkeletonTwoSided = true;
56  static constexpr bool isLinear = JacobianLOP::isLinear;
57 
58  BlockSORPreconditionerLocalOperator(JacobianLOP& jacobianlop,
59  BlockOffDiagonalLOP& blockOffDiagonalLOP,
60  const GridFunctionSpace& gridFunctionSpace,
61  const double omega=1.0)
62  : _jacobianlop(jacobianlop)
63  , _blockOffDiagonalLOP(blockOffDiagonalLOP)
64  , _omega(omega)
65  , _a_i(gridFunctionSpace.ordering().maxLocalSize())
66  , _b_i(gridFunctionSpace.ordering().maxLocalSize())
67  {}
68 
73  bool requireSetup()
74  {
75  return _jacobianlop.requireSetup();
76  }
77  void setRequireSetup(bool v)
78  {
79  _jacobianlop.setRequireSetup(v);
80  }
81 
83  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
84  void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
85  {
86  _jacobianlop.alpha_volume(eg,lfsu,x,lfsv,y);
87  }
88 
90  template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
91  void alpha_skeleton(const IG& ig,
92  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
93  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
94  Y& y_s, Y& y_n) const
95  {}
96 
98  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
99  void alpha_volume_post_skeleton(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
100  {}
101 
103  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
104  void jacobian_apply_volume(const EG& eg,
105  const LFSU& lfsu, const X& x,
106  const LFSV& lfsv, Y& y) const
107  {
108  //====================================
109  // How does this preconditioner work
110  //====================================
111  //
112  // When jacobian_apply is called on the grid operator the assembly
113  // machine will iterate over the grid. Since we have set twosided to
114  // true the order of methods will be:
115  //
116  // jacobian_apply_volume
117  // jacobian_apply_skeleton for each intersection
118  // jacobian_apply_volume_post_skeleton
119  //
120  // In the volume part we set _a_i to zero. During the skeleton part we
121  // apply the off-diagonal blocks to the old and new solution, step (1)
122  // of the algorithm description at the top of this class. In the
123  // post_skeleton section we will do steps (2) and (3) of the algorithm.
124  _a_i = 0.0;
125  }
126 
128  template<typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
129  void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const
130  {
131  _a_i = 0.0;
132  }
133 
135  template<typename IG, typename LFSU, typename Z, typename LFSV, typename Y>
136  void jacobian_apply_skeleton(const IG& ig,
137  const LFSU& lfsu_s, const Z& z_s, const LFSV& lfsv_s,
138  const LFSU& lfsu_n, const Z& z_n, const LFSV& lfsv_n,
139  Y& y_s, Y& y_n) const
140  {
141  // This step is bit tricky.
142  //
143  // In the apply method of the GridOperatorPreconditioner class
144  // jacobian_apply(d, v) is called on the preconditioner grid
145  // operator.
146  //
147  // z_s and z_n will be the local entries of the block vector d.
148  //
149  // y_s will always be the solution of the previous step, e.i
150  // $v_i^{(k-1)} in the above algorithm. This will be zero since that is
151  // the value you get from ISTL for v^{(0)} and we do only one
152  // preconditioner step.
153  //
154  // For cells you have not yet visited y_n will also be the old solution
155  // $v_j^{(k-1)}. If you have already visited the neigbor cell it will
156  // instead be the new solution $v_j^{(k)}$.
157  //
158  // For this reason we pass the local vectors of y_s and y_n as
159  // coefficients to the block off-diagonal local operator. The result
160  // will be accumulated in _a_i. Following the above algorithm we will
161  // have
162  //
163  // _a_i = \sum_{j<i} A_ij v_j^{(k)} + \sum_{j>i} A_ij v_j^{(k-1)}
164  //
165  // after visiting all the intersections.
166  //
167  // NOTE: This only work for the FastDGGridOperator!
168  //
169  // This works perfectly for the FastDGGridOperator: The fast-DG
170  // grid operator directly works on the degrees of freedom vector, so if
171  // we set some values during the jacobian apply skeleton we will get
172  // the correct result in the jacobian apply post skeleton.
173  //
174  // For the regular GridOperator we instead write data into a local
175  // degrees of freedom vectors and the results will only written back to
176  // the global data structure during unbinding of the LFS. This means
177  // you get the same old coefficients in the
178  // jacobian_apply_volume_post_skeleton. There are two ways around this issue:
179  //
180  // 1. Copy data, e.g. by storing the values as members of this
181  // class. This does not require changing pdelab, but: Copying data is
182  // slow and the whole point of matrix-free solvers is to avoid
183  // unnecessary memory work. In the end all these solvers are intended
184  // to be used together with the FastDGGridOperator.
185  //
186  // 2. It would be possible to change the GridOperator assembler. This
187  // is not really appealing since it's behavior is perfect for the
188  // regular use case and having an additional unbind/bind would slow
189  // down those regular use cases.
190  //
191  // It is of course possible to achieve both: A fast implementation in
192  // case of FastDGGridOperator and a slower one for regular GridOperator
193  // involving data movement but this would make this whole
194  // implementation even more complicated without big benefits.
195  W _a_i_view(_a_i, y_s.weight());
196  if (ig.inside().partitionType() == Dune::InteriorEntity){
197  // Note: Since the block off diagonal local operator works two sided
198  // we can pass _a_i_view two times to the jacobian_apply_skeleton but
199  // will only accumulate once.
200 
201  // TODO: Only works for FastDGGridOperator (y_* being an AliasedVectorView)
202  _blockOffDiagonalLOP.jacobian_apply_skeleton(ig, lfsu_s, y_s, lfsv_s, lfsu_n, y_n, lfsv_s, _a_i_view, _a_i_view);
203  }
204  }
205 
207  template<typename IG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
208  void jacobian_apply_skeleton(const IG& ig,
209  const LFSU& lfsu_s, const X& x_s, const Z& z_s, const LFSV& lfsv_s,
210  const LFSU& lfsu_n, const X& x_n, const Z& z_n, const LFSV& lfsv_n,
211  Y& y_s, Y& y_n) const
212  {
213  // See documentation above!
214  W _a_i_view(_a_i, y_s.weight());
215  if (ig.inside().partitionType() == Dune::InteriorEntity)
216  _blockOffDiagonalLOP.jacobian_apply_skeleton(ig, lfsu_s, x_s, y_s, lfsv_s, lfsu_n, x_n, y_n, lfsv_s, _a_i_view, _a_i_view);
217  }
218 
220  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
222  const LFSU& lfsu, const X& x,
223  const LFSV& lfsv, Y& y) const
224  {
225  // Subtracting x: _a_i -= x_e. Looking at the algorithm above this
226  // finishes step (1) since the coefficient x is the vector d
227  // above. After this we have r_tmp = - a_i
228  std::transform(_a_i.data(),
229  _a_i.data() + _a_i.size(),
230  x.data(),
231  _a_i.data(),
232  std::minus<value_type>{});
233 
234  // Solve the block diagonal system. This is step (2) of the
235  // algorithm. After this we have _b_i = -b_i.
236  _b_i = 0.0;
237  W _b_i_view(_b_i, y.weight());
238  _jacobianlop.jacobian_apply_volume(eg, lfsu, _a_i, lfsv, _b_i_view);
239 
240 
241  // Do the Update step (3) of the algorithm
242  std::transform(y.data(),
243  y.data() + y.size(),
244  _b_i.data(),
245  y.data(),
246  [=](const double &a,
247  const double& b)
248  {
249  return (1-_omega)*a - _omega*b;
250  });
251  }
252 
254  template<typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
255  void jacobian_apply_volume_post_skeleton(const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const
256  {
257  // static checks
258  //----------------------
259  // static_assert(std::is_same<decltype(x),decltype(_a_i)>::value,"Both types have to be the same for nonlinear Jacobian apply");
260  // static_assert(not std::is_same<decltype(x),decltype(_a_i)>::value,"Both types have to be the same for nonlinear Jacobian apply");
261  //
262  // dynamic crashes
263  //----------------------
264  // x.bar(); // type ConstAliasedVectorView
265  // _a_i.bar(); // type LocalVector<double, >
266 
267  // Subtract x_e: _a_i -= x_e
268  std::transform(_a_i.data(),_a_i.data()+_a_i.size(),
269  z.data(),
270  _a_i.data(),
271  std::minus<value_type>{});
272  // Divide by block diagonal, _b_i = D_{ee}^{-1} _a_i
273  _b_i = 0.0;
274  W _b_i_view(_b_i,y.weight());
275  _jacobianlop.jacobian_apply_volume(eg, lfsu, x, _a_i, lfsv, _b_i_view);
276  // Update vector r_e -= omega*_b_i
277  // <=> r_e = (1-\omega)r_e + \omega D_{ee}^{-1}(x_e - \sum_{e'}A_{e,e'}r_{e'})
278  std::transform(y.data(),
279  y.data()+y.size(),
280  _b_i.data(),
281  y.data(),
282  [=](const double& a,
283  const double& b)
284  {
285  return (1-_omega)*a - _omega*b;
286  });
287  }
288 
289  private :
290 
291  JacobianLOP& _jacobianlop;
292  BlockOffDiagonalLOP& _blockOffDiagonalLOP;
293  const double _omega;
294  mutable LocalVector _a_i;
295  mutable LocalVector _b_i;
296  }; // end class BlockSORPreconditionerLocalOperator
297 
298  } // namespace PDELab
299 } // namespace Dune
300 
301 #endif
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Turn a matrix-free Jacobi-type local preconditioner to a SOR one.
Definition: blocksorpreconditioner.hh:41
void jacobian_apply_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const Z &z_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const Z &z_n, const LFSV &lfsv_n, Y &y_s, Y &y_n) const
Gather off-block-diagonals in Gauss-Seidel process of linearized operator.
Definition: blocksorpreconditioner.hh:208
static constexpr bool doAlphaVolume
Definition: blocksorpreconditioner.hh:52
void setRequireSetup(bool v)
Definition: blocksorpreconditioner.hh:77
static constexpr bool doAlphaVolumePostSkeleton
Definition: blocksorpreconditioner.hh:54
BlockSORPreconditionerLocalOperator(JacobianLOP &jacobianlop, BlockOffDiagonalLOP &blockOffDiagonalLOP, const GridFunctionSpace &gridFunctionSpace, const double omega=1.0)
Definition: blocksorpreconditioner.hh:58
void jacobian_apply_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const Z &z, const LFSV &lfsv, Y &y) const
apply preconditioner after skeleton terms, linearized version
Definition: blocksorpreconditioner.hh:255
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Prepare underlying diagonal block preconditioner.
Definition: blocksorpreconditioner.hh:84
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, Y &y_s, Y &y_n) const
Provide this method, but it actually does not nothing.
Definition: blocksorpreconditioner.hh:91
bool requireSetup()
Definition: blocksorpreconditioner.hh:73
static constexpr bool doSkeletonTwoSided
Definition: blocksorpreconditioner.hh:55
void jacobian_apply_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Apply preconditioner after skeleton terms, linear version.
Definition: blocksorpreconditioner.hh:221
static constexpr bool doPatternVolumePostSkeleton
Definition: blocksorpreconditioner.hh:51
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Linear operator application, volume terms.
Definition: blocksorpreconditioner.hh:104
static constexpr bool doAlphaSkeleton
Definition: blocksorpreconditioner.hh:53
static constexpr bool isLinear
Definition: blocksorpreconditioner.hh:56
static constexpr bool doPatternVolume
Definition: blocksorpreconditioner.hh:49
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const Z &z, const LFSV &lfsv, Y &y) const
linearized operator application, volume terms
Definition: blocksorpreconditioner.hh:129
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Provide this method, but it actually does nothing.
Definition: blocksorpreconditioner.hh:99
static constexpr bool doPatternSkeleton
Definition: blocksorpreconditioner.hh:50
void jacobian_apply_skeleton(const IG &ig, const LFSU &lfsu_s, const Z &z_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const Z &z_n, const LFSV &lfsv_n, Y &y_s, Y &y_n) const
Gather off-block-diagonals in Gauss-Seidel process of linear operator.
Definition: blocksorpreconditioner.hh:136
A grid function space.
Definition: gridfunctionspace.hh:191
size_type size() const
The size of the container.
Definition: localvector.hh:317
WeightedVectorAccumulationView< LocalVector > WeightedAccumulationView
An accumulate-only view of this container that automatically applies a weight to all contributions.
Definition: localvector.hh:211
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