dune-pdelab  2.7-git
linesearch.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_LINESEARCH_HH
4 #define DUNE_PDELAB_SOLVER_LINESEARCH_HH
5 
6 
7 namespace Dune::PDELab
8 {
9 
10  class LineSearchError : public Exception {};
11 
13  template <typename Domain>
15  {
16  public:
18  virtual ~LineSearchInterface () {}
19 
21  virtual void lineSearch(Domain&, const Domain&) = 0;
22 
24  virtual void setParameters(const ParameterTree&) = 0;
25 
27  virtual void printParameters() const
28  {
29  std::cout << "NewtonMethod::_lineSearch->printParameters() is not implemented." << std::endl;
30  }
31  };
32 
33 
35  template <typename Solver>
36  class LineSearchNone : public LineSearchInterface<typename Solver::Domain>
37  {
38  public:
39  using Domain = typename Solver::Domain;
40  using Real = typename Solver::Real;
41 
42  LineSearchNone(Solver& solver) : _solver(solver) {}
43 
45  virtual void lineSearch(Domain& solution, const Domain& correction) override
46  {
47  solution.axpy(-1.0, correction);
48  _solver.updateDefect(solution);
49  }
50 
51  virtual void setParameters(const ParameterTree&) override {}
52 
53  // print line search type
54  virtual void printParameters() const override
55  {
56  std::cout << "LineSearch.Type........... None" << std::endl;
57  }
58 
59  private:
60  Solver& _solver;
61  };
62 
63 
70  template <typename Solver>
71  class LineSearchHackbuschReusken : public LineSearchInterface<typename Solver::Domain>
72  {
73  public:
74  using Domain = typename Solver::Domain;
75  using Real = typename Solver::Real;
76 
77  LineSearchHackbuschReusken(Solver& solver, bool forceAcceptBest = false) :
78  _solver(solver), _forceAcceptBest(forceAcceptBest) {}
79 
81  virtual void lineSearch(Domain& solution, const Domain& correction) override
82  {
83  if ((_solver.result().defect < _solver.getAbsoluteLimit())){
84  solution.axpy(-1.0, correction);
85  _solver.updateDefect(solution);
86  return;
87  }
88 
89  auto verbosity = _solver.getVerbosityLevel();
90 
91  if (verbosity >= 4)
92  std::cout << " Performing line search..." << std::endl;
93  Real lambda = 1.0;
94  Real bestLambda = 0.0;
95  Real bestDefect = _solver.result().defect;
96  Real previousDefect = _solver.result().defect;
97  bool converged = false;
98 
99  if (not _previousSolution)
100  _previousSolution = std::make_shared<Domain>(solution);
101  else
102  *_previousSolution = solution;
103 
104  for (unsigned int iteration = 0; iteration < _lineSearchMaxIterations; ++iteration){
105  if (verbosity >= 4)
106  std::cout << " trying line search damping factor: "
107  << std::setw(12) << std::setprecision(4) << std::scientific
108  << lambda
109  << std::endl;
110 
111  solution.axpy(-lambda, correction);
112  _solver.updateDefect(solution);
113  if (verbosity >= 4){
114  if (not std::isfinite(_solver.result().defect))
115  std::cout << " NaNs detected" << std::endl;
116  } // ignore NaNs and try again with lower lambda
117 
118  if (_solver.result().defect <= (1.0 - lambda/4) * previousDefect){
119  if (verbosity >= 4)
120  std::cout << " line search converged" << std::endl;
121  converged = true;
122  break;
123  }
124 
125  if (_solver.result().defect < bestDefect){
126  bestDefect = _solver.result().defect;
127  bestLambda = lambda;
128  }
129 
130  lambda *= _lineSearchDampingFactor;
131  solution = *_previousSolution;
132  }
133 
134  if (not converged){
135  if (verbosity >= 4)
136  std::cout << " max line search iterations exceeded" << std::endl;
137 
138  if (not (_acceptBest or _forceAcceptBest)){
139  solution = *_previousSolution;
140  _solver.updateDefect(solution);
141  DUNE_THROW( LineSearchError,
142  "LineSearch::lineSearch(): line search failed, "
143  "max iteration count reached, "
144  "defect did not improve enough");
145  }
146  else{
147  if (bestLambda == 0.0){
148  solution = *_previousSolution;
149  _solver.updateDefect(solution);
150  DUNE_THROW(LineSearchError,
151  "LineSearch::lineSearch(): line search failed, "
152  "max iteration count reached, "
153  "defect did not improve in any of the iterations");
154  }
155  if (bestLambda != lambda){
156  solution = *_previousSolution;
157  solution.axpy(-bestLambda, correction);
158  _solver.updateDefect(solution);
159  converged = true;
160  }
161  }
162  }
163 
164  if (converged)
165  if (verbosity >= 4)
166  std::cout << " line search damping factor: "
167  << std::setw(12) << std::setprecision(4) << std::scientific
168  << lambda << std::endl;
169  }
170 
171 
172  /* \brief Set parameters
173  *
174  * Possible parameters are:
175  *
176  * - MaxIterations: Maximum number of line search iterations.
177  *
178  * - DampingFactor: Multiplier to line search parameter after each iteration.
179  *
180  * - AcceptBest: Accept the best line search parameter if
181  * there was any improvement, even if the convergence criterion was not
182  * reached.
183  */
184  virtual void setParameters(const ParameterTree& parameterTree) override
185  {
186  _lineSearchMaxIterations = parameterTree.get<unsigned int>("MaxIterations",
187  _lineSearchMaxIterations);
188  _lineSearchDampingFactor = parameterTree.get<Real>("DampingFactor",
189  _lineSearchDampingFactor);
190  _acceptBest = parameterTree.get<bool>("AcceptBest", _acceptBest);
191  }
192 
193  virtual void printParameters() const override
194  {
195  std::cout << "LineSearch.Type........... Hackbusch-Reusken" << std::endl;
196  std::cout << "LineSearch.MaxIterations.. " << _lineSearchMaxIterations << std::endl;
197  std::cout << "LineSearch.DampingFactor.. " << _lineSearchDampingFactor << std::endl;
198  std::cout << "LineSearch.AcceptBest..... " << (_acceptBest or _forceAcceptBest) << std::endl;
199  }
200 
201  private:
202  Solver& _solver;
203  std::shared_ptr<Domain> _previousSolution;
204 
205  // Line search parameters
206  unsigned int _lineSearchMaxIterations = 10;
207  Real _lineSearchDampingFactor = 0.5;
208  bool _acceptBest = false;
209  bool _forceAcceptBest;
210  };
211 
214  {
215  noLineSearch,
218  };
219 
220  // we put this into an emty namespace, so that we don't violate the one-definition-rule
221  namespace {
228  inline
229  LineSearchStrategy lineSearchStrategyFromString (const std::string& name)
230  {
231  if (name == "noLineSearch")
233  if (name == "hackbuschReusken")
235  if (name == "hackbuschReuskenAcceptBest")
237  DUNE_THROW(Exception,"Unkown line search strategy: " << name);
238  }
239  }
240 
241 
252  template <typename Solver>
253  std::shared_ptr<LineSearchInterface<typename Solver::Domain>>
254  createLineSearch(Solver& solver, LineSearchStrategy strategy)
255  {
256  if (strategy == LineSearchStrategy::noLineSearch){
257  auto lineSearch = std::make_shared<LineSearchNone<Solver>> (solver);
258  return lineSearch;
259  }
260  if (strategy == LineSearchStrategy::hackbuschReusken){
261  auto lineSearch = std::make_shared<LineSearchHackbuschReusken<Solver>> (solver);
262  return lineSearch;
263  }
265  auto lineSearch = std::make_shared<LineSearchHackbuschReusken<Solver>> (solver, true);
266  std::cout << "Warning: linesearch hackbuschReuskenAcceptBest is deprecated and will be removed after PDELab 2.7.\n"
267  << " Please use 'hackbuschReusken' and add the parameter 'LineSearchAcceptBest : true'";
268  return lineSearch;
269  }
270  DUNE_THROW(Exception,"Unkown line search strategy");
271  }
272 
273 } // namespace Dune::PDELab
274 
275 #endif
Definition: adaptivity.hh:29
std::shared_ptr< LineSearchInterface< typename Solver::Domain > > createLineSearch(Solver &solver, LineSearchStrategy strategy)
fectory function to create an instace of a line-search
Definition: linesearch.hh:254
LineSearchStrategy
Flags for different line search strategies.
Definition: linesearch.hh:214
Base class for all PDELab exceptions.
Definition: exceptions.hh:19
Definition: linesearch.hh:10
Abstract base class describing the line search interface.
Definition: linesearch.hh:15
virtual ~LineSearchInterface()
Every abstract base class should have a virtual destructor.
Definition: linesearch.hh:18
virtual void printParameters() const
Print paramters.
Definition: linesearch.hh:27
virtual void setParameters(const ParameterTree &)=0
Set parameters.
virtual void lineSearch(Domain &, const Domain &)=0
Do line search.
Class for simply updating the solution without line search.
Definition: linesearch.hh:37
virtual void printParameters() const override
Print paramters.
Definition: linesearch.hh:54
typename Solver::Domain Domain
Definition: linesearch.hh:39
LineSearchNone(Solver &solver)
Definition: linesearch.hh:42
virtual void lineSearch(Domain &solution, const Domain &correction) override
Do line search (in this case just update the solution)
Definition: linesearch.hh:45
typename Solver::Real Real
Definition: linesearch.hh:40
virtual void setParameters(const ParameterTree &) override
Set parameters.
Definition: linesearch.hh:51
Hackbusch-Reusken line search.
Definition: linesearch.hh:72
typename Solver::Real Real
Definition: linesearch.hh:75
virtual void setParameters(const ParameterTree &parameterTree) override
Set parameters.
Definition: linesearch.hh:184
virtual void printParameters() const override
Print paramters.
Definition: linesearch.hh:193
virtual void lineSearch(Domain &solution, const Domain &correction) override
Do line search.
Definition: linesearch.hh:81
LineSearchHackbuschReusken(Solver &solver, bool forceAcceptBest=false)
Definition: linesearch.hh:77
typename Solver::Domain Domain
Definition: linesearch.hh:74