dune-pdelab  2.7-git
solver/newton.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_PDELAB_SOLVER_NEWTON_HH
4 #define DUNE_PDELAB_SOLVER_NEWTON_HH
5 
6 #include <dune/common/exceptions.hh>
7 #include <dune/common/ios_state.hh>
8 
13 
14 namespace Dune::PDELab
15 {
16  namespace Impl
17  {
18  // Some SFinae magic to execute setReuse(bool) on a backend
19  template<typename T1, typename = void>
20  struct HasSetReuse
21  : std::false_type
22  {};
23 
24  template<typename T>
25  struct HasSetReuse<T, decltype(std::declval<T>().setReuse(true), void())>
26  : std::true_type
27  {};
28 
29  template<typename T>
30  inline void setLinearSystemReuse(T& solver_backend, bool reuse, std::true_type)
31  {
32  if (!solver_backend.getReuse() && reuse)
33  dwarn << "WARNING: Newton needed to override your choice to reuse the linear system in order to work!" << std::endl;
34  solver_backend.setReuse(reuse);
35  }
36 
37  template<typename T>
38  inline void setLinearSystemReuse(T&, bool, std::false_type)
39  {}
40 
41  template<typename T>
42  inline void setLinearSystemReuse(T& solver_backend, bool reuse)
43  {
44  setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
45  }
46  }
47 
48 
61  template <typename GridOperator_, typename LinearSolver_>
63  {
64  public:
66  using GridOperator = GridOperator_;
67 
69  using LinearSolver = LinearSolver_;
70 
73 
76 
79 
81  using Real = typename Dune::FieldTraits<typename Domain::ElementType>::real_type;
82 
85 
87  const Result& result() const
88  {
89  if (not _resultValid)
90  DUNE_THROW(NewtonError, "NewtonMethod::result() called before NewtonMethod::solve()");
91  return _result;
92  }
93 
94  virtual void prepareStep(Domain& solution)
95  {
96  _reassembled = false;
97  if (_result.defect/_previousDefect > _reassembleThreshold){
98  if (_hangingNodeModifications){
99  auto dirichletValues = solution;
100  // Set all non dirichlet values to zero
101  Dune::PDELab::set_shifted_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, dirichletValues);
102  // Set all constrained DOFs to zero in solution
103  Dune::PDELab::set_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, solution);
104  // Copy correct Dirichlet values back into solution vector
105  Dune::PDELab::copy_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), dirichletValues, solution);
106  // Interpolate periodic constraints / hanging nodes
107  _gridOperator.localAssembler().backtransform(solution);
108  }
109  if (_verbosity>=3)
110  std::cout << " Reassembling matrix..." << std::endl;
111  *_jacobian = Real(0.0);
112  _gridOperator.jacobian(solution, *_jacobian);
113  _reassembled = true;
114  }
115 
116  _linearReduction = _minLinearReduction;
117  if (not _fixedLinearReduction){
118  // Determine maximum defect, where Newton is converged.
119  using std::min;
120  using std::max;
121  auto stop_defect = max(_result.first_defect*_reduction, _absoluteLimit);
122 
123  // To achieve second order convergence of newton we need a linear
124  // reduction of at least current_defect^2/prev_defect^2. For the
125  // last newton step a linear reduction of
126  // 1/10*end_defect/current_defect is sufficient for convergence.
127  if (stop_defect/(10*_result.defect) > _result.defect*_result.defect/(_previousDefect*_previousDefect))
128  _linearReduction = stop_defect/(10*_result.defect);
129  else
130  _linearReduction = min(_minLinearReduction, _result.defect*_result.defect/(_previousDefect*_previousDefect));
131  }
132 
133  if (_verbosity >= 3)
134  std::cout << " requested linear reduction: "
135  << std::setw(12) << std::setprecision(4) << std::scientific
136  << _linearReduction << std::endl;
137  }
138 
139  virtual void linearSolve()
140  {
141  if (_verbosity >= 4)
142  std::cout << " Solving linear system..." << std::endl;
143 
144  // If the jacobian was not reassembled we might save some work in the solver backend
145  Impl::setLinearSystemReuse(_linearSolver, not _reassembled);
146 
147  // Solve the linear system
148  _correction = 0.0;
149  _linearSolver.apply(*_jacobian, _correction, _residual, _linearReduction);
150 
151  if (not _linearSolver.result().converged)
152  DUNE_THROW(NewtonLinearSolverError,
153  "NewtonMethod::linearSolve(): Linear solver did not converge "
154  "in " << _linearSolver.result().iterations << " iterations");
155  if (_verbosity >= 4)
156  std::cout << " linear solver iterations: "
157  << std::setw(12) << _linearSolver.result().iterations << std::endl
158  << " linear defect reduction: "
159  << std::setw(12) << std::setprecision(4) << std::scientific
160  << _linearSolver.result().reduction << std::endl;
161  }
162 
164  virtual void apply(Domain& solution)
165  {
166  // Reset solver statistics
167  _result.clear();
168  _resultValid = true;
169 
170  // Store old ios flags (will be reset when this goes out of scope)
171  ios_base_all_saver restorer(std::cout);
172 
173  // Prepare time measuring
174  using Clock = std::chrono::steady_clock;
175  using Duration = Clock::duration;
176  auto assembler_time = Duration::zero();
177  auto linear_solver_time = Duration::zero();
178  auto line_search_time = Duration::zero();
179  auto to_seconds =
180  [](Duration duration){
181  return std::chrono::duration<double>(duration).count();
182  };
183  auto start_solve = Clock::now();
184 
185  //=========================
186  // Calculate initial defect
187  //=========================
188  updateDefect(solution);
189  _result.first_defect = _result.defect;
190  _previousDefect = _result.defect;
191 
192  if (_verbosity >= 2)
193  std::cout << " Initial defect: "
194  << std::setw(12) << std::setprecision(4) << std::scientific
195  << _result.defect << std::endl;
196 
197  //==========================
198  // Calculate Jacobian matrix
199  //==========================
200  if (not _jacobian)
201  _jacobian = std::make_shared<Jacobian>(_gridOperator);
202 
203  //=========================
204  // Nonlinear iteration loop
205  //=========================
206  while (not _terminate->terminate()){
207  if(_verbosity >= 3)
208  std::cout << " Newton iteration " << _result.iterations
209  << " --------------------------------" << std::endl;
210 
211  //=============
212  // Prepare step
213  //=============
214  auto start = Clock::now();
215  prepareStep(solution);
216  auto end = Clock::now();
217  assembler_time += end -start;
218 
219  // Store defect
220  _previousDefect = _result.defect;
221 
222  //====================
223  // Solve linear system
224  //====================
225  start = Clock::now();
226  linearSolve();
227  end = Clock::now();
228  linear_solver_time += end -start;
229 
230  //===================================
231  // Do line search and update solution
232  //===================================
233  start = Clock::now();
234  _lineSearch->lineSearch(solution, _correction);
235  // lineSearch(solution);
236  end = Clock::now();
237  line_search_time += end -start;
238 
239  //========================================
240  // Store statistics and create some output
241  //========================================
242  _result.reduction = _result.defect/_result.first_defect;
243  _result.iterations++;
244  _result.conv_rate = std::pow(_result.reduction, 1.0/_result.iterations);
245 
246  if (_verbosity >= 3)
247  std::cout << " linear solver time: "
248  << std::setw(12) << std::setprecision(4) << std::scientific
249  << to_seconds(linear_solver_time) << std::endl
250  << " defect reduction (this iteration):"
251  << std::setw(12) << std::setprecision(4) << std::scientific
252  << _result.defect/_previousDefect << std::endl
253  << " defect reduction (total): "
254  << std::setw(12) << std::setprecision(4) << std::scientific
255  << _result.reduction << std::endl
256  << " new defect: "
257  << std::setw(12) << std::setprecision(4) << std::scientific
258  << _result.defect << std::endl;
259  if (_verbosity == 2)
260  std::cout << " Newton iteration "
261  << std::setw(2)
262  << _result.iterations
263  << ". New defect: "
264  << std::setw(12) << std::setprecision(4) << std::scientific
265  << _result.defect
266  << ". Reduction (this): "
267  << std::setw(12) << std::setprecision(4) << std::scientific
268  << _result.defect/_previousDefect
269  << ". Reduction (total): "
270  << std::setw(12) << std::setprecision(4) << std::scientific
271  << _result.reduction
272  << std::endl;
273  } // end while loop of nonlinear iterations
274 
275  auto end_solve = Clock::now();
276  auto solve_time = end_solve - start_solve;
277  _result.elapsed = to_seconds(solve_time);
278 
279  if (_verbosity == 1)
280  std::cout << " Newton converged after "
281  << std::setw(2)
282  << _result.iterations
283  << " iterations. Reduction: "
284  << std::setw(12) << std::setprecision(4) << std::scientific
285  << _result.reduction
286  << " (" << std::setprecision(4) << _result.elapsed << "s)"
287  << std::endl;
288 
289  if (not _keepMatrix)
290  _jacobian.reset();
291  }
292 
294  virtual void updateDefect(Domain& solution)
295  {
296  if (_hangingNodeModifications){
297  auto dirichletValues = solution;
298  // Set all non dirichlet values to zero
299  Dune::PDELab::set_shifted_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, dirichletValues);
300  // Set all constrained DOFs to zero in solution
301  Dune::PDELab::set_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, solution);
302  // Copy correct Dirichlet values back into solution vector
303  Dune::PDELab::copy_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), dirichletValues, solution);
304  // Interpolate periodic constraints / hanging nodes
305  _gridOperator.localAssembler().backtransform(solution);
306  }
307 
308  _residual = 0.0;
309  _gridOperator.residual(solution, _residual);
310 
311  // Use the maximum norm as a stopping criterion. This helps loosen the tolerance
312  // when solving for stationary solutions of nonlinear time-dependent problems.
313  // The default is to use the Euclidean norm (in the else-block) as before
314  if (_useMaxNorm){
315  auto rankMax = Backend::native(_residual).infinity_norm();
316  _result.defect = _gridOperator.testGridFunctionSpace().gridView().comm().max(rankMax);
317  }
318  else
319  _result.defect = _linearSolver.norm(_residual);
320  }
321 
323  void setVerbosityLevel(unsigned int verbosity)
324  {
325  if (_gridOperator.trialGridFunctionSpace().gridView().comm().rank()>0)
326  _verbosity = 0;
327  else
328  _verbosity = verbosity;
329  }
330 
332  unsigned int getVerbosityLevel() const
333  {
334  return _verbosity;
335  }
336 
338  void setReduction(Real reduction)
339  {
340  _reduction = reduction;
341  }
342 
345  {
346  return _reduction;
347  }
348 
350  void setAbsoluteLimit(Real absoluteLimit)
351  {
352  _absoluteLimit = absoluteLimit;
353  }
354 
356  {
357  return _absoluteLimit;
358  }
359 
361  void setKeepMatrix(bool b)
362  {
363  _keepMatrix = b;
364  }
365 
367  void setUseMaxNorm(bool b)
368  {
369  _useMaxNorm = b;
370  }
371 
374  {
375  _hangingNodeModifications = b;
376  }
377 
378 
380  bool keepMatrix() const
381  {
382  return _keepMatrix;
383  }
384 
387  {
388  if(_jacobian)
389  _jacobian.reset();
390  }
391 
399  void setMinLinearReduction(Real minLinearReduction)
400  {
401  _minLinearReduction = minLinearReduction;
402  }
403 
409  void setFixedLinearReduction(bool fixedLinearReduction)
410  {
411  _fixedLinearReduction = fixedLinearReduction;
412  }
413 
420  void setReassembleThreshold(Real reassembleThreshold)
421  {
422  _reassembleThreshold = reassembleThreshold;
423  }
424 
458  void setParameters(const ParameterTree& parameterTree){
459  if (parameterTree.hasKey("verbosity"))
460  setVerbosityLevel(parameterTree.get<unsigned int>("verbosity"));
461  if (parameterTree.hasKey("reduction"))
462  setReduction(parameterTree.get<Real>("reduction"));
463  if (parameterTree.hasKey("absolute_limit"))
464  setAbsoluteLimit(parameterTree.get<Real>("absolute_limit"));
465  if (parameterTree.hasKey("keeep_matrix"))
466  setKeepMatrix(parameterTree.get<bool>("keep_matrix"));
467  if (parameterTree.hasKey("use_max_norm"))
468  setUseMaxNorm(parameterTree.get<bool>("use_max_norm"));
469  if (parameterTree.hasKey("hanging_node_modifications"))
470  _hangingNodeModifications = parameterTree.get<bool>("hanging_node_modifications");
471 
472  if (parameterTree.hasKey("min_linear_reduction"))
473  setMinLinearReduction(parameterTree.get<Real>("min_linear_reduction"));
474  if (parameterTree.hasKey("fixed_linear_reduction"))
475  setFixedLinearReduction(parameterTree.get<bool>("fixed_linear_reduction"));
476  if (parameterTree.hasKey("reassemble_threshold"))
477  setReassembleThreshold(parameterTree.get<Real>("reassemble_threshold"));
478 
479  _terminate->setParameters(parameterTree.sub("terminate"));
480  _lineSearch->setParameters(parameterTree.sub("line_search"));
481  }
482 
484  void setTerminate(std::shared_ptr<TerminateInterface> terminate)
485  {
486  _terminate = terminate;
487  }
488 
493  void setLineSearch(std::shared_ptr<LineSearchInterface<Domain>> lineSearch)
494  {
495  _lineSearch = lineSearch;
496  }
497 
500  const GridOperator& gridOperator,
501  LinearSolver& linearSolver,
502  const std::string& lineSearchStrategy="hackbusch_reusken")
503  : _gridOperator(gridOperator)
504  , _linearSolver(linearSolver)
505  , _residual(gridOperator.testGridFunctionSpace())
506  , _correction(gridOperator.trialGridFunctionSpace())
507  {
508  _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
509  _lineSearch = getLineSearch(*this, lineSearchStrategy);
510  }
511 
514  const GridOperator& gridOperator,
515  LinearSolver& linearSolver,
516  const ParameterTree& parameterTree,
517  const std::string& lineSearchStrategy="hackbusch_reusken")
518  : _gridOperator(gridOperator)
519  , _linearSolver(linearSolver)
520  , _residual(gridOperator.testGridFunctionSpace())
521  , _correction(gridOperator.trialGridFunctionSpace())
522 
523  {
524  setParameters(parameterTree);
525  _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
526  _lineSearch = getLineSearch(*this, lineSearchStrategy);
527  }
528 
529  private:
530  const GridOperator& _gridOperator;
531  LinearSolver& _linearSolver;
532 
533  // Vectors and Jacobi matrix we set up only once
534  Range _residual;
535  Domain _correction;
536  std::shared_ptr<Jacobian> _jacobian;
537  std::shared_ptr<Domain> _previousSolution;
538 
539  std::shared_ptr<TerminateInterface> _terminate;
540  std::shared_ptr<LineSearchInterface<Domain>> _lineSearch;
541 
542  // Class for storing results
543  Result _result;
544  bool _resultValid = false; // result class only valid after calling apply
545  Real _previousDefect = 0.0;
546 
547  // Remember if jacobian was reassembled in prepareStep
548  bool _reassembled = true; // will be set in prepare step
549  Real _linearReduction = 0.0; // will be set in prepare step
550 
551  // User parameters
552  unsigned int _verbosity = 0;
553  Real _reduction = 1e-8;
554  Real _absoluteLimit = 1e-12;
555  bool _keepMatrix = true;
556  bool _useMaxNorm = false;
557 
558  // Special treatment if we have hanging nodes
559  bool _hangingNodeModifications = false;
560 
561  // User parameters for prepareStep()
562  Real _minLinearReduction = 1e-3;
563  bool _fixedLinearReduction = false;
564  Real _reassembleThreshold = 0.0;
565  };
566 
567 } // namespace Dune::PDELab
568 
569 #endif
const Entity & e
Definition: localfunctionspace.hh:121
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1014
void copy_constrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:936
Definition: adaptivity.hh:29
std::shared_ptr< LineSearchInterface< typename Newton::Domain > > getLineSearch(Newton &newton, const std::string &name)
Get a pointer to a line search.
Definition: newtonlinesearch.hh:225
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
RFType conv_rate
Definition: solver.hh:36
RFType reduction
Definition: solver.hh:35
unsigned int iterations
Definition: solver.hh:33
double elapsed
Definition: solver.hh:34
Dune::PDELab::Backend::Matrix< MB, Domain, Range, JF > Jacobian
The type of the jacobian.
Definition: gridoperatorutilities.hh:72
Dune::PDELab::Backend::Vector< GFSV, RF > Range
The type of the range (residual).
Definition: gridoperatorutilities.hh:65
Dune::PDELab::Backend::Vector< GFSU, DF > Domain
The type of the domain (solution).
Definition: gridoperatorutilities.hh:58
Standard grid operator implementation.
Definition: gridoperator.hh:36
Newton solver for solving non-linear problems.
Definition: solver/newton.hh:63
const Result & result() const
Return results.
Definition: solver/newton.hh:87
typename GridOperator::Traits::Jacobian Jacobian
Type of the Jacobian matrix.
Definition: solver/newton.hh:78
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver, const ParameterTree &parameterTree, const std::string &lineSearchStrategy="hackbusch_reusken")
Construct Newton passing a parameter tree.
Definition: solver/newton.hh:513
void setMinLinearReduction(Real minLinearReduction)
Set the minimal reduction in the linear solver.
Definition: solver/newton.hh:399
void setTerminate(std::shared_ptr< TerminateInterface > terminate)
Set the termination criterion.
Definition: solver/newton.hh:484
virtual void linearSolve()
Definition: solver/newton.hh:139
Real getAbsoluteLimit() const
Definition: solver/newton.hh:355
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: solver/newton.hh:380
unsigned int getVerbosityLevel() const
Get verbosity level.
Definition: solver/newton.hh:332
typename GridOperator::Traits::Domain Domain
Type of the domain (solution)
Definition: solver/newton.hh:72
void setParameters(const ParameterTree &parameterTree)
Interpret a parameter tree as a set of options for the newton solver.
Definition: solver/newton.hh:458
virtual void prepareStep(Domain &solution)
Definition: solver/newton.hh:94
void setHangingNodeModifications(bool b)
Does the problem have hanging nodes.
Definition: solver/newton.hh:373
void setFixedLinearReduction(bool fixedLinearReduction)
Set wether to use a fixed reduction in the linear solver.
Definition: solver/newton.hh:409
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: solver/newton.hh:386
GridOperator_ GridOperator
Type of the grid operator.
Definition: solver/newton.hh:66
void setReduction(Real reduction)
Set reduction Newton needs to achieve.
Definition: solver/newton.hh:338
void setLineSearch(std::shared_ptr< LineSearchInterface< Domain >> lineSearch)
Set the line search.
Definition: solver/newton.hh:493
typename Dune::FieldTraits< typename Domain::ElementType >::real_type Real
Number type.
Definition: solver/newton.hh:81
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: solver/newton.hh:361
LinearSolver_ LinearSolver
Type of the linear solver.
Definition: solver/newton.hh:69
void setUseMaxNorm(bool b)
Set whether to use the maximum norm for stopping criteria.
Definition: solver/newton.hh:367
typename GridOperator::Traits::Range Range
Type of the range (residual)
Definition: solver/newton.hh:75
void setAbsoluteLimit(Real absoluteLimit)
Set absolute convergence limit.
Definition: solver/newton.hh:350
virtual void updateDefect(Domain &solution)
Update _residual and defect in _result.
Definition: solver/newton.hh:294
virtual void apply(Domain &solution)
Solve the nonlinear problem using solution as initial guess and for storing the result.
Definition: solver/newton.hh:164
PDESolverResult< Real > Result
Type of results.
Definition: solver/newton.hh:84
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver, const std::string &lineSearchStrategy="hackbusch_reusken")
Construct Newton using default parameters.
Definition: solver/newton.hh:499
void setReassembleThreshold(Real reassembleThreshold)
Set a threshold, when the linear operator is reassembled.
Definition: solver/newton.hh:420
void setVerbosityLevel(unsigned int verbosity)
Set how much output you get.
Definition: solver/newton.hh:323
Real getReduction() const
Get reduction.
Definition: solver/newton.hh:344
Abstract base class describing the line search interface.
Definition: newtonlinesearch.hh:13
RFType defect
Definition: solver/utility.hh:11
RFType first_defect
Definition: solver/utility.hh:10
void clear()
Definition: solver/utility.hh:21