dune-pdelab  2.0.0
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 
10 #include <dune/common/shared_ptr.hh>
11 #include <dune/common/timer.hh>
12 
14 
15 namespace Dune {
16  namespace PDELab {
17 
19 
29  template<class GOS, class SB, class Coeff>
31  {
32  typedef typename GOS::template MatrixContainer<Coeff>::Type Matrix;
34  <typename GOS::Traits::TrialGridFunctionSpace, Coeff>::Type VectorU;
36  <typename GOS::Traits::TestGridFunctionSpace, Coeff>::Type VectorV;
37 
38  const GOS& gos;
39  SB& sb;
40  shared_ptr<Matrix> m;
41  Coeff reduction;
42  Coeff mindefect;
43 
44  public:
45  StationaryMatrixLinearSolver(const GOS& gos_, SB& sb_, Coeff reduction_,
46  Coeff mindefect_ = 1e-99) :
47  gos(gos_), sb(sb_), reduction(reduction_), mindefect(mindefect_)
48  { }
49 
50  void apply (VectorU& x) {
51  Dune::Timer watch;
52  double timing;
53 
54  if(!m) {
55  // setup new matrix from sparsity pattern
56  watch.reset();
57 
58  m.reset(new Matrix(gos));
59 
60  timing = watch.elapsed();
61  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
62  std::cout << "=== matrix setup " << timing << " s" << std::endl;
63 
64  // assemble matrix
65  watch.reset();
66 
67  *m = 0.0;
68  gos.jacobian(x,*m);
69 
70  timing = watch.elapsed();
71  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
72  std::cout << "=== matrix assembly " << timing << " s" << std::endl;
73  }
74  else {
75  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
76  std::cout << "=== matrix setup skipped" << std::endl
77  << "=== matrix assembly skipped" << std::endl;
78  }
79 
80  // assemble residual
81  watch.reset();
82 
83  VectorV r(gos.testGridFunctionSpace(),0.0);
84  gos.residual(x,r); // residual is additive
85 
86  timing = watch.elapsed();
87  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
88  std::cout << "=== residual assembly " << timing << " s" << std::endl;
89 
90  Coeff defect = sb.norm(r);
91 
92  // compute correction
93  watch.reset();
94  VectorU z(gos.trialGridFunctionSpace(),0.0);
95  Coeff red = std::min(reduction,defect/mindefect);
96  sb.apply(*m,z,r,red); // solver makes right hand side consistent
97  timing = watch.elapsed();
98 
99  if (gos.trialGridFunctionSpace().gridView().comm().rank()==0)
100  std::cout << "=== solving (reduction: " << red << ") "
101  << timing << " s" << std::endl;
102 
103  // and update
104  x -= z;
105  }
106 
107  };
108 
109  } // namespace PDELab
110 } // namespace Dune
111 
112 #endif // DUNE_PDELAB_LINEARPROBLEM_STATIONARYMATRIX_HH
StationaryMatrixLinearSolver(const GOS &gos_, SB &sb_, Coeff reduction_, Coeff mindefect_=1e-99)
Definition: stationarymatrix.hh:45
Definition: backendselector.hh:11
void apply(VectorU &x)
Definition: stationarymatrix.hh:50
A class for solving linear problems with stationary matrices.
Definition: stationarymatrix.hh:30
const E & e
Definition: interpolate.hh:172