dune-pdelab  2.5-dev
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_NEWTON_NEWTON_HH
4 #define DUNE_PDELAB_NEWTON_NEWTON_HH
5 
6 #include <iostream>
7 #include <iomanip>
8 #include <cmath>
9 #include <memory>
10 
11 #include <type_traits>
12 
13 #include <math.h>
14 
15 #include <dune/common/debugstream.hh>
16 #include <dune/common/exceptions.hh>
17 #include <dune/common/ios_state.hh>
18 #include <dune/common/timer.hh>
19 #include <dune/common/parametertree.hh>
20 
22 
23 namespace Dune
24 {
25  namespace PDELab
26  {
27  namespace Impl
28  {
29  // Some SFinae magic to execute setReuse(bool) on a backend
30  template<typename T1, typename = void>
31  struct HasSetReuse
32  : std::false_type
33  {};
34 
35  template<typename T>
36  struct HasSetReuse<T, decltype(std::declval<T>().setReuse(true), void())>
37  : std::true_type
38  {};
39 
40  template<typename T>
41  inline void setLinearSystemReuse(T& solver_backend, bool reuse, std::true_type)
42  {
43  if (!solver_backend.getReuse() && reuse)
44  dwarn << "WARNING: Newton needed to override your choice to reuse the linear system in order to work!" << std::endl;
45  solver_backend.setReuse(reuse);
46  }
47 
48  template<typename T>
49  inline void setLinearSystemReuse(T&, bool, std::false_type)
50  {}
51 
52  template<typename T>
53  inline void setLinearSystemReuse(T& solver_backend, bool reuse)
54  {
55  setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
56  }
57  }
58 
59  // Exception classes used in NewtonSolver
60  class NewtonError : public Exception {};
61  class NewtonDefectError : public NewtonError {};
64  class NewtonNotConverged : public NewtonError {};
65 
66  // Status information of Newton's method
67  template<class RFType>
69  {
70  RFType first_defect; // the first defect
71  RFType defect; // the final defect
72  double assembler_time; // Cumulative time for matrix assembly
73  double linear_solver_time; // Cumulative time for linear solver
74  int linear_solver_iterations; // Total number of linear iterations
75 
77  first_defect(0.0), defect(0.0), assembler_time(0.0), linear_solver_time(0.0),
78  linear_solver_iterations(0) {}
79  };
80 
81  template<class GOS, class TrlV, class TstV>
82  class NewtonBase
83  {
84  typedef GOS GridOperator;
85  typedef TrlV TrialVector;
86  typedef TstV TestVector;
87 
88  typedef typename TestVector::ElementType RFType;
89  typedef typename GOS::Traits::Jacobian Matrix;
90 
91 
92  public:
93  // export result type
95 
96  void setVerbosityLevel(unsigned int verbosity_level)
97  {
98  if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
99  verbosity_level_ = 0;
100  else
101  verbosity_level_ = verbosity_level;
102  }
103 
105  void setKeepMatrix(bool b)
106  {
107  keep_matrix_ = b;
108  }
109 
111  bool keepMatrix() const
112  {
113  return keep_matrix_;
114  }
115 
118  {
119  if(A_)
120  A_.reset();
121  }
122 
123  protected:
124  const GridOperator& gridoperator_;
125  TrialVector *u_;
126  std::shared_ptr<TrialVector> z_;
127  std::shared_ptr<TestVector> r_;
128  std::shared_ptr<Matrix> A_;
129  Result res_;
130  unsigned int verbosity_level_;
131  RFType prev_defect_;
134  RFType reduction_;
135  RFType abs_limit_;
137 
138  NewtonBase(const GridOperator& go, TrialVector& u)
139  : gridoperator_(go)
140  , u_(&u)
141  , verbosity_level_(1)
142  , keep_matrix_(true)
143  {
144  if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
145  verbosity_level_ = 0;
146  }
147 
148  NewtonBase(const GridOperator& go)
149  : gridoperator_(go)
150  , u_(0)
151  , verbosity_level_(1)
152  , keep_matrix_(true)
153  {
154  if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
155  verbosity_level_ = 0;
156  }
157 
158  virtual ~NewtonBase() { }
159 
160  virtual bool terminate() = 0;
161  virtual void prepare_step(Matrix& A, TestVector& r) = 0;
162  virtual void line_search(TrialVector& z, TestVector& r) = 0;
163  virtual void defect(TestVector& r) = 0;
164  }; // end class NewtonBase
165 
166  template<class GOS, class S, class TrlV, class TstV>
167  class NewtonSolver : public virtual NewtonBase<GOS,TrlV,TstV>
168  {
169  typedef S Solver;
170  typedef GOS GridOperator;
171  typedef TrlV TrialVector;
172  typedef TstV TestVector;
173 
174  typedef typename TestVector::ElementType RFType;
175  typedef typename GOS::Traits::Jacobian Matrix;
176 
177  public:
179 
180  NewtonSolver(const GridOperator& go, TrialVector& u_, Solver& solver)
181  : NewtonBase<GOS,TrlV,TstV>(go,u_)
182  , solver_(solver)
183  , result_valid_(false)
184  {}
185 
186  NewtonSolver(const GridOperator& go, Solver& solver)
187  : NewtonBase<GOS,TrlV,TstV>(go)
188  , solver_(solver)
189  , result_valid_(false)
190  {}
191 
192  void apply();
193 
194  void apply(TrialVector& u_);
195 
196  const Result& result() const
197  {
198  if (!result_valid_)
199  DUNE_THROW(NewtonError,
200  "NewtonSolver::result() called before NewtonSolver::solve()");
201  return this->res_;
202  }
203 
204  protected:
205  virtual void defect(TestVector& r)
206  {
207  r = 0.0;
208  this->gridoperator_.residual(*this->u_, r);
209  this->res_.defect = this->solver_.norm(r);
210  if (!std::isfinite(this->res_.defect))
211  DUNE_THROW(NewtonDefectError,
212  "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
213  }
214 
215 
216  private:
217  void linearSolve(Matrix& A, TrialVector& z, TestVector& r) const
218  {
219  if (this->verbosity_level_ >= 4)
220  std::cout << " Solving linear system..." << std::endl;
221  z = 0.0;
222  // If possible tell solver backend to reuse linear system when we did not reassemble.
223  Impl::setLinearSystemReuse(this->solver_, this->reassembled_);
224  this->solver_.apply(A, z, r, this->linear_reduction_);
225 
226  ios_base_all_saver restorer(std::cout); // store old ios flags
227 
228  if (!this->solver_.result().converged)
229  DUNE_THROW(NewtonLinearSolverError,
230  "NewtonSolver::linearSolve(): Linear solver did not converge "
231  "in " << this->solver_.result().iterations << " iterations");
232  if (this->verbosity_level_ >= 4)
233  std::cout << " linear solver iterations: "
234  << std::setw(12) << solver_.result().iterations << std::endl
235  << " linear defect reduction: "
236  << std::setw(12) << std::setprecision(4) << std::scientific
237  << solver_.result().reduction << std::endl;
238  }
239 
240  Solver& solver_;
241  bool result_valid_;
242  }; // end class NewtonSolver
243 
244  template<class GOS, class S, class TrlV, class TstV>
246  {
247  this->u_ = &u;
248  apply();
249  }
250 
251  template<class GOS, class S, class TrlV, class TstV>
253  {
254  this->res_.iterations = 0;
255  this->res_.converged = false;
256  this->res_.reduction = 1.0;
257  this->res_.conv_rate = 1.0;
258  this->res_.elapsed = 0.0;
259  this->res_.assembler_time = 0.0;
260  this->res_.linear_solver_time = 0.0;
261  this->res_.linear_solver_iterations = 0;
262  result_valid_ = true;
263  Timer timer;
264 
265  try
266  {
267  if(!this->r_) {
268  // std::cout << "=== Setting up residual vector ..." << std::endl;
269  this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
270  }
271  // residual calculation in member function "defect":
272  //--------------------------------------------------
273  // - set residual vector to zero
274  // - calculate new residual
275  // - store norm of residual in "this->res_.defect"
276  this->defect(*this->r_);
277  this->res_.first_defect = this->res_.defect;
278  this->prev_defect_ = this->res_.defect;
279 
280  if (this->verbosity_level_ >= 2)
281  {
282  // store old ios flags
283  ios_base_all_saver restorer(std::cout);
284  std::cout << " Initial defect: "
285  << std::setw(12) << std::setprecision(4) << std::scientific
286  << this->res_.defect << std::endl;
287  }
288 
289  if(!this->A_) {
290  // std::cout << "==== Setting up jacobian matrix ... " << std::endl;
291  this->A_ = std::make_shared<Matrix>(this->gridoperator_);
292  }
293  if(!this->z_) {
294  // std::cout << "==== Setting up correction vector ... " << std::endl;
295  this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
296  }
297 
298  while (!this->terminate())
299  {
300  if (this->verbosity_level_ >= 3)
301  std::cout << " Newton iteration " << this->res_.iterations
302  << " --------------------------------" << std::endl;
303 
304  Timer assembler_timer;
305  try
306  {
307  // jacobian calculation in member function "prepare_step"
308  //-------------------------------------------------------
309  // - if above reassemble threshold
310  // - set jacobian to zero
311  // - calculate new jacobian
312  // - set linear reduction
313  this->prepare_step(*this->A_,*this->r_);
314  }
315  catch (...)
316  {
317  this->res_.assembler_time += assembler_timer.elapsed();
318  throw;
319  }
320  double assembler_time = assembler_timer.elapsed();
321  this->res_.assembler_time += assembler_time;
322  if (this->verbosity_level_ >= 3)
323  std::cout << " matrix assembly time: "
324  << std::setw(12) << std::setprecision(4) << std::scientific
325  << assembler_time << std::endl;
326 
327  Timer linear_solver_timer;
328  try
329  {
330  // solution of linear system in member function "linearSolve"
331  //-----------------------------------------------------------
332  // - set initial guess for correction z to zero
333  // - call linear solver
334  this->linearSolve(*this->A_, *this->z_, *this->r_);
335  }
336  catch (...)
337  {
338  this->res_.linear_solver_time += linear_solver_timer.elapsed();
339  this->res_.linear_solver_iterations += this->solver_.result().iterations;
340  throw;
341  }
342  double linear_solver_time = linear_solver_timer.elapsed();
343  this->res_.linear_solver_time += linear_solver_time;
344  this->res_.linear_solver_iterations += this->solver_.result().iterations;
345 
346  try
347  {
348  // line search with correction z
349  // the undamped version is also integrated in here
350  this->line_search(*this->z_, *this->r_);
351  }
352  catch (NewtonLineSearchError)
353  {
354  if (this->reassembled_)
355  throw;
356  if (this->verbosity_level_ >= 3)
357  std::cout << " line search failed - trying again with reassembled matrix" << std::endl;
358  continue;
359  }
360 
361  this->res_.reduction = this->res_.defect/this->res_.first_defect;
362  this->res_.iterations++;
363  this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
364 
365  // store old ios flags
366  ios_base_all_saver restorer(std::cout);
367 
368  if (this->verbosity_level_ >= 3)
369  std::cout << " linear solver time: "
370  << std::setw(12) << std::setprecision(4) << std::scientific
371  << linear_solver_time << std::endl
372  << " defect reduction (this iteration):"
373  << std::setw(12) << std::setprecision(4) << std::scientific
374  << this->res_.defect/this->prev_defect_ << std::endl
375  << " defect reduction (total): "
376  << std::setw(12) << std::setprecision(4) << std::scientific
377  << this->res_.reduction << std::endl
378  << " new defect: "
379  << std::setw(12) << std::setprecision(4) << std::scientific
380  << this->res_.defect << std::endl;
381  if (this->verbosity_level_ == 2)
382  std::cout << " Newton iteration " << std::setw(2) << this->res_.iterations
383  << ". New defect: "
384  << std::setw(12) << std::setprecision(4) << std::scientific
385  << this->res_.defect
386  << ". Reduction (this): "
387  << std::setw(12) << std::setprecision(4) << std::scientific
388  << this->res_.defect/this->prev_defect_
389  << ". Reduction (total): "
390  << std::setw(12) << std::setprecision(4) << std::scientific
391  << this->res_.reduction << std::endl;
392  } // end while
393  } // end try
394  catch(...)
395  {
396  this->res_.elapsed = timer.elapsed();
397  throw;
398  }
399  this->res_.elapsed = timer.elapsed();
400 
401  ios_base_all_saver restorer(std::cout); // store old ios flags
402 
403  if (this->verbosity_level_ == 1)
404  std::cout << " Newton converged after " << std::setw(2) << this->res_.iterations
405  << " iterations. Reduction: "
406  << std::setw(12) << std::setprecision(4) << std::scientific
407  << this->res_.reduction
408  << " (" << std::setprecision(4) << this->res_.elapsed << "s)"
409  << std::endl;
410 
411  if(!this->keep_matrix_)
412  this->A_.reset();
413  } // end apply
414 
415  template<class GOS, class TrlV, class TstV>
416  class NewtonTerminate : public virtual NewtonBase<GOS,TrlV,TstV>
417  {
418  typedef GOS GridOperator;
419  typedef TrlV TrialVector;
420 
421  typedef typename TstV::ElementType RFType;
422 
423  public:
424  NewtonTerminate(const GridOperator& go, TrialVector& u_)
425  : NewtonBase<GOS,TrlV,TstV>(go,u_)
426  , maxit_(40)
427  , force_iteration_(false)
428  {
429  this->reduction_ = 1e-8;
430  this->abs_limit_ = 1e-12;
431  }
432 
433  NewtonTerminate(const GridOperator& go)
434  : NewtonBase<GOS,TrlV,TstV>(go)
435  , maxit_(40)
436  , force_iteration_(false)
437  {
438  this->reduction_ = 1e-8;
439  this->abs_limit_ = 1e-12;
440  }
441 
442  void setReduction(RFType reduction)
443  {
444  this->reduction_ = reduction;
445  }
446 
447  void setMaxIterations(unsigned int maxit)
448  {
449  maxit_ = maxit;
450  }
451 
452  void setForceIteration(bool force_iteration)
453  {
454  force_iteration_ = force_iteration;
455  }
456 
457  void setAbsoluteLimit(RFType abs_limit_)
458  {
459  this->abs_limit_ = abs_limit_;
460  }
461 
462  virtual bool terminate()
463  {
464  if (force_iteration_ && this->res_.iterations == 0)
465  return false;
466  this->res_.converged = this->res_.defect < this->abs_limit_
467  || this->res_.defect < this->res_.first_defect * this->reduction_;
468  if (this->res_.iterations >= maxit_ && !this->res_.converged)
469  DUNE_THROW(NewtonNotConverged,
470  "NewtonTerminate::terminate(): Maximum iteration count reached");
471  return this->res_.converged;
472  }
473 
474  private:
475  unsigned int maxit_;
476  bool force_iteration_;
477  }; // end class NewtonTerminate
478 
479  template<class GOS, class TrlV, class TstV>
480  class NewtonPrepareStep : public virtual NewtonBase<GOS,TrlV,TstV>
481  {
482  typedef GOS GridOperator;
483  typedef TrlV TrialVector;
484 
485  typedef typename TstV::ElementType RFType;
486  typedef typename GOS::Traits::Jacobian Matrix;
487 
488  public:
489  NewtonPrepareStep(const GridOperator& go, TrialVector& u_)
490  : NewtonBase<GOS,TrlV,TstV>(go,u_)
491  , min_linear_reduction_(1e-3)
492  , fixed_linear_reduction_(0.0)
493  , reassemble_threshold_(0.0)
494  {}
495 
496  NewtonPrepareStep(const GridOperator& go)
497  : NewtonBase<GOS,TrlV,TstV>(go)
498  , min_linear_reduction_(1e-3)
499  , fixed_linear_reduction_(0.0)
500  , reassemble_threshold_(0.0)
501  {}
502 
509  void setMinLinearReduction(RFType min_linear_reduction)
510  {
511  min_linear_reduction_ = min_linear_reduction;
512  }
513 
518  void setFixedLinearReduction(bool fixed_linear_reduction)
519  {
520  fixed_linear_reduction_ = fixed_linear_reduction;
521  }
522 
530  void setReassembleThreshold(RFType reassemble_threshold)
531  {
532  reassemble_threshold_ = reassemble_threshold;
533  }
534 
535  virtual void prepare_step(Matrix& A, TstV& )
536  {
537  this->reassembled_ = false;
538  if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
539  {
540  if (this->verbosity_level_ >= 3)
541  std::cout << " Reassembling matrix..." << std::endl;
542  A = 0.0;
543  this->gridoperator_.jacobian(*this->u_, A);
544  this->reassembled_ = true;
545  }
546 
547  if (fixed_linear_reduction_ == true)
548  this->linear_reduction_ = min_linear_reduction_;
549  else {
550  // determine maximum defect, where Newton is converged.
551  RFType stop_defect =
552  std::max(this->res_.first_defect * this->reduction_,
553  this->abs_limit_);
554 
555  /*
556  To achieve second order convergence of newton
557  we need a linear reduction of at least
558  current_defect^2/prev_defect^2.
559  For the last newton step a linear reduction of
560  1/10*end_defect/current_defect
561  is sufficient for convergence.
562  */
563  if ( stop_defect/(10*this->res_.defect) >
564  this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_) )
565  this->linear_reduction_ =
566  stop_defect/(10*this->res_.defect);
567  else
568  this->linear_reduction_ =
569  std::min(min_linear_reduction_,this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_));
570  }
571 
572  this->prev_defect_ = this->res_.defect;
573 
574  ios_base_all_saver restorer(std::cout); // store old ios flags
575 
576  if (this->verbosity_level_ >= 3)
577  std::cout << " requested linear reduction: "
578  << std::setw(12) << std::setprecision(4) << std::scientific
579  << this->linear_reduction_ << std::endl;
580  }
581 
582  private:
583  RFType min_linear_reduction_;
584  bool fixed_linear_reduction_;
585  RFType reassemble_threshold_;
586  }; // end class NewtonPrepareStep
587 
588  template<class GOS, class TrlV, class TstV>
589  class NewtonLineSearch : public virtual NewtonBase<GOS,TrlV,TstV>
590  {
591  typedef GOS GridOperator;
592  typedef TrlV TrialVector;
593  typedef TstV TestVector;
594 
595  typedef typename TestVector::ElementType RFType;
596 
597  public:
598  enum Strategy {
606  hackbuschReuskenAcceptBest };
607 
608  NewtonLineSearch(const GridOperator& go, TrialVector& u_)
609  : NewtonBase<GOS,TrlV,TstV>(go,u_)
610  , strategy_(hackbuschReusken)
611  , maxit_(10)
612  , damping_factor_(0.5)
613  {}
614 
615  NewtonLineSearch(const GridOperator& go)
616  : NewtonBase<GOS,TrlV,TstV>(go)
617  , strategy_(hackbuschReusken)
618  , maxit_(10)
619  , damping_factor_(0.5)
620  {}
621 
623  {
624  strategy_ = strategy;
625  }
626 
627  void setLineSearchStrategy(std::string strategy)
628  {
629  strategy_ = strategyFromName(strategy);
630  }
631 
632  void setLineSearchMaxIterations(unsigned int maxit)
633  {
634  maxit_ = maxit;
635  }
636 
637  void setLineSearchDampingFactor(RFType damping_factor)
638  {
639  damping_factor_ = damping_factor;
640  }
641 
642  virtual void line_search(TrialVector& z, TestVector& r)
643  {
644  if (strategy_ == noLineSearch)
645  {
646  this->u_->axpy(-1.0, z);
647  this->defect(r);
648  return;
649  }
650 
651  if (this->verbosity_level_ >= 4)
652  std::cout << " Performing line search..." << std::endl;
653  RFType lambda = 1.0;
654  RFType best_lambda = 0.0;
655  RFType best_defect = this->res_.defect;
656  TrialVector prev_u(*this->u_);
657  unsigned int i = 0;
658  ios_base_all_saver restorer(std::cout); // store old ios flags
659 
660  while (1)
661  {
662  if (this->verbosity_level_ >= 4)
663  std::cout << " trying line search damping factor: "
664  << std::setw(12) << std::setprecision(4) << std::scientific
665  << lambda
666  << std::endl;
667 
668  this->u_->axpy(-lambda, z);
669  try {
670  this->defect(r);
671  }
672  catch (NewtonDefectError)
673  {
674  if (this->verbosity_level_ >= 4)
675  std::cout << " NaNs detected" << std::endl;
676  } // ignore NaNs and try again with lower lambda
677 
678  if (this->res_.defect <= (1.0 - lambda/4) * this->prev_defect_)
679  {
680  if (this->verbosity_level_ >= 4)
681  std::cout << " line search converged" << std::endl;
682  break;
683  }
684 
685  if (this->res_.defect < best_defect)
686  {
687  best_defect = this->res_.defect;
688  best_lambda = lambda;
689  }
690 
691  if (++i >= maxit_)
692  {
693  if (this->verbosity_level_ >= 4)
694  std::cout << " max line search iterations exceeded" << std::endl;
695  switch (strategy_)
696  {
697  case hackbuschReusken:
698  *this->u_ = prev_u;
699  this->defect(r);
700  DUNE_THROW(NewtonLineSearchError,
701  "NewtonLineSearch::line_search(): line search failed, "
702  "max iteration count reached, "
703  "defect did not improve enough");
704  case hackbuschReuskenAcceptBest:
705  if (best_lambda == 0.0)
706  {
707  *this->u_ = prev_u;
708  this->defect(r);
709  DUNE_THROW(NewtonLineSearchError,
710  "NewtonLineSearch::line_search(): line search failed, "
711  "max iteration count reached, "
712  "defect did not improve in any of the iterations");
713  }
714  if (best_lambda != lambda)
715  {
716  *this->u_ = prev_u;
717  this->u_->axpy(-best_lambda, z);
718  this->defect(r);
719  }
720  break;
721  case noLineSearch:
722  break;
723  }
724  break;
725  }
726 
727  lambda *= damping_factor_;
728  *this->u_ = prev_u;
729  }
730  if (this->verbosity_level_ >= 4)
731  std::cout << " line search damping factor: "
732  << std::setw(12) << std::setprecision(4) << std::scientific
733  << lambda << std::endl;
734  } // end line_search
735 
736  protected:
738  Strategy strategyFromName(const std::string & s) {
739  if (s == "noLineSearch")
740  return noLineSearch;
741  if (s == "hackbuschReusken")
742  return hackbuschReusken;
743  if (s == "hackbuschReuskenAcceptBest")
744  return hackbuschReuskenAcceptBest;
745  DUNE_THROW(Exception, "unknown line search strategy" << s);
746  }
747 
748  private:
749  Strategy strategy_;
750  unsigned int maxit_;
751  RFType damping_factor_;
752  }; // end class NewtonLineSearch
753 
754  template<class GOS, class S, class TrlV, class TstV = TrlV>
755  class Newton : public NewtonSolver<GOS,S,TrlV,TstV>
756  , public NewtonTerminate<GOS,TrlV,TstV>
757  , public NewtonLineSearch<GOS,TrlV,TstV>
758  , public NewtonPrepareStep<GOS,TrlV,TstV>
759  {
760  typedef GOS GridOperator;
761  typedef S Solver;
762  typedef TrlV TrialVector;
763 
764  public:
765  Newton(const GridOperator& go, TrialVector& u_, Solver& solver_)
766  : NewtonBase<GOS,TrlV,TstV>(go,u_)
767  , NewtonSolver<GOS,S,TrlV,TstV>(go,u_,solver_)
768  , NewtonTerminate<GOS,TrlV,TstV>(go,u_)
769  , NewtonLineSearch<GOS,TrlV,TstV>(go,u_)
770  , NewtonPrepareStep<GOS,TrlV,TstV>(go,u_)
771  {}
772  Newton(const GridOperator& go, Solver& solver_)
773  : NewtonBase<GOS,TrlV,TstV>(go)
774  , NewtonSolver<GOS,S,TrlV,TstV>(go,solver_)
775  , NewtonTerminate<GOS,TrlV,TstV>(go)
776  , NewtonLineSearch<GOS,TrlV,TstV>(go)
777  , NewtonPrepareStep<GOS,TrlV,TstV>(go)
778  {}
780 
801  void setParameters(const Dune::ParameterTree & param)
802  {
803  typedef typename TstV::ElementType RFType;
804  if (param.hasKey("VerbosityLevel"))
805  this->setVerbosityLevel(
806  param.get<unsigned int>("VerbosityLevel"));
807  if (param.hasKey("Reduction"))
808  this->setReduction(
809  param.get<RFType>("Reduction"));
810  if (param.hasKey("MaxIterations"))
811  this->setMaxIterations(
812  param.get<unsigned int>("MaxIterations"));
813  if (param.hasKey("ForceIteration"))
814  this->setForceIteration(
815  param.get<bool>("ForceIteration"));
816  if (param.hasKey("AbsoluteLimit"))
817  this->setAbsoluteLimit(
818  param.get<RFType>("AbsoluteLimit"));
819  if (param.hasKey("MinLinearReduction"))
820  this->setMinLinearReduction(
821  param.get<RFType>("MinLinearReduction"));
822  if (param.hasKey("FixedLinearReduction"))
823  this->setFixedLinearReduction(
824  param.get<bool>("FixedLinearReduction"));
825  if (param.hasKey("ReassembleThreshold"))
826  this->setReassembleThreshold(
827  param.get<RFType>("ReassembleThreshold"));
828  if (param.hasKey("LineSearchStrategy"))
829  this->setLineSearchStrategy(
830  param.get<std::string>("LineSearchStrategy"));
831  if (param.hasKey("LineSearchMaxIterations"))
832  this->setLineSearchMaxIterations(
833  param.get<unsigned int>("LineSearchMaxIterations"));
834  if (param.hasKey("LineSearchDampingFactor"))
835  this->setLineSearchDampingFactor(
836  param.get<RFType>("LineSearchDampingFactor"));
837  if (param.hasKey("KeepMatrix"))
838  this->setKeepMatrix(
839  param.get<bool>("KeepMatrix"));
840  }
841  }; // end class Newton
842  } // end namespace PDELab
843 } // end namespace Dune
844 
845 #endif // DUNE_PDELAB_NEWTON_NEWTON_HH
NewtonResult< RFType > Result
Definition: newton.hh:178
void setLinearSystemReuse(T &solver_backend, bool reuse, std::true_type)
Definition: newton.hh:41
Newton(const GridOperator &go, Solver &solver_)
Definition: newton.hh:772
Strategy
Definition: newton.hh:598
RFType first_defect
Definition: newton.hh:70
const GridOperator & gridoperator_
Definition: newton.hh:124
void setMinLinearReduction(RFType min_linear_reduction)
set the minimal reduction in the linear solver
Definition: newton.hh:509
const Entity & e
Definition: localfunctionspace.hh:111
typename impl::BackendMatrixSelector< Backend, VU, VV, E >::Type Matrix
alias of the return type of BackendMatrixSelector
Definition: backend/interface.hh:127
NewtonLineSearch(const GridOperator &go)
Definition: newton.hh:615
Definition: solver.hh:30
virtual void line_search(TrialVector &z, TestVector &r)
Definition: newton.hh:642
const std::string s
Definition: function.hh:830
Newton(const GridOperator &go, TrialVector &u_, Solver &solver_)
Definition: newton.hh:765
NewtonPrepareStep(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:489
void setVerbosityLevel(unsigned int verbosity_level)
Definition: newton.hh:96
Result res_
Definition: newton.hh:129
perform a linear search for the optimal damping parameter with multiples of damping ...
Definition: newton.hh:604
virtual ~NewtonBase()
Definition: newton.hh:158
NewtonSolver(const GridOperator &go, Solver &solver)
Definition: newton.hh:186
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: newton.hh:117
Definition: newton.hh:60
virtual void prepare_step(Matrix &A, TstV &)
Definition: newton.hh:535
void setLineSearchStrategy(std::string strategy)
Definition: newton.hh:627
const Result & result() const
Definition: newton.hh:196
NewtonTerminate(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:424
NewtonBase(const GridOperator &go, TrialVector &u)
Definition: newton.hh:138
Definition: newton.hh:68
void setReduction(RFType reduction)
Definition: newton.hh:442
void setMaxIterations(unsigned int maxit)
Definition: newton.hh:447
Definition: newton.hh:755
double assembler_time
Definition: newton.hh:72
TrialVector * u_
Definition: newton.hh:125
NewtonLineSearch(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:608
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
NewtonResult< RFType > Result
Definition: newton.hh:94
void setReassembleThreshold(RFType reassemble_threshold)
set a threshold, when the linear operator is reassembled
Definition: newton.hh:530
RFType linear_reduction_
Definition: newton.hh:132
RFType reduction_
Definition: newton.hh:134
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
NewtonResult()
Definition: newton.hh:76
std::shared_ptr< TrialVector > z_
Definition: newton.hh:126
NewtonPrepareStep(const GridOperator &go)
Definition: newton.hh:496
bool reassembled_
Definition: newton.hh:133
std::shared_ptr< TestVector > r_
Definition: newton.hh:127
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: newton.hh:111
unsigned int verbosity_level_
Definition: newton.hh:130
int linear_solver_iterations
Definition: newton.hh:74
std::shared_ptr< Matrix > A_
Definition: newton.hh:128
void setLineSearchStrategy(Strategy strategy)
Definition: newton.hh:622
Definition: newton.hh:167
RFType prev_defect_
Definition: newton.hh:131
Definition: newton.hh:64
Definition: newton.hh:61
RFType abs_limit_
Definition: newton.hh:135
STL namespace.
Definition: newton.hh:31
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: newton.hh:105
don&#39;t do any line search or damping
Definition: newton.hh:600
double linear_solver_time
Definition: newton.hh:73
Definition: newton.hh:480
virtual void defect(TestVector &r)
Definition: newton.hh:205
Strategy strategyFromName(const std::string &s)
Definition: newton.hh:738
NewtonBase(const GridOperator &go)
Definition: newton.hh:148
Definition: newton.hh:82
void setAbsoluteLimit(RFType abs_limit_)
Definition: newton.hh:457
void setForceIteration(bool force_iteration)
Definition: newton.hh:452
Definition: newton.hh:63
void setLineSearchDampingFactor(RFType damping_factor)
Definition: newton.hh:637
virtual bool terminate()
Definition: newton.hh:462
bool keep_matrix_
Definition: newton.hh:136
void setLineSearchMaxIterations(unsigned int maxit)
Definition: newton.hh:632
Definition: newton.hh:416
void setParameters(const Dune::ParameterTree &param)
interpret a parameter tree as a set of options for the newton solver
Definition: newton.hh:801
RFType defect
Definition: newton.hh:71
void setFixedLinearReduction(bool fixed_linear_reduction)
set a fixed reduction in the linear solver (overwrites setMinLinearReduction)
Definition: newton.hh:518
void apply()
Definition: newton.hh:252
NewtonSolver(const GridOperator &go, TrialVector &u_, Solver &solver)
Definition: newton.hh:180
Definition: newton.hh:589
NewtonTerminate(const GridOperator &go)
Definition: newton.hh:433