dune-pdelab  2.4.1
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/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 sovler
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)
233  std::cout << "=== solving (reduction: " << red << ") " << std::endl;
234  _ls.apply(*_jacobian,z,r,red); // solver makes right hand side consistent
235  _linear_solver_result = _ls.result();
236  timing = watch.elapsed();
237  // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
238  if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
239  std::cout << timing << " s" << std::endl;
240  _res.linear_solver_time = timing;
241 
242  _res.converged = _linear_solver_result.converged;
243  _res.iterations = _linear_solver_result.iterations;
244  _res.elapsed = _linear_solver_result.elapsed;
245  _res.reduction = _linear_solver_result.reduction;
246  _res.conv_rate = _linear_solver_result.conv_rate;
247  _res.first_defect = static_cast<double>(defect);
248  _res.defect = static_cast<double>(defect)*_linear_solver_result.reduction;
249  _res.linear_solver_iterations = _linear_solver_result.iterations;
250 
251  // and update
252  if (_hanging_node_modifications)
253  Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
254  *_x -= z;
255  if (_hanging_node_modifications)
256  _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
257 
258  if (!_keep_matrix)
259  _jacobian.reset();
260  }
261 
264  {
265  if(_jacobian)
266  _jacobian.reset();
267  }
268 
270  return _linear_solver_result;
271  }
272 
273  Real reduction() const
274  {
275  return _reduction;
276  }
277 
279  {
280  _reduction = reduction;
281  }
282 
283 
284  private:
285  const GO& _go;
286  LS& _ls;
287  V* _x;
288  shared_ptr<M> _jacobian;
289  Real _reduction;
290  Real _min_defect;
291  Dune::PDELab::LinearSolverResult<double> _linear_solver_result;
292  Result _res;
293  bool _hanging_node_modifications;
294  bool _keep_matrix;
295  int _verbose;
296  };
297 
298  } // namespace PDELab
299 } // namespace Dune
300 
301 #endif
bool hangingNodeModifications() const
Return whether the solver performs the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:142
StationaryLinearProblemSolverResult< double > Result
Definition: linearproblem.hh:53
Real reduction() const
Definition: linearproblem.hh:273
int linear_solver_iterations
Definition: linearproblem.hh:31
StationaryLinearProblemSolverResult()
Definition: linearproblem.hh:33
void setHangingNodeModifications(bool b)
Set whether the solver should apply the necessary transformations for calculations on hanging nodes...
Definition: linearproblem.hh:136
StationaryLinearProblemSolver(const GO &go, LS &ls, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:66
Definition: linearproblem.hh:44
Definition: adaptivity.hh:27
void setReduction(Real reduction)
Definition: linearproblem.hh:278
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: linearproblem.hh:148
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: linearproblem.hh:263
void apply(bool reuse_matrix=false)
Definition: linearproblem.hh:169
const Entity & e
Definition: localfunctionspace.hh:111
RFType defect
Definition: linearproblem.hh:28
double linear_solver_time
Definition: linearproblem.hh:30
RFType reduction
Definition: solver.hh:35
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
RFType first_defect
Definition: linearproblem.hh:27
double assembler_time
Definition: linearproblem.hh:29
Definition: solver.hh:30
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
const Result & result() const
Definition: linearproblem.hh:159
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:113
const Dune::PDELab::LinearSolverResult< double > & ls_result() const
Definition: linearproblem.hh:269
void apply(V &x, bool reuse_matrix=false)
Definition: linearproblem.hh:164
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1016
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: linearproblem.hh:154
StationaryLinearProblemSolver(const GO &go, LS &ls, V &x, Real reduction, Real min_defect=1e-99, int verbose=1)
Definition: linearproblem.hh:55