dune-pdelab  2.0.0
linearproblem.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_STATIONARYLINEARPROBLEM_HH
2 #define DUNE_PDELAB_STATIONARYLINEARPROBLEM_HH
3 
4 #include <iostream>
5 
6 #include <dune/common/timer.hh>
7 #include <dune/common/deprecated.hh>
8 #include <dune/common/parametertree.hh>
9 
13 
14 namespace Dune {
15  namespace PDELab {
16 
17  //===============================================================
18  // A class for solving linear stationary problems.
19  // It assembles the matrix, computes the right hand side and
20  // solves the problem.
21  // This is only a first vanilla implementation which has to be improved.
22  //===============================================================
23 
24  // Status information of linear problem solver
25  template<class RFType>
27  {
28  RFType first_defect; // the first defect
29  RFType defect; // the final defect
30  double assembler_time; // Cumulative time for matrix assembly
31  double linear_solver_time; // Cumulative time for linear sovler
32  int linear_solver_iterations; // Total number of linear iterations
33 
35  : first_defect(0.0)
36  , defect(0.0)
37  , assembler_time(0.0)
38  , linear_solver_time(0.0)
40  {}
41 
42  };
43 
44  template<typename GO, typename LS, typename V>
46  {
47  typedef typename V::ElementType Real;
48  typedef typename GO::Traits::Jacobian M;
49  typedef typename GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace;
51  typedef GO GridOperator;
52 
53  public:
55 
56  StationaryLinearProblemSolver(const GO& go, V& x, LS& ls, typename V::ElementType reduction, typename V::ElementType min_defect = 1e-99, int verbose=1) DUNE_DEPRECATED_MSG("Use StationaryLinearProblemSolver(const GO&, LS&, V&, ...) instead.")
57  : _go(go)
58  , _ls(ls)
59  , _x(&x)
60  , _reduction(reduction)
61  , _min_defect(min_defect)
62  , _hanging_node_modifications(false)
63  , _keep_matrix(true)
64  , _verbose(verbose)
65  {}
66 
67  StationaryLinearProblemSolver(const GO& go, LS& ls, V& x, typename V::ElementType reduction, typename V::ElementType min_defect = 1e-99, int verbose=1)
68  : _go(go)
69  , _ls(ls)
70  , _x(&x)
71  , _reduction(reduction)
72  , _min_defect(min_defect)
73  , _hanging_node_modifications(false)
74  , _keep_matrix(true)
75  , _verbose(verbose)
76  {}
77 
78  StationaryLinearProblemSolver (const GO& go, LS& ls, typename V::ElementType reduction, typename V::ElementType min_defect = 1e-99, int verbose=1)
79  : _go(go)
80  , _ls(ls)
81  , _x()
82  , _reduction(reduction)
83  , _min_defect(min_defect)
84  , _hanging_node_modifications(false)
85  , _keep_matrix(true)
86  , _verbose(verbose)
87  {}
88 
90 
107  StationaryLinearProblemSolver(const GO& go, LS& ls, V& x, const ParameterTree& params)
108  : _go(go)
109  , _ls(ls)
110  , _x(&x)
111  , _reduction(params.get<typename V::ElementType>("reduction"))
112  , _min_defect(params.get<typename V::ElementType>("min_defect",1e-99))
113  , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
114  , _keep_matrix(params.get<bool>("keep_matrix",true))
115  , _verbose(params.get<int>("verbosity",1))
116  {}
117 
119 
136  StationaryLinearProblemSolver(const GO& go, LS& ls, const ParameterTree& params)
137  : _go(go)
138  , _ls(ls)
139  , _x()
140  , _reduction(params.get<typename V::ElementType>("reduction"))
141  , _min_defect(params.get<typename V::ElementType>("min_defect",1e-99))
142  , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
143  , _keep_matrix(params.get<bool>("keep_matrix",true))
144  , _verbose(params.get<int>("verbosity",1))
145  {}
146 
149  {
150  _hanging_node_modifications=b;
151  }
152 
155  {
156  return _hanging_node_modifications;
157  }
158 
160  void setKeepMatrix(bool b)
161  {
162  _keep_matrix = b;
163  }
164 
166  bool keepMatrix() const
167  {
168  return _keep_matrix;
169  }
170 
171  const Result& result() const
172  {
173  return _res;
174  }
175 
176  void apply(V& x) {
177  _x = &x;
178  apply();
179  }
180 
181  void apply ()
182  {
183  Dune::Timer watch;
184  double timing,assembler_time=0;
185 
186  // assemble matrix; optional: assemble only on demand!
187  watch.reset();
188 
189  if (!_jacobian)
190  {
191  _jacobian = make_shared<M>(_go);
192  timing = watch.elapsed();
193  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
194  std::cout << "=== matrix setup (max) " << timing << " s" << std::endl;
195  watch.reset();
196  assembler_time += timing;
197  }
198  else if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
199  std::cout << "=== matrix setup skipped (matrix already allocated)" << std::endl;
200 
201  (*_jacobian) = Real(0.0);
202  if (_hanging_node_modifications)
203  {
204  Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
205  _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
206  }
207  _go.jacobian(*_x,*_jacobian);
208 
209  timing = watch.elapsed();
210  // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
211  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
212  std::cout << "=== matrix assembly (max) " << timing << " s" << std::endl;
213  assembler_time += timing;
214 
215  // assemble residual
216  watch.reset();
217 
218  W r(_go.testGridFunctionSpace(),0.0);
219  _go.residual(*_x,r); // residual is additive
220 
221  timing = watch.elapsed();
222  // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
223  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
224  std::cout << "=== residual assembly (max) " << timing << " s" << std::endl;
225  assembler_time += timing;
226  _res.assembler_time = assembler_time;
227 
228  typename V::ElementType defect = _ls.norm(r);
229 
230  // compute correction
231  watch.reset();
232  V z(_go.trialGridFunctionSpace(),0.0);
233  typename V::ElementType red = std::max(_reduction,_min_defect/defect);
234  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0)
235  std::cout << "=== solving (reduction: " << red << ") ";
236  _ls.apply(*_jacobian,z,r,red); // solver makes right hand side consistent
237  _linear_solver_result = _ls.result();
238  timing = watch.elapsed();
239  // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
240  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
241  std::cout << timing << " s" << std::endl;
242  _res.linear_solver_time = timing;
243 
244  _res.converged = _linear_solver_result.converged;
245  _res.iterations = _linear_solver_result.iterations;
246  _res.elapsed = _linear_solver_result.elapsed;
247  _res.reduction = _linear_solver_result.reduction;
248  _res.conv_rate = _linear_solver_result.conv_rate;
249  _res.first_defect = static_cast<double>(defect);
250  _res.defect = static_cast<double>(defect)*_linear_solver_result.reduction;
251  _res.linear_solver_iterations = _linear_solver_result.iterations;
252 
253  // and update
254  if (_hanging_node_modifications)
255  Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
256  *_x -= z;
257  if (_hanging_node_modifications)
258  _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
259 
260  if (!_keep_matrix)
261  _jacobian.reset();
262  }
263 
266  {
267  if(_jacobian)
268  _jacobian.reset();
269  }
270 
272  return _linear_solver_result;
273  }
274 
275  private:
276  const GO& _go;
277  LS& _ls;
278  V* _x;
279  shared_ptr<M> _jacobian;
280  typename V::ElementType _reduction;
281  typename V::ElementType _min_defect;
282  Dune::PDELab::LinearSolverResult<double> _linear_solver_result;
283  Result _res;
284  bool _hanging_node_modifications;
285  bool _keep_matrix;
286  int _verbose;
287  };
288 
289  } // namespace PDELab
290 } // namespace Dune
291 
292 #endif
StationaryLinearProblemSolver(const GO &go, V &x, LS &ls, typename V::ElementType reduction, typename V::ElementType min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:56
bool converged
Definition: solver.hh:32
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: linearproblem.hh:265
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: common/constraints.hh:1009
double elapsed
Definition: solver.hh:34
void setHangingNodeModifications(bool b)
Set whether the solver should apply the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:148
StationaryLinearProblemSolverResult< double > Result
Definition: linearproblem.hh:54
Definition: solver.hh:30
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: linearproblem.hh:160
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: linearproblem.hh:166
const Dune::PDELab::LinearSolverResult< double > & ls_result() const
Definition: linearproblem.hh:271
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, typename V::ElementType reduction, typename V::ElementType min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:67
bool hangingNodeModifications() const
Return whether the solver performs the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:154
RFType reduction
Definition: solver.hh:35
void apply(V &x)
Definition: linearproblem.hh:176
BackendVectorSelectorHelper< Backend, GridFunctionSpace, FieldType >::Type Type
Definition: backendselector.hh:14
double linear_solver_time
Definition: linearproblem.hh:31
const Result & result() const
Definition: linearproblem.hh:171
Definition: linearproblem.hh:45
StationaryLinearProblemSolverResult()
Definition: linearproblem.hh:34
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, const ParameterTree &params)
Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterT...
Definition: linearproblem.hh:107
StationaryLinearProblemSolver(const GO &go, LS &ls, typename V::ElementType reduction, typename V::ElementType min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:78
const E & e
Definition: interpolate.hh:172
StationaryLinearProblemSolver(const GO &go, LS &ls, const ParameterTree &params)
Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterT...
Definition: linearproblem.hh:136
RFType conv_rate
Definition: solver.hh:36
RFType defect
Definition: linearproblem.hh:29
void apply()
Definition: linearproblem.hh:181
RFType first_defect
Definition: linearproblem.hh:28
int linear_solver_iterations
Definition: linearproblem.hh:32
double assembler_time
Definition: linearproblem.hh:30
unsigned int iterations
Definition: solver.hh:33