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