DUNE PDELab (git)

assembledblockjacobipreconditioner.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_ASSEMBLEDBLOCKJACOBIPREONDITIONER_HH
4#define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ASSEMBLEDBLOCKJACOBIPREONDITIONER_HH
5
6#if HAVE_EIGEN
7
9
10#include <dune/pdelab/localoperator/blockdiagonalwrapper.hh>
11
12#include<Eigen/Dense>
13
14
15namespace Dune {
16 namespace PDELab {
17
18 namespace impl {
19 // Below we will assemble a block of the matrix as Eigen matrix. This
20 // accumulation view makes sure that we apply the weights.
21 template<typename C>
22 class WeightedEigenMatrixAccumulationView
23 {
24 public :
25 typedef C Container;
26
27 typedef typename Container::Scalar value_type;
28 using size_type = typename Container::StorageIndex;
29 using weight_type = value_type;
30
31 const weight_type& weight()
32 {
33 return _weight;
34 }
35
36 auto data()
37 {
38 return _container.data();
39 }
40
41 auto data() const
42 {
43 return _container.data();
44 }
45
46 auto& container()
47 {
48 return _container;
49 }
50
51 template <typename LFSV, typename LFSU>
52 void accumulate(const LFSV& lfsv, size_type i, const LFSU& lfsu, size_type j, value_type value)
53 {
54 _container(i,j) += _weight * value;
55 }
56
57 WeightedEigenMatrixAccumulationView(Container& container, weight_type weight)
58 : _container(container)
59 , _weight(weight)
60 {}
61
62 private :
63 Container& _container;
64 weight_type _weight;
65 };
66 }
67
68
96 template<class BlockDiagonalLocalOperator, class GridFunctionSpace, class DomainField>
97 class AssembledBlockJacobiPreconditionerLocalOperator
100 {
101 // Extract some types
102 using GridView = typename GridFunctionSpace::Traits::GridViewType;
103 using Grid = typename GridView::Traits::Grid;
104 using EntityType = typename Grid::template Codim<0>::Entity;
105 using value_type = DomainField;
106 using DiagonalBlock = ::Eigen::Matrix<value_type, ::Eigen::Dynamic, ::Eigen::Dynamic, ::Eigen::RowMajor>;
107 using PermutationIndices = ::Eigen::Matrix<int, ::Eigen::Dynamic, 1>;
108 using TupleType = std::tuple<DiagonalBlock, PermutationIndices>;
109
110 public :
111
112 // Since this class implements something like D^{-1}v for a diagonal
113 // block matrix D we will only have volume methods. The underlying local
114 // operator that describes the block diagonal might of course have
115 // skeleton or boundary parts.
116 static constexpr bool doPatternVolume = true;
117 static constexpr bool doAlphaVolume = true;
118 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
119
130 AssembledBlockJacobiPreconditionerLocalOperator(const BlockDiagonalLocalOperator& blockDiagonalLocalOperator,
131 const GridFunctionSpace& gridFunctionSpace,
132 const bool verbose=0)
133 : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
134 , _gridFunctionSpace(gridFunctionSpace)
135 , _mapper(gridFunctionSpace.gridView())
136 , _precCache(_mapper.size())
137 , _verbose(verbose)
138 , _requireSetup(true)
139 {
140 }
141
142
143 // Before we can call the jacobian_apply methods we need to compute the
144 // LU decomposition of the diagonal block. We use the bool _requireSetup
145 // to make sure that we have computed this decomposition before we try to
146 // use it.
147 bool requireSetup()
148 {
149 return _requireSetup;
150 }
151 void setRequireSetup(bool requireSetup)
152 {
153 _requireSetup = requireSetup;
154 }
155
156
162 template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
163 void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
164 {
165 // Matrix for storing the LU decomposition of the block diagonal
166 const auto size = lfsu.size();
167 assert(lfsv.size() == size);
168
169 // Assemble local block diagonal
170 DiagonalBlock diagonalBlock = DiagonalBlock::Constant(size, size, 0.0);
171 impl::WeightedEigenMatrixAccumulationView<DiagonalBlock> viewDiagonalBlock(diagonalBlock, y.weight());
172 assembleLocalDiagonalBlock(_blockDiagonalLocalOperator, eg, lfsu, x, lfsv, viewDiagonalBlock);
173
174 // Compute the LU-factorization directly inside the matrix
175 ::Eigen::PartialPivLU<::Eigen::Ref<DiagonalBlock>> lu(diagonalBlock);
176 auto inversePermutationIndices = lu.permutationP().indices();
177
178 // We need the inverse of the permutation matrix that Eigen gives
179 PermutationIndices permutationIndices(size);
180 for(int i=0; i<inversePermutationIndices.size(); ++i)
181 permutationIndices[inversePermutationIndices[i]] = i;
182
183 // Fill up correct entry of preconditioner cache
184 std::size_t cache_idx = _mapper.index(eg.entity());
185 _precCache[cache_idx] = std::make_tuple(diagonalBlock, permutationIndices);
186 }
187
188 template<typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
189 void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const
190 {
191 // We don't need the linearization point as the preconditioner is
192 // already set up
193 jacobian_apply_volume(eg,lfsu,z,lfsv,y);
194 }
195
196
202 template<typename EG, typename LFSU, typename Z, typename LFSV, typename Y>
203 void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const Z& z, const LFSV& lfsv, Y& y) const
204 {
205 assert(not _requireSetup);
206
207 // Get correct LU decomposition of the block diagonal from the cache
208 std::size_t cache_idx = _mapper.index(eg.entity());
209 const auto& cacheEntry = _precCache[cache_idx];
210 auto diagonalBlockLU = std::get<0>(cacheEntry);
211 auto permutationIndices = std::get<1>(cacheEntry);
212
213 // Apply the preconditioner
214 const int size = lfsu.size();
215 const auto& z_container = accessBaseContainer(z);
216 std::vector<value_type> ztmp_container(size);
217 auto& y_container = accessBaseContainer(y);
218
219 // Forward solve with lower triangle
220 for(int i=0; i<size; ++i){
221 value_type rhs(z_container[permutationIndices[i]]);
222 for(int j=0; j<i; ++j)
223 rhs -= diagonalBlockLU(i, j) * ztmp_container[j];
224 ztmp_container[i] = rhs; // L has ones on the diagonal
225 }
226 // Backward solve with upper triangular
227 for(int i=size-1; i>=0; i--){
228 value_type rhs(ztmp_container[i]);
229 for(int j=size-1; j>i; j--)
230 rhs -= diagonalBlockLU(i, j) * y_container[j];
231 y_container[i] = rhs/diagonalBlockLU(i,i);
232 }
233 }
234
235 private :
236 BlockDiagonalLocalOperator _blockDiagonalLocalOperator;
237 const GridFunctionSpace& _gridFunctionSpace;
239 mutable std::vector<TupleType> _precCache;
240 const int _verbose;
241 bool _requireSetup;
242 };
243
244 } // namespace PDELab
245} // namespace Dune
246
247#endif // HAVE_EIGEN
248
249#endif
sparsity pattern generator
Definition: pattern.hh:14
Default flags for all local operators.
Definition: flags.hh:19
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
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
Mapper classes are used to attach data to a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)