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