dune-pdelab  2.7-git
backends.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_BACKENDS_HH
5 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_BACKENDS_HH
6 
7 #include<dune/grid/common/rangegenerators.hh>
8 #include<dune/istl/preconditioner.hh>
9 #include<dune/istl/operators.hh>
10 #include<dune/istl/solvers.hh>
11 
16 
17 
18 namespace Dune {
19  namespace PDELab {
20 
21  namespace impl{
22 
23  // Can be used to check if a local operator is a
24  // BlockSORPreconditionerLocalOperator
25  template <typename>
26  struct isBlockSORPreconditionerLocalOperator : std::false_type {};
27 
28  template <typename JacobianLOP, typename BlockOffDiagonalLOP, typename GridFunctionSpace>
30  BlockOffDiagonalLOP,
32  : std::true_type {};
33 
34  // Can be used to check if a grid operator is a FastDGGridOperator
35  template <typename>
36  struct isFastDGGridOperator : std::false_type {};
37 
38  template<typename GFSU, typename GFSV, typename LOP,
39  typename MB, typename DF, typename RF, typename JF,
40  typename CU, typename CV>
41  struct isFastDGGridOperator<FastDGGridOperator<GFSU, GFSV, LOP, MB, DF, RF, JF, CU, CV >>
42  : std::true_type {};
43 
44  }
45 
62  template<class GO, class PrecGO,
63  template<class> class Solver>
67  {
68  using V = typename GO::Traits::Domain;
69  using W = typename GO::Traits::Range;
72 
73  public :
74  ISTLBackend_SEQ_MatrixFree_Base(const GO& go, const PrecGO& precgo,
75  unsigned maxiter=5000, int verbose=1)
76  : opa_(go), seqprec_(precgo)
77  , maxiter_(maxiter), verbose_(verbose)
78  {
79  // First lets brake that down: The case we want to avoid is having SOR
80  // with regular grid operator. We check for SOR and not fast grid
81  // operator and assert that this is not the case.
82  //
83  // For more information why this is not working have a look at the
84  // documentation of the class BlockSORPreconditionerLocalOperator
86  typename PrecGO::Traits::LocalAssembler::LocalOperator>::value
88  "If you use the BlockSORPreconditioner you need to use FastDGGridOperator!");
89 
90 
91  // We need to assert that the local operator uses the new interface and not the old one
92  static_assert(!decltype(Dune::PDELab::impl::hasOldLOPInterface(go.localAssembler().localOperator()))::value,
93  "\n"
94  "In order to use the matrix-free solvers your grid operator must follow the new\n"
95  "local operator interface for nonlinear jacobian applys, where the current\n"
96  "solution and the linearization point can have different types. This is best\n"
97  "explained by an example.\n\n"
98  "old: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const X& z, const LFSV& lfsv, Y& y) const {}\n"
99  "new: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const {}\n\n"
100  "Note that in the new interface the type of z is Z and doesn't have to be X.\n"
101  "You need to adjust all the nonlinear jacobian apply methods in your local operator\n"
102  "in a similar way.");
103  static_assert(!decltype(Dune::PDELab::impl::hasOldLOPInterface(precgo.localAssembler().localOperator()))::value,
104  "\n"
105  "In order to use the matrix-free solvers your grid operator must follow the new\n"
106  "local operator interface for nonlinear jacobian applys, where the current\n"
107  "solution and the linearization point can have different types. This is best\n"
108  "explained by an example.\n\n"
109  "old: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const X& z, const LFSV& lfsv, Y& y) const {}\n"
110  "new: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const {}\n\n"
111  "Note that in the new interface the type of z is Z and doesn't have to be X.\n"
112  "You need to adjust all the nonlinear jacobian apply methods in your local operator\n"
113  "in a similar way.");
114  }
115 
116  void apply(V& z, W& r, typename V::ElementType reduction)
117  {
118  Solver<V> solver(opa_,seqprec_,reduction,maxiter_,verbose_);
119  Dune::InverseOperatorResult stat;
120  solver.apply(z, r, stat);
121  res.converged = stat.converged;
122  res.iterations = stat.iterations;
123  res.elapsed = stat.elapsed;
124  res.reduction = stat.reduction;
125  res.conv_rate = stat.conv_rate;
126  }
127 
130  void setLinearizationPoint(const V& u)
131  {
132  opa_.setLinearizationPoint(u);
133  seqprec_.setLinearizationPoint(u);
134  }
135 
136  constexpr static bool isMatrixFree {true};
137 
138  private :
139  Operator opa_;
140  SeqPrec seqprec_;
141  unsigned maxiter_;
142  int verbose_;
143  }; // end class ISTLBackend_SEQ_MatrixFree_Base
144 
145 
146  // A matrix based SOR backend for comparing the matrix-free version
148  : public ISTLBackend_SEQ_Base<Dune::SeqSOR, Dune::BiCGSTABSolver>
149  {
150  public:
156  explicit ISTLBackend_SEQ_BCGS_SOR (unsigned maxiter_=5000, int verbose_=1)
157  : ISTLBackend_SEQ_Base<Dune::SeqSOR, Dune::BiCGSTABSolver>(maxiter_, verbose_)
158  {}
159  };
160 
161  } // namespace PDELab
162 } // namespace Dune
163 #endif
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
constexpr auto hasOldLOPInterface(T &t) -> typename std::enable_if<(decltype(hasOldOrNewJacobianApplyVolume(t))::value &&!decltype(hasNewJacobianApplyVolume(t))::value)||(decltype(hasOldOrNewJacobianApplyVolumePostSkeleton(t))::value &&!decltype(hasNewJacobianApplyVolumePostSkeleton(t))::value)||(decltype(hasOldOrNewJacobianApplySkeleton(t))::value &&!decltype(hasNewJacobianApplySkeleton(t))::value)||(decltype(hasOldOrNewJacobianApplyBoundary(t))::value &&!decltype(hasNewJacobianApplyBoundary(t))::value), std::true_type >::type
Definition: checklopinterface.hh:153
Sequential matrix-free solver backend.
Definition: backends.hh:67
constexpr static bool isMatrixFree
Definition: backends.hh:136
ISTLBackend_SEQ_MatrixFree_Base(const GO &go, const PrecGO &precgo, unsigned maxiter=5000, int verbose=1)
Definition: backends.hh:74
void setLinearizationPoint(const V &u)
Definition: backends.hh:130
void apply(V &z, W &r, typename V::ElementType reduction)
Definition: backends.hh:116
Definition: backends.hh:149
ISTLBackend_SEQ_BCGS_SOR(unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: backends.hh:156
Turn a matrix-free Jacobi-type local preconditioner to a SOR one.
Definition: blocksorpreconditioner.hh:41
Turn a grid operator that represents a preconditioner into an ISTL preconditioner.
Definition: gridoperatorpreconditioner.hh:21
void setLinearizationPoint(const Domain &u)
Definition: gridoperatorpreconditioner.hh:39
void setLinearizationPoint(const X &u)
Definition: seqistlsolverbackend.hh:61
Definition: seqistlsolverbackend.hh:211
Definition: solver.hh:17
RFType conv_rate
Definition: solver.hh:36
bool converged
Definition: solver.hh:32
RFType reduction
Definition: solver.hh:35
unsigned int iterations
Definition: solver.hh:33
double elapsed
Definition: solver.hh:34
Definition: solver.hh:54
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:63
A grid function space.
Definition: gridfunctionspace.hh:191
Definition: fastdg.hh:34
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139