dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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
16namespace 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;
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);
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.
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 weight_type & weight()
Definition: iterativeblockjacobipreconditioner.hh:163
const auto data() const
Definition: iterativeblockjacobipreconditioner.hh:173
WeightedPointDiagonalAccumulationView(Container &container, weight_type weight)
Definition: iterativeblockjacobipreconditioner.hh:184
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
Class for collecting statistics over several invocations.
Definition: solverstatistics.hh:39
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