dune-pdelab  2.7-git
gridoperatorpreconditioner.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 
4 #ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_GRIDOPERATORPRECONDITIONER_HH
5 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_GRIDOPERATORPRECONDITIONER_HH
6 
7 #include<dune/grid/common/rangegenerators.hh>
8 #include<dune/istl/preconditioner.hh>
9 #include<dune/istl/solvercategory.hh>
10 
11 namespace Dune::PDELab
12 {
13 
18  template<class PrecGO>
20  : public Dune::Preconditioner<typename PrecGO::Traits::Domain,typename PrecGO::Traits::Range>
21  {
22  using Domain = typename PrecGO::Traits::Domain;
23  using Range = typename PrecGO::Traits::Range;
24 
25  public :
26  static constexpr bool isLinear = PrecGO::LocalAssembler::isLinear();
27  Dune::SolverCategory::Category category() const override
28  {
29  return Dune::SolverCategory::sequential;
30  }
31 
32  GridOperatorPreconditioner(const PrecGO& precgo)
33  : _precgo(precgo)
34  , _u(nullptr)
35  {}
36 
39  void setLinearizationPoint(const Domain& u)
40  {
41  _u = &u;
42  }
43 
45  void pre(Domain& v, Range& d) override
46  {
47  if (not isLinear and _u == nullptr)
48  DUNE_THROW(Dune::InvalidStateException, "You seem to apply a preconditioner to a linearized problem but haven't set the linearization point explicitly!");
49 
50  auto& lop = _precgo.localAssembler().localOperator();
51  if (lop.requireSetup()){
52  // We use the residual methods of the preconditioner grid operator to
53  // do some set up work, if required (this could eg be the assembly and
54  // inversion of the block diagonal). This was done through this
55  // interface since the dimension match nicely and it avoids rewriting
56  // assembler code like iterating over the grid, binding local function
57  // spaces, etc.
58  //
59  // In the linear case the Jacobian does not depend on the current
60  // solution (*_u in the nonlinear case). This means we can just pass a
61  // dummy object of the right type, since it won't be used in the local
62  // operator.
63  //
64  // Note:
65  // - We do not have the current solution in the linear case so
66  // passing u is not possible
67  // - Having access to the current solution here (for the linear case) is not
68  // difficult but needs changes in the solver backends and the linear problem
69  // solver.
70  if (isLinear)
71  _precgo.residual(d, d);
72  else
73  _precgo.residual(*_u, d);
74  lop.setRequireSetup(false);
75  }
76  }
77 
78  void apply(Domain& v, const Range& d) override
79  {
80  if(isLinear)
81  _precgo.jacobian_apply(d, v);
82  else
83  _precgo.jacobian_apply(*_u, d, v);
84  }
85 
86  void post(Domain& v) override {}
87 
88  private :
89  const PrecGO& _precgo;
90  const Domain* _u;
91  }; // end class GridOperatorPreconditioner
92 
93 } // namespace Dune::PDELab
94 #endif
Definition: adaptivity.hh:29
Turn a grid operator that represents a preconditioner into an ISTL preconditioner.
Definition: gridoperatorpreconditioner.hh:21
void post(Domain &v) override
Definition: gridoperatorpreconditioner.hh:86
static constexpr bool isLinear
Definition: gridoperatorpreconditioner.hh:26
void pre(Domain &v, Range &d) override
prepare preconditioner
Definition: gridoperatorpreconditioner.hh:45
Dune::SolverCategory::Category category() const override
Definition: gridoperatorpreconditioner.hh:27
GridOperatorPreconditioner(const PrecGO &precgo)
Definition: gridoperatorpreconditioner.hh:32
void apply(Domain &v, const Range &d) override
Definition: gridoperatorpreconditioner.hh:78
void setLinearizationPoint(const Domain &u)
Definition: gridoperatorpreconditioner.hh:39