dune-pdelab  2.4-dev
stationarymatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_LINEARPROBLEM_STATIONARYMATRIX_HH
5 #define DUNE_PDELAB_LINEARPROBLEM_STATIONARYMATRIX_HH
6 
7 #include <iostream>
8 #include <ostream>
9 #include <memory>
10 
11 #include <dune/common/timer.hh>
12 #include <dune/common/deprecated.hh>
13 
14 #warning <dune/pdelab/linearsolver/stationarymatrix.hh> and StationaryMatrixLinearSolver are deprecated and will be removed after PDELab 2.4. Please use LinearProblemSolver instead.
15 
17 
18 namespace Dune {
19  namespace PDELab {
20 
22 
32  template<class GOS, class SB, class Coeff>
33  class
34  DUNE_DEPRECATED_MSG("StationaryMatrixLinearSolver is deprecated and will be removed after PDELab 2.4. Please use LinearProblemSolver instead.")
36  {
37  typedef typename GOS::template MatrixContainer<Coeff>::Type Matrix;
38  using VectorU = Dune::PDELab::Backend::Vector
39  <typename GOS::Traits::TrialGridFunctionSpace, Coeff>;
40  using VectorV = Dune::PDELab::Backend::Vector
41  <typename GOS::Traits::TestGridFunctionSpace, Coeff>;
42 
43  const GOS& gos;
44  SB& sb;
45  std::shared_ptr<Matrix> m;
46  Coeff reduction;
47  Coeff mindefect;
48 
49  public:
50  StationaryMatrixLinearSolver(const GOS& gos_, SB& sb_, Coeff reduction_,
51  Coeff mindefect_ = 1e-99) :
52  gos(gos_), sb(sb_), reduction(reduction_), mindefect(mindefect_)
53  { }
54 
55  void apply (VectorU& x) {
56  Dune::Timer watch;
57  double timing;
58 
59  if(!m) {
60  // setup new matrix from sparsity pattern
61  watch.reset();
62 
63  m.reset(new Matrix(gos));
64 
65  timing = watch.elapsed();
66  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
67  std::cout << "=== matrix setup " << timing << " s" << std::endl;
68 
69  // assemble matrix
70  watch.reset();
71 
72  *m = 0.0;
73  gos.jacobian(x,*m);
74 
75  timing = watch.elapsed();
76  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
77  std::cout << "=== matrix assembly " << timing << " s" << std::endl;
78  }
79  else {
80  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
81  std::cout << "=== matrix setup skipped" << std::endl
82  << "=== matrix assembly skipped" << std::endl;
83  }
84 
85  // assemble residual
86  watch.reset();
87 
88  VectorV r(gos.testGridFunctionSpace(),0.0);
89  gos.residual(x,r); // residual is additive
90 
91  timing = watch.elapsed();
92  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
93  std::cout << "=== residual assembly " << timing << " s" << std::endl;
94 
95  Coeff defect = sb.norm(r);
96 
97  // compute correction
98  watch.reset();
99  VectorU z(gos.trialGridFunctionSpace(),0.0);
100  Coeff red = std::min(reduction,defect/mindefect);
101  sb.apply(*m,z,r,red); // solver makes right hand side consistent
102  timing = watch.elapsed();
103 
104  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
105  std::cout << "=== solving (reduction: " << red << ") "
106  << timing << " s" << std::endl;
107 
108  // and update
109  x -= z;
110  }
111 
112  };
113 
114  } // namespace PDELab
115 } // namespace Dune
116 
117 #endif // DUNE_PDELAB_LINEARPROBLEM_STATIONARYMATRIX_HH
StationaryMatrixLinearSolver(const GOS &gos_, SB &sb_, Coeff reduction_, Coeff mindefect_=1e-99)
Definition: stationarymatrix.hh:50
void apply(VectorU &x)
Definition: stationarymatrix.hh:55
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:112
const E & e
Definition: interpolate.hh:172
Definition: adaptivity.hh:27
typename impl::BackendMatrixSelector< Backend, VU, VV, E >::Type Matrix
alias of the return type of BackendMatrixSelector
Definition: backend/interface.hh:133
A class for solving linear problems with stationary matrices.
Definition: stationarymatrix.hh:33