DUNE PDELab (git)

iterativeblockjacobipreconditioner.hh
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
9
10#include <dune/pdelab/localoperator/pointdiagonalwrapper.hh>
11#include <dune/pdelab/localoperator/blockdiagonalwrapper.hh>
12
14#include<dune/pdelab/backend/istl/matrixfree/assembledblockjacobipreconditioner.hh>
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 {
39 }
40
47 LocalPointJacobiPreconditioner(const InvDiagonal& invDiagonal,
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 [this](const value_type& a,
65 const value_type& b){
66 return _diagonalWeight*a*b;
67 });
68 }
69 else{
70 std::copy(d.data(),d.data()+d.size(),v.data());
71 }
72 }
73
74 void post (domain_type& x) override {}
75
76 private :
77 const bool _precondition;
78 const InvDiagonal& _invDiagonal;
79 const value_type _diagonalWeight;
80 };
81
82
92 template<typename BlockDiagonalLocalOperator, typename W, typename XView, typename EG, typename LFSU, typename LFSV>
94 : public Dune::LinearOperator<W,W>
95 {
97 using field_type = typename Base::field_type;
98 using weight_type = typename W::WeightedAccumulationView::weight_type;
99 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
100
102 {
104 }
105
106 BlockDiagonalOperator(const BlockDiagonalLocalOperator& blockDiagonalLocalOperator,
107 const EG& eg,
108 const LFSU& lfsu,
109 const LFSV& lfsv)
110 : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
111 , _eg(eg)
112 , _lfsu(lfsu)
113 , _lfsv(lfsv)
114 , _tmp(lfsu.size(),0.0)
115 , _u(nullptr)
116 {}
117
118 void apply(const W& z, W& y) const override
119 {
120 y = 0.0;
121 typename W::WeightedAccumulationView y_view(y, _weight);
122 if (isLinear)
123 applyLocalDiagonalBlock(_blockDiagonalLocalOperator, _eg, _lfsu, z, z, _lfsv, y_view);
124 else
125 applyLocalDiagonalBlock(_blockDiagonalLocalOperator, _eg, _lfsu, *_u, z, _lfsv, y_view);
126 }
127
128 void applyscaleadd(field_type alpha, const W& z, W& y) const override
129 {
130 apply(z, _tmp);
131 y.axpy(alpha, _tmp);
132 }
133
134 void setLinearizationPoint(const XView& u)
135 {
136 _u = &u;
137 }
138
139 void setWeight(const weight_type& weight)
140 {
141 _weight = weight;
142 }
143
144 private :
145 const BlockDiagonalLocalOperator& _blockDiagonalLocalOperator;
146 const EG& _eg;
147 const LFSU& _lfsu;
148 const LFSV& _lfsv;
149 mutable W _tmp;
150 const XView* _u;
151 weight_type _weight;
152 };
153
154 template<typename C>
155 class WeightedPointDiagonalAccumulationView
156 {
157 public :
158 using Container = C;
159
160 using value_type = typename Container::value_type;
161 using size_type = typename Container::size_type;
162 using weight_type = value_type;
163
164 const weight_type& weight()
165 {
166 return _weight;
167 }
168
169 auto data()
170 {
171 return _container.data();
172 }
173
174 auto data() const
175 {
176 return _container.data();
177 }
178
179 template<typename LFSV>
180 void accumulate(const LFSV& lfsv, size_type i, value_type value)
181 {
182 _container.data()[lfsv.localIndex(i)] += _weight * value;
183 }
184
185 WeightedPointDiagonalAccumulationView(Container& container, weight_type weight)
186 : _container(container)
187 , _weight(weight)
188 {}
189
190 private :
191 Container& _container;
192 weight_type _weight;
193 };
194
195 } // namespace impl
196
197
203 {
212 BlockSolverOptions(const double resreduction=1e-5,
213 const size_t maxiter=100,
214 const bool precondition=true,
215 const double weight=1.0,
216 const int verbose=0)
217 : _resreduction(resreduction)
218 , _maxiter(maxiter)
219 , _precondition(precondition)
220 , _weight(weight)
221 , _verbose(verbose)
222 {}
226 size_t _maxiter;
230 double _weight;
233 }; // end struct BlockSolverOptions
234
235
262 template<typename BlockDiagonalLocalOperator,
263 typename PointDiagonalLocalOperator,
264 typename GridFunctionSpace,
265 typename DomainField,
266 template<typename> class IterativeSolver>
270 {
271 // Extract some types
272 using GridView = typename GridFunctionSpace::Traits::GridViewType;
273 using Grid = typename GridView::Traits::Grid;
274 using EntityType = typename Grid::template Codim<0>::Entity;
275 using value_type = DomainField;
277 using InvDiagonal = typename LocalVector::BaseContainer;
278
279 public:
280 // Since this class implements something like D^{-1}v for a diagonal
281 // block matrix D we will only have volume methods. The underlying local
282 // operator that describes the block diagonal might of course have
283 // skeleton or boundary parts.
284 static constexpr bool doPatternVolume = true;
285 static constexpr bool doAlphaVolume = true;
286 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
287
302 IterativeBlockJacobiPreconditionerLocalOperator(const BlockDiagonalLocalOperator& blockDiagonalLocalOperator,
303 const PointDiagonalLocalOperator& pointDiagonalLocalOperator,
304 const GridFunctionSpace& gridFunctionSpace,
305 SolverStatistics<int>& solverStatiscits,
306 BlockSolverOptions solveroptions,
307 const bool verbose=0)
308 : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
309 , _pointDiagonalLocalOperator(pointDiagonalLocalOperator)
310 , _gridFunctionSpace(gridFunctionSpace)
311 , _solverStatistics(solverStatiscits)
312 , _mapper(gridFunctionSpace.gridView())
313 , _invDiagonalCache(_mapper.size())
314 , _solveroptions(solveroptions)
315 , _verbose(verbose)
316 , _requireSetup(true)
317 {}
318
319 // Before we can call the jacobian_apply methods we need to assemble the
320 // point diagonal of the diagonal block. This will be used as a preconditioner
321 // for the iterative matrix free solver on the diagonal block.
322 bool requireSetup()
323 {
324 return _requireSetup;
325 }
326 void setRequireSetup(bool v)
327 {
328 _requireSetup = v;
329 }
330
337 template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
338 void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
339 {
340 const std::size_t size = lfsu.size();
341 assert(lfsv.size() == size);
342
343 // Assemble point diagonal
344 std::size_t cache_idx = _mapper.index(eg.entity());
345 _invDiagonalCache[cache_idx].resize(lfsu.size());
346 using TmpView = impl::WeightedPointDiagonalAccumulationView<InvDiagonal>;
347 TmpView view(_invDiagonalCache[cache_idx], y.weight());
348 assembleLocalPointDiagonal(_pointDiagonalLocalOperator, eg, lfsu, x, lfsv, view);
349
350 // Invert this and store it (will later be used as a preconditioner)
351 for(auto& val : _invDiagonalCache[cache_idx]){
352 val = 1. / val;
353 }
354 }
355
356
358 template<typename EG, typename LFSU, typename Z, typename LFSV, typename Y>
359 void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const Z& z, const LFSV& lfsv, Y& y) const
360 {
361 assert(not _requireSetup);
362
363 // Get correct point diagonal from the cache
364 std::size_t cache_idx = _mapper.index(eg.entity());
365 const auto& invDiagonal = _invDiagonalCache[cache_idx];
366
367 // Create iterative solver with point Jacobi preconditioner for solving
368 // the local block diagonal system
369 impl::LocalPointJacobiPreconditioner<LocalVector> pointJacobi(invDiagonal,
370 _solveroptions._weight,
371 _solveroptions._precondition);
373 op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
374 op.setWeight(y.weight());
376 pointJacobi,
377 _solveroptions._resreduction,
378 _solveroptions._maxiter,
379 _solveroptions._verbose);
381
382 // Set up local right hand side
383 LocalVector z_tmp(lfsu.size(), 0.0);
384 std::copy(z.data(), z.data()+z.size(), z_tmp.data());
385
386 // Solve local blocks iteratively
387 LocalVector y_tmp(lfsv.size(), 0.0);
388 solver.apply(y_tmp, z_tmp, stat);
389 _solverStatistics.append(stat.iterations);
390 std::transform(y.data(),
391 y.data()+y.size(),
392 y_tmp.data(),
393 y.data(),
394 std::plus<value_type>{});
395 } // end jacobian_apply_volume
396
397
403 template<typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
404 void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const
405 {
406 assert(not _requireSetup);
407
408 // Get correct point diagonal from the cache
409 std::size_t cache_idx = _mapper.index(eg.entity());
410 const auto& invDiagonal = _invDiagonalCache[cache_idx];
411
412 // Create iterative solver with point Jacobi preconditioner for solving
413 // the local block diagonal system
414 impl::LocalPointJacobiPreconditioner<LocalVector> pointJacobi(invDiagonal,
415 _solveroptions._weight,
416 _solveroptions._precondition);
418 op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
419 op.setLinearizationPoint(x);
420 op.setWeight(y.weight());
422 pointJacobi,
423 _solveroptions._resreduction,
424 _solveroptions._maxiter,
425 _solveroptions._verbose);
427
428 // Set up local right hand side
429 LocalVector z_tmp(lfsu.size(), 0.0);
430 std::copy(z.data(), z.data()+z.size(), z_tmp.data());
431
432 // Solve local blocks iteratively
433 LocalVector y_tmp(lfsv.size(), 0.0);
434 solver.apply(y_tmp, z_tmp, stat);
435 _solverStatistics.append(stat.iterations);
436 std::transform(y.data(),
437 y.data()+y.size(),
438 y_tmp.data(),
439 y.data(),
440 std::plus<value_type>{});
441 } // end jacobian_apply_volume
442
443 private :
444 BlockDiagonalLocalOperator _blockDiagonalLocalOperator;
445 PointDiagonalLocalOperator _pointDiagonalLocalOperator;
446 const GridFunctionSpace& _gridFunctionSpace;
447 SolverStatistics<int>& _solverStatistics;
449 mutable std::vector<InvDiagonal> _invDiagonalCache;
450 mutable BlockSolverOptions _solveroptions;
451 const int _verbose;
452 bool _requireSetup;
453 }; // end class IterativeBlockJacobiPreconditionerLocalOperator
454
455 } // namespace PDELab
456} // namespace Dune
457#endif
Base class for all implementations of iterative solvers.
Definition: solver.hh:205
void apply(X &x, X &b, double reduction, InverseOperatorResult &res) override
Apply inverse operator with given reduction factor.
Definition: solver.hh:376
A linear operator.
Definition: operators.hh:68
X::field_type field_type
The field type of the operator.
Definition: operators.hh:75
sparsity pattern generator
Definition: pattern.hh:14
A grid function space.
Definition: gridfunctionspace.hh:191
Local operator that can be used to create a fully matrix-free Jacobi preconditioner.
Definition: iterativeblockjacobipreconditioner.hh:270
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:359
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:338
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:404
IterativeBlockJacobiPreconditionerLocalOperator(const BlockDiagonalLocalOperator &blockDiagonalLocalOperator, const PointDiagonalLocalOperator &pointDiagonalLocalOperator, const GridFunctionSpace &gridFunctionSpace, SolverStatistics< int > &solverStatiscits, BlockSolverOptions solveroptions, const bool verbose=0)
Constructor.
Definition: iterativeblockjacobipreconditioner.hh:302
Default flags for all local operators.
Definition: flags.hh:19
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
void append(const T x)
Add new data point.
Definition: solverstatistics.hh:52
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
Index index(const EntityType &e) const
Map entity to array index.
Definition: scsgmapper.hh:71
Implementations of the inverse operator interface.
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Define general, extensible interface for operators. The available implementation wraps a matrix.
Mapper classes are used to attach data to a grid.
Provides a class for collecting statistics on the number of block-solves.
Static tag representing a codimension.
Definition: dimension.hh:24
Statistics about the application of an inverse operator.
Definition: solver.hh:50
int iterations
Number of iterations.
Definition: solver.hh:69
Options for IterativeBlockJacobiPreconditionerLocalOperator.
Definition: iterativeblockjacobipreconditioner.hh:203
double _weight
Weight for point-jacobi.
Definition: iterativeblockjacobipreconditioner.hh:230
size_t _maxiter
Maximal number of iterations.
Definition: iterativeblockjacobipreconditioner.hh:226
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:212
double _resreduction
Residual reduction, i.e. solver accuracy.
Definition: iterativeblockjacobipreconditioner.hh:224
int _verbose
verbosity level
Definition: iterativeblockjacobipreconditioner.hh:232
bool _precondition
Precondition with point-Jacobi?
Definition: iterativeblockjacobipreconditioner.hh:228
Create ISTL operator that applies a local block diagonal.
Definition: iterativeblockjacobipreconditioner.hh:95
Dune::SolverCategory::Category category() const override
Category of the linear operator (see SolverCategory::Category)
Definition: iterativeblockjacobipreconditioner.hh:101
void applyscaleadd(field_type alpha, const W &z, W &y) const override
apply operator to x, scale and add:
Definition: iterativeblockjacobipreconditioner.hh:128
void apply(const W &z, W &y) const override
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
Definition: iterativeblockjacobipreconditioner.hh:118
Category
Definition: solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:25
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)