3#ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ITERATIVEBLOCKJACOBIPRECONDITIONER_HH
4#define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ITERATIVEBLOCKJACOBIPRECONDITIONER_HH
6#include<dune/grid/common/scsgmapper.hh>
7#include<dune/istl/solvers.hh>
8#include<dune/istl/operators.hh>
36 Dune::SolverCategory::Category
category()
const override
38 return Dune::SolverCategory::sequential;
49 const bool precondition=
true)
50 : _precondition(precondition)
51 , _invDiagonal(invDiagonal)
52 , _diagonalWeight(diagonalWeight)
60 std::transform(d.data(),
65 return _diagonalWeight*a*b;
69 std::copy(d.data(),d.data()+d.size(),v.data());
76 const bool _precondition;
91 template<
typename BlockDiagonalLocalOperator,
typename W,
typename XView,
typename EG,
typename LFSU,
typename LFSV>
93 :
public Dune::LinearOperator<W,W>
95 using Base = Dune::LinearOperator<W,W>;
97 using weight_type =
typename W::WeightedAccumulationView::weight_type;
98 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
100 Dune::SolverCategory::Category
category()
const override
102 return Dune::SolverCategory::sequential;
109 : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
113 , _tmp(lfsu.size(),0.0)
117 void apply(
const W& z, W& y)
const override
120 typename W::WeightedAccumulationView y_view(y, _weight);
144 const BlockDiagonalLocalOperator& _blockDiagonalLocalOperator;
170 return _container.data();
175 return _container.data();
178 template<
typename LFSV>
181 _container.data()[lfsv.localIndex(i)] += _weight *
value;
185 : _container(container)
212 const size_t maxiter=100,
213 const bool precondition=
true,
214 const double weight=1.0,
261 template<
typename BlockDiagonalLocalOperator,
262 typename PointDiagonalLocalOperator,
264 typename DomainField,
265 template<
typename>
class IterativeSolver>
272 using Grid =
typename GridView::Traits::Grid;
273 using EntityType =
typename Grid::template Codim<0>::Entity;
274 using value_type = DomainField;
285 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
302 const PointDiagonalLocalOperator& pointDiagonalLocalOperator,
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)
315 , _requireSetup(true)
323 return _requireSetup;
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
339 const int size = lfsu.size();
340 assert(lfsv.size() == size);
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());
350 for(
auto& val : _invDiagonalCache[cache_idx]){
357 template<
typename EG,
typename LFSU,
typename Z,
typename LFSV,
typename Y>
360 assert(not _requireSetup);
363 std::size_t cache_idx = _mapper.index(eg.entity());
364 const auto& invDiagonal = _invDiagonalCache[cache_idx];
372 op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
374 IterativeSolver<LocalVector> solver(op,
379 Dune::InverseOperatorResult stat;
383 std::copy(z.data(), z.data()+z.size(), z_tmp.
data());
387 solver.apply(y_tmp, z_tmp, stat);
388 _solverStatistics.
append(stat.iterations);
389 std::transform(y.data(),
393 std::plus<value_type>{});
402 template<
typename EG,
typename LFSU,
typename X,
typename Z,
typename LFSV,
typename Y>
405 assert(not _requireSetup);
408 std::size_t cache_idx = _mapper.index(eg.entity());
409 const auto& invDiagonal = _invDiagonalCache[cache_idx];
417 op(_blockDiagonalLocalOperator, eg, lfsu, lfsv);
420 IterativeSolver<LocalVector> solver(op,
425 Dune::InverseOperatorResult stat;
429 std::copy(z.data(), z.data()+z.size(), z_tmp.
data());
433 solver.apply(y_tmp, z_tmp, stat);
434 _solverStatistics.
append(stat.iterations);
435 std::transform(y.data(),
439 std::plus<value_type>{});
443 BlockDiagonalLocalOperator _blockDiagonalLocalOperator;
444 PointDiagonalLocalOperator _pointDiagonalLocalOperator;
447 typename Dune::SingleCodimSingleGeomTypeMapper<GridView, 0> _mapper;
448 mutable std::vector<InvDiagonal> _invDiagonalCache;
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