3#ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ITERATIVEBLOCKJACOBIPRECONDITIONER_HH
4#define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ITERATIVEBLOCKJACOBIPRECONDITIONER_HH
10#include <dune/pdelab/localoperator/pointdiagonalwrapper.hh>
11#include <dune/pdelab/localoperator/blockdiagonalwrapper.hh>
14#include<dune/pdelab/backend/istl/matrixfree/assembledblockjacobipreconditioner.hh>
31 typedef X domain_type;
33 typedef typename X::BaseContainer InvDiagonal;
34 typedef typename X::value_type value_type;
47 LocalPointJacobiPreconditioner(
const InvDiagonal& invDiagonal,
48 const value_type diagonalWeight,
49 const bool precondition=
true)
50 : _precondition(precondition)
51 , _invDiagonal(invDiagonal)
52 , _diagonalWeight(diagonalWeight)
55 void pre (domain_type& x, range_type& b)
override {}
57 void apply (domain_type& v,
const range_type& d)
override
60 std::transform(d.data(),
64 [
this](
const value_type& a,
66 return _diagonalWeight*a*b;
70 std::copy(d.data(),d.data()+d.size(),v.data());
74 void post (domain_type& x)
override {}
77 const bool _precondition;
78 const InvDiagonal& _invDiagonal;
79 const value_type _diagonalWeight;
92 template<
typename BlockDiagonalLocalOperator,
typename W,
typename XView,
typename EG,
typename LFSU,
typename LFSV>
98 using weight_type =
typename W::WeightedAccumulationView::weight_type;
99 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
110 : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
114 , _tmp(lfsu.
size(),0.0)
118 void apply(
const W& z, W& y)
const override
121 typename W::WeightedAccumulationView y_view(y, _weight);
123 applyLocalDiagonalBlock(_blockDiagonalLocalOperator, _eg, _lfsu, z, z, _lfsv, y_view);
125 applyLocalDiagonalBlock(_blockDiagonalLocalOperator, _eg, _lfsu, *_u, z, _lfsv, y_view);
134 void setLinearizationPoint(
const XView& u)
139 void setWeight(
const weight_type& weight)
145 const BlockDiagonalLocalOperator& _blockDiagonalLocalOperator;
155 class WeightedPointDiagonalAccumulationView
160 using value_type =
typename Container::value_type;
161 using size_type =
typename Container::size_type;
162 using weight_type = value_type;
164 const weight_type& weight()
171 return _container.data();
174 const auto data()
const
176 return _container.data();
179 template<
typename LFSV>
180 void accumulate(
const LFSV& lfsv, size_type i, value_type value)
182 _container.data()[lfsv.localIndex(i)] += _weight * value;
185 WeightedPointDiagonalAccumulationView(Container& container, weight_type weight)
186 : _container(container)
191 Container& _container;
213 const size_t maxiter=100,
214 const bool precondition=
true,
215 const double weight=1.0,
262 template<
typename BlockDiagonalLocalOperator,
263 typename PointDiagonalLocalOperator,
265 typename DomainField,
272 using GridView =
typename GridFunctionSpace::Traits::GridViewType;
273 using Grid =
typename GridView::Traits::Grid;
275 using value_type = DomainField;
284 static constexpr bool doPatternVolume =
true;
285 static constexpr bool doAlphaVolume =
true;
286 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
303 const PointDiagonalLocalOperator& pointDiagonalLocalOperator,
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)
316 , _requireSetup(true)
324 return _requireSetup;
326 void setRequireSetup(
bool v)
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
340 const std::size_t
size = lfsu.size();
341 assert(lfsv.size() ==
size);
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);
351 for(
auto& val : _invDiagonalCache[cache_idx]){
358 template<
typename EG,
typename LFSU,
typename Z,
typename LFSV,
typename Y>
361 assert(not _requireSetup);
364 std::size_t cache_idx = _mapper.
index(eg.entity());
365 const auto& invDiagonal = _invDiagonalCache[cache_idx];
369 impl::LocalPointJacobiPreconditioner<LocalVector> pointJacobi(invDiagonal,
373 op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
374 op.setWeight(y.weight());
384 std::copy(z.data(), z.data()+z.size(), z_tmp.
data());
388 solver.
apply(y_tmp, z_tmp, stat);
390 std::transform(y.data(),
394 std::plus<value_type>{});
403 template<
typename EG,
typename LFSU,
typename X,
typename Z,
typename LFSV,
typename Y>
406 assert(not _requireSetup);
409 std::size_t cache_idx = _mapper.
index(eg.entity());
410 const auto& invDiagonal = _invDiagonalCache[cache_idx];
414 impl::LocalPointJacobiPreconditioner<LocalVector> pointJacobi(invDiagonal,
418 op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
419 op.setLinearizationPoint(x);
420 op.setWeight(y.weight());
430 std::copy(z.data(), z.data()+z.size(), z_tmp.
data());
434 solver.
apply(y_tmp, z_tmp, stat);
436 std::transform(y.data(),
440 std::plus<value_type>{});
444 BlockDiagonalLocalOperator _blockDiagonalLocalOperator;
445 PointDiagonalLocalOperator _pointDiagonalLocalOperator;
449 mutable std::vector<InvDiagonal> _invDiagonalCache;
Base class for all implementations of iterative solvers.
Definition: solver.hh:205
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solver.hh:376
A linear operator.
Definition: operators.hh:69
X::field_type field_type
The field type of the operator.
Definition: operators.hh:76
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