dune-pdelab  2.5-dev
linearproblem.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
2 #define DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
3 
4 #include <iostream>
5 
6 #include <dune/common/timer.hh>
7 #include <dune/common/parametertree.hh>
8 
12 
13 namespace Dune {
14  namespace PDELab {
15 
16  //===============================================================
17  // A class for solving linear stationary problems.
18  // It assembles the matrix, computes the right hand side and
19  // solves the problem.
20  // This is only a first vanilla implementation which has to be improved.
21  //===============================================================
22 
23  // Status information of linear problem solver
24  template<class RFType>
26  {
27  RFType first_defect; // the first defect
28  RFType defect; // the final defect
29  double assembler_time; // Cumulative time for matrix assembly
30  double linear_solver_time; // Cumulative time for linear solver
31  int linear_solver_iterations; // Total number of linear iterations
32 
34  : first_defect(0.0)
35  , defect(0.0)
36  , assembler_time(0.0)
37  , linear_solver_time(0.0)
38  , linear_solver_iterations(0)
39  {}
40 
41  };
42 
43  template<typename GO, typename LS, typename V>
45  {
46  typedef typename Dune::template FieldTraits<typename V::ElementType >::real_type Real;
47  typedef typename GO::Traits::Jacobian M;
48  typedef typename GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace;
50  typedef GO GridOperator;
51 
52  public:
54 
55  StationaryLinearProblemSolver(const GO& go, LS& ls, V& x, Real reduction, Real min_defect = 1e-99, int verbose=1)
56  : _go(go)
57  , _ls(ls)
58  , _x(&x)
59  , _reduction(reduction)
60  , _min_defect(min_defect)
61  , _hanging_node_modifications(false)
62  , _keep_matrix(true)
63  , _verbose(verbose)
64  {}
65 
66  StationaryLinearProblemSolver (const GO& go, LS& ls, Real reduction, Real min_defect = 1e-99, int verbose=1)
67  : _go(go)
68  , _ls(ls)
69  , _x()
70  , _reduction(reduction)
71  , _min_defect(min_defect)
72  , _hanging_node_modifications(false)
73  , _keep_matrix(true)
74  , _verbose(verbose)
75  {}
76 
78 
95  StationaryLinearProblemSolver(const GO& go, LS& ls, V& x, const ParameterTree& params)
96  : _go(go)
97  , _ls(ls)
98  , _x(&x)
99  , _reduction(params.get<Real>("reduction"))
100  , _min_defect(params.get<Real>("min_defect",1e-99))
101  , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
102  , _keep_matrix(params.get<bool>("keep_matrix",true))
103  , _verbose(params.get<int>("verbosity",1))
104  {}
105 
107 
124  StationaryLinearProblemSolver(const GO& go, LS& ls, const ParameterTree& params)
125  : _go(go)
126  , _ls(ls)
127  , _x()
128  , _reduction(params.get<Real>("reduction"))
129  , _min_defect(params.get<Real>("min_defect",1e-99))
130  , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
131  , _keep_matrix(params.get<bool>("keep_matrix",true))
132  , _verbose(params.get<int>("verbosity",1))
133  {}
134 
137  {
138  _hanging_node_modifications=b;
139  }
140 
143  {
144  return _hanging_node_modifications;
145  }
146 
148  void setKeepMatrix(bool b)
149  {
150  _keep_matrix = b;
151  }
152 
154  bool keepMatrix() const
155  {
156  return _keep_matrix;
157  }
158 
159  const Result& result() const
160  {
161  return _res;
162  }
163 
164  void apply(V& x, bool reuse_matrix = false) {
165  _x = &x;
166  apply(reuse_matrix);
167  }
168 
169  void apply (bool reuse_matrix = false)
170  {
171  Dune::Timer watch;
172  double timing,assembler_time=0;
173 
174  // assemble matrix; optional: assemble only on demand!
175  watch.reset();
176 
177  if (!_jacobian)
178  {
179  _jacobian = std::make_shared<M>(_go);
180  timing = watch.elapsed();
181  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
182  std::cout << "=== matrix setup (max) " << timing << " s" << std::endl;
183  watch.reset();
184  assembler_time += timing;
185  }
186  else if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
187  std::cout << "=== matrix setup skipped (matrix already allocated)" << std::endl;
188 
189  if (_hanging_node_modifications)
190  {
191  Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
192  _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
193  }
194 
195  if (!reuse_matrix)
196  {
197  (*_jacobian) = Real(0.0);
198  _go.jacobian(*_x,*_jacobian);
199  }
200 
201  timing = watch.elapsed();
202  // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
203  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
204  {
205  if (reuse_matrix)
206  std::cout << "=== matrix assembly SKIPPED" << std::endl;
207  else
208  std::cout << "=== matrix assembly (max) " << timing << " s" << std::endl;
209  }
210 
211  assembler_time += timing;
212 
213  // assemble residual
214  watch.reset();
215 
216  W r(_go.testGridFunctionSpace(),0.0);
217  _go.residual(*_x,r); // residual is additive
218 
219  timing = watch.elapsed();
220  // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
221  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
222  std::cout << "=== residual assembly (max) " << timing << " s" << std::endl;
223  assembler_time += timing;
224  _res.assembler_time = assembler_time;
225 
226  auto defect = _ls.norm(r);
227 
228  // compute correction
229  watch.reset();
230  V z(_go.trialGridFunctionSpace(),0.0);
231  auto red = std::max(_reduction,_min_defect/defect);
232  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
233  {
234  std::cout << "=== solving (reduction: " << red << ") ";
235  if (_verbose>=1)
236  std::cout << std::flush;
237  else
238  std::cout << std::endl;
239  }
240  _ls.apply(*_jacobian,z,r,red); // solver makes right hand side consistent
241  _linear_solver_result = _ls.result();
242  timing = watch.elapsed();
243  // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
244  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
245  std::cout << timing << " s" << std::endl;
246  _res.linear_solver_time = timing;
247 
248  _res.converged = _linear_solver_result.converged;
249  _res.iterations = _linear_solver_result.iterations;
250  _res.elapsed = _linear_solver_result.elapsed;
251  _res.reduction = _linear_solver_result.reduction;
252  _res.conv_rate = _linear_solver_result.conv_rate;
253  _res.first_defect = static_cast<double>(defect);
254  _res.defect = static_cast<double>(defect)*_linear_solver_result.reduction;
255  _res.linear_solver_iterations = _linear_solver_result.iterations;
256 
257  // and update
258  if (_hanging_node_modifications)
259  Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
260  *_x -= z;
261  if (_hanging_node_modifications)
262  _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
263 
264  if (!_keep_matrix)
265  _jacobian.reset();
266  }
267 
270  {
271  if(_jacobian)
272  _jacobian.reset();
273  }
274 
276  return _linear_solver_result;
277  }
278 
279  Real reduction() const
280  {
281  return _reduction;
282  }
283 
285  {
286  _reduction = reduction;
287  }
288 
289 
290  private:
291  const GO& _go;
292  LS& _ls;
293  V* _x;
294  shared_ptr<M> _jacobian;
295  Real _reduction;
296  Real _min_defect;
297  Dune::PDELab::LinearSolverResult<double> _linear_solver_result;
298  Result _res;
299  bool _hanging_node_modifications;
300  bool _keep_matrix;
301  int _verbose;
302  };
303 
304  } // namespace PDELab
305 } // namespace Dune
306 
307 #endif // DUNE_PDELAB_STATIONARY_LINEARPROBLEM_HH
Definition: solver.hh:30
RFType reduction
Definition: solver.hh:35
const Entity & e
Definition: localfunctionspace.hh:120
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: linearproblem.hh:269
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: linearproblem.hh:148
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:124
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
bool hangingNodeModifications() const
Return whether the solver performs the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:142
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1016
RFType first_defect
Definition: linearproblem.hh:27
double linear_solver_time
Definition: linearproblem.hh:30
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:55
Definition: linearproblem.hh:44
StationaryLinearProblemSolver(const GO &go, LS &ls, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:66
const Result & result() const
Definition: linearproblem.hh:159
double assembler_time
Definition: linearproblem.hh:29
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
void apply(bool reuse_matrix=false)
Definition: linearproblem.hh:169
const Dune::PDELab::LinearSolverResult< double > & ls_result() const
Definition: linearproblem.hh:275
RFType defect
Definition: linearproblem.hh:28
void setHangingNodeModifications(bool b)
Set whether the solver should apply the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:136
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: linearproblem.hh:154
int linear_solver_iterations
Definition: linearproblem.hh:31
void apply(V &x, bool reuse_matrix=false)
Definition: linearproblem.hh:164
Real reduction() const
Definition: linearproblem.hh:279
void setReduction(Real reduction)
Definition: linearproblem.hh:284
StationaryLinearProblemSolverResult()
Definition: linearproblem.hh:33
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:95
StationaryLinearProblemSolverResult< double > Result
Definition: linearproblem.hh:53