dune-fem  2.8-git
newtoninverseoperator.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_NEWTONINVERSEOPERATOR_HH
2 #define DUNE_FEM_NEWTONINVERSEOPERATOR_HH
3 
4 #include <cassert>
5 #include <cmath>
6 
7 #include <iostream>
8 #include <limits>
9 #include <memory>
10 #include <string>
11 #include <utility>
12 
14 #include <dune/fem/io/parameter.hh>
17 
18 namespace Dune
19 {
20 
21  namespace Fem
22  {
23 
24  // NewtonParameter
25  // ---------------
26 
27  template <class SolverParam = SolverParameter>
29  : public Dune::Fem::LocalParameter< NewtonParameter<SolverParam>, NewtonParameter<SolverParam> >
30  {
31  protected:
32 
33  std::shared_ptr<SolverParam> baseParam_;
34  // key prefix, default is fem.solver.newton. (can be overloaded by user)
35  const std::string keyPrefix_;
36 
38 
39  public:
40  NewtonParameter( const SolverParam& baseParameter, const std::string keyPrefix = "fem.solver.newton." )
41  : baseParam_( static_cast< SolverParam* > (baseParameter.clone()) ),
42  keyPrefix_( keyPrefix ),
43  parameter_( baseParameter.parameter() )
44  {}
45 
46  template <class Parameter, std::enable_if_t<!std::is_base_of<SolverParam,Parameter>::value && !std::is_same<Parameter,ParameterReader>::value,int> i=0>
47  NewtonParameter( const Parameter& solverParameter, const std::string keyPrefix = "fem.solver.newton." )
48  : baseParam_( new SolverParam(solverParameter) ),
49  keyPrefix_( keyPrefix ),
51  {}
52 
53  template <class ParamReader, std::enable_if_t<!std::is_same<ParamReader,SolverParam>::value && std::is_same<ParamReader,ParameterReader>::value,int> i=0>
54  NewtonParameter( const ParamReader &parameter, const std::string keyPrefix = "fem.solver.newton." )
55  : baseParam_( std::make_shared<SolverParam>( keyPrefix + "linear.", parameter) ),
56  keyPrefix_( keyPrefix),
58  {}
59 
60  const ParameterReader &parameter () const { return parameter_; }
61  const SolverParam& solverParameter () const { return *baseParam_; }
62  const SolverParam& linear () const { return *baseParam_; }
63 
64  virtual void reset ()
65  {
66  baseParam_->reset();
67  tolerance_ = -1;
68  verbose_ = -1;
69  maxIterations_ = -1;
70  maxLinearIterations_ = -1;
71  maxLineSearchIterations_ = -1;
72  }
73 
74  //These methods affect the nonlinear solver
75  virtual double tolerance () const
76  {
77  if(tolerance_ < 0)
78  tolerance_ = parameter_.getValue< double >( keyPrefix_ + "tolerance", 1e-6 );
79  return tolerance_;
80  }
81 
82  virtual void setTolerance ( const double tol )
83  {
84  assert( tol > 0 );
85  tolerance_ = tol;
86  }
87 
88  virtual bool verbose () const
89  {
90  if(verbose_ < 0)
91  {
92  // the following causes problems with different default values
93  // used if baseParam_->keyPrefix is not default but the default is
94  // also used in the program
95  // const bool v = baseParam_? baseParam_->verbose() : false;
96  const bool v = false;
97  verbose_ = parameter_.getValue< bool >(keyPrefix_ + "verbose", v ) ? 1 : 0 ;
98  }
99  return verbose_;
100  }
101 
102  virtual void setVerbose( bool verb)
103  {
104  verbose_ = verb ? 1 : 0;
105  }
106 
107  virtual int maxIterations () const
108  {
109  if(maxIterations_ < 0)
110  maxIterations_ = parameter_.getValue< int >( keyPrefix_ + "maxiterations", std::numeric_limits< int >::max() );
111  return maxIterations_;
112  }
113 
114  virtual void setMaxIterations ( const int maxIter )
115  {
116  assert(maxIter >= 0);
117  maxIterations_ = maxIter;
118  }
119 
120  //Maximum Linear Iterations in total
122  virtual int maxLinearIterations () const
123  {
124  if(maxLinearIterations_ < 0)
125  maxLinearIterations_ = linear().maxIterations();
126  return maxLinearIterations_;
127  }
128 
129  virtual void setMaxLinearIterations ( const int maxLinearIter )
130  {
131  assert(maxLinearIter >=0);
132  maxLinearIterations_ = maxLinearIter;
133  }
134 
135  virtual int maxLineSearchIterations () const
136  {
137  if(maxLineSearchIterations_ < 0)
138  maxLineSearchIterations_ = parameter_.getValue< int >( keyPrefix_ + "maxlinesearchiterations", std::numeric_limits< int >::max() );
139  return maxLineSearchIterations_;
140  }
141 
142  virtual void setMaxLineSearchIterations ( const int maxLineSearchIter )
143  {
144  assert( maxLineSearchIter >= 0);
145  maxLineSearchIterations_ = maxLineSearchIter;
146  }
147 
148  enum class LineSearchMethod {
149  none = 0,
150  simple = 1
151  };
152 
153  virtual LineSearchMethod lineSearch () const
154  {
155  const std::string lineSearchMethods[] = { "none", "simple" };
156  return static_cast< LineSearchMethod>( parameter_.getEnum( keyPrefix_ + "lineSearch", lineSearchMethods, 0 ) );
157  }
158 
159  virtual void setLineSearch ( const LineSearchMethod method )
160  {
161  const std::string lineSearchMethods[] = { "none", "simple" };
162  Parameter::append( keyPrefix_ + "lineSearch", lineSearchMethods[int(method)], true );
163  }
164 
165  private:
166  mutable double tolerance_ = -1;
167  mutable int verbose_ = -1;
168  mutable int maxIterations_ = -1;
169  mutable int maxLinearIterations_ = -1;
170  mutable int maxLineSearchIterations_ = -1;
171  };
172 
173 
174 
175  // NewtonFailure
176  // -------------
177 
178  enum class NewtonFailure
179  // : public int
180  {
181  Success = 0,
182  InvalidResidual = 1,
183  IterationsExceeded = 2,
185  LineSearchFailed = 4,
186  TooManyIterations = 5,
189  };
190 
191 
192 
193  // NewtonInverseOperator
194  // ---------------------
195 
206  template< class JacobianOperator, class LInvOp >
208  : public Operator< typename JacobianOperator::RangeFunctionType, typename JacobianOperator::DomainFunctionType >
209  {
212 
213  public:
215  typedef JacobianOperator JacobianOperatorType;
216 
219 
222 
225 
227 
229 
230  typedef std::function< bool ( const RangeFunctionType &w, const RangeFunctionType &dw, double residualNorm ) > ErrorMeasureType;
231 
249  // main constructor
251  : verbose_( parameter.verbose() && Parameter::verbose() ),
252  maxLineSearchIterations_( parameter.maxLineSearchIterations() ),
253  jInv_( std::move( jInv ) ),
254  parameter_(parameter),
255  lsMethod_( parameter.lineSearch() ),
256  finished_( [ epsilon ] ( const RangeFunctionType &w, const RangeFunctionType &dw, double res ) { return res < epsilon; } )
257  {}
258 
259 
265  /*
266  explicit NewtonInverseOperator ( const ParameterType &parameter )
267  : NewtonInverseOperator( parameter.tolerance(), parameter )
268  {}
269  */
270 
272  : NewtonInverseOperator( parameter.tolerance(), parameter )
273  {
274  // std::cout << "in Newton inv op should use:" << parameter.linear().solverMethod({SolverParameter::gmres,SolverParameter::bicgstab,SolverParameter::minres}) << std::endl;
275  }
276 
281  NewtonInverseOperator ( const DomainFieldType &epsilon, const ParameterType &parameter )
283  LinearInverseOperatorType( parameter.linear() ),
284  epsilon, parameter )
285  {}
286 
288  const ParameterReader &parameter = Parameter::container() )
289  : NewtonInverseOperator( epsilon, ParameterType( parameter ) )
290  {}
291 
292 
301  void setErrorMeasure ( ErrorMeasureType finished ) { finished_ = std::move( finished ); }
302 
303  void bind ( const OperatorType &op ) { op_ = &op; }
304 
305  void unbind () { op_ = nullptr; }
306 
307  virtual void operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const;
308 
309  int iterations () const { return iterations_; }
310  void setMaxIterations ( int maxIterations ) { parameter_.setMaxIterations( maxIterations ); }
311  int linearIterations () const { return linearIterations_; }
312  void setMaxLinearIterations ( int maxLinearIterations ) { parameter_.setMaxLinearIterations( maxLinearIterations ); }
313  bool verbose() const { return verbose_; }
314 
316  {
317  // check for finite |residual| - this also works for -ffinite-math-only (gcc)
318  // nan test not working with optimization flags...
319  if( !(delta_ < std::numeric_limits< DomainFieldType >::max()) || std::isnan( delta_ ) )
321  else if( iterations_ >= parameter_.maxIterations() )
323  else if( linearIterations_ >= parameter_.maxLinearIterations() )
325  else if( linearIterations_ < 0)
327  else if( !stepCompleted_ )
329  else
330  return NewtonFailure::Success;
331  }
332 
333  bool converged () const { return failed() == NewtonFailure::Success; }
334 
336  const DomainFunctionType &u, DomainFunctionType &residual) const
337  {
338  double deltaOld = delta_;
339  delta_ = std::sqrt( residual.scalarProductDofs( residual ) );
340  if (lsMethod_ == ParameterType::LineSearchMethod::none)
341  return 0;
343  {
344  double test = dw.scalarProductDofs( dw );
346  !std::isnan(test)) )
347  delta_ = 2.*deltaOld; // enter line search
348  }
349  double factor = 1.0;
350  int noLineSearch = (delta_ < deltaOld)?1:0;
351  int lineSearchIteration = 0;
352  while (delta_ >= deltaOld)
353  {
354  double deltaPrev = delta_;
355  factor *= 0.5;
356  if( verbose() )
357  std::cerr << " line search:" << delta_ << ">" << deltaOld << std::endl;
358  if (std::abs(delta_-deltaOld) < 1e-5*delta_) // || !converged()) // line search not working
359  return -1; // failed
360  w.axpy(factor,dw);
361  (*op_)( w, residual );
362  residual -= u;
363  delta_ = std::sqrt( residual.scalarProductDofs( residual ) );
364  if (std::abs(delta_-deltaPrev) < 1e-15)
365  return -1;
367  delta_ = 2.*deltaOld; // remain in line search
368 
369  ++lineSearchIteration;
370  if( lineSearchIteration >= maxLineSearchIterations_ )
371  return -1; // failed
372  }
373  return noLineSearch;
374  }
375 
376  protected:
377  // hold pointer to jacobian operator, if memory reallocation is needed, the operator should know how to handle this.
378  template< class ... Args>
379  JacobianOperatorType& jacobian ( Args && ... args ) const
380  {
381  if( !jOp_ )
382  jOp_.reset( new JacobianOperatorType( std::forward< Args >( args ) ...) ); //, parameter_.parameter() ) );
383  return *jOp_;
384  }
385 
386  private:
387  const OperatorType *op_ = nullptr;
388 
389  const bool verbose_;
390  const int maxLineSearchIterations_;
391 
392  mutable DomainFieldType delta_;
393  mutable int iterations_;
394  mutable int linearIterations_;
395  mutable LinearInverseOperatorType jInv_;
396  mutable std::unique_ptr< JacobianOperatorType > jOp_;
397  ParameterType parameter_;
398  mutable int stepCompleted_;
399  typename ParameterType::LineSearchMethod lsMethod_;
400  ErrorMeasureType finished_;
401  };
402 
403 
404  // Implementation of NewtonInverseOperator
405  // ---------------------------------------
406 
407  template< class JacobianOperator, class LInvOp >
410  {
411  assert( op_ );
412 
413  DomainFunctionType residual( u );
414  RangeFunctionType dw( w );
415  JacobianOperatorType& jOp = jacobian( "jacobianOperator", dw.space(), u.space()); // , parameter_.parameter() );
416 
417  stepCompleted_ = true;
418  iterations_ = 0;
419  linearIterations_ = 0;
420  // compute initial residual
421  (*op_)( w, residual );
422  residual -= u;
423  delta_ = std::sqrt( residual.scalarProductDofs( residual ) );
424 
425  if( verbose() )
426  std::cerr << "Newton iteration " << iterations_ << ": |residual| = " << delta_;
427  while( true )
428  {
429  if( verbose() )
430  std::cerr << std::endl;
431  // evaluate operator's jacobian
432  (*op_).jacobian( w, jOp );
433 
434  // David: With this factor, the tolerance of CGInverseOp is the absolute
435  // rather than the relative error
436  // (see also dune-fem/dune/fem/solver/krylovinverseoperators.hh)
437  jInv_.bind( jOp );
438  if ( parameter_.maxLinearIterations() - linearIterations_ <= 0 )
439  break;
440  jInv_.setMaxIterations( parameter_.maxLinearIterations() - linearIterations_ );
441 
442  dw.clear();
443  jInv_( residual, dw );
444  if (jInv_.iterations() < 0) // iterations are negative if solver didn't converge
445  {
446  linearIterations_ = jInv_.iterations();
447  break;
448  }
449  linearIterations_ += jInv_.iterations();
450  w -= dw;
451 
452  (*op_)( w, residual );
453  residual -= u;
454  int ls = lineSearch(w,dw,u,residual);
455  stepCompleted_ = ls >= 0;
456  ++iterations_;
457  if( verbose() )
458  std::cerr << "Newton iteration " << iterations_ << ": |residual| = " << delta_ << std::flush;
459  // if ( (ls==1 && finished_(w, dw, delta_)) || !converged())
460  if ( (finished_(w, dw, delta_)) || !converged())
461  {
462  if( verbose() )
463  std::cerr << std::endl;
464  break;
465  }
466  }
467  if( verbose() )
468  std::cerr << std::endl;
469 
470  jInv_.unbind();
471  }
472 
473  } // namespace Fem
474 
475 } // namespace Dune
476 
477 #endif // #ifndef DUNE_FEM_NEWTONINVERSEOPERATOR_HH
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:871
static double sqrt(const Double &v)
Definition: double.hh:886
NewtonFailure
Definition: newtoninverseoperator.hh:180
static constexpr T max(T a)
Definition: utility.hh:77
int getEnum(const std::string &key, const std::string(&values)[n]) const
Definition: reader.hh:225
T getValue(const std::string &key) const
get mandatory parameter
Definition: reader.hh:159
Container for User Specified Parameters.
Definition: io/parameter.hh:191
static void append(int &argc, char **argv)
add parameters from the command line RangeType gRight;
Definition: io/parameter.hh:208
static ParameterContainer & container()
Definition: io/parameter.hh:193
Definition: io/parameter.hh:551
abstract differentiable operator
Definition: differentiableoperator.hh:29
abstract operator
Definition: operator.hh:34
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
DomainFunction::RangeFieldType DomainFieldType
field type of the operator's domain
Definition: operator.hh:41
RangeFunction RangeFunctionType
type of discrete function in the operator's range
Definition: operator.hh:38
Definition: newtoninverseoperator.hh:30
virtual void setVerbose(bool verb)
Definition: newtoninverseoperator.hh:102
const SolverParam & linear() const
Definition: newtoninverseoperator.hh:62
virtual void reset()
Definition: newtoninverseoperator.hh:64
virtual int maxLinearIterations() const
= max iterations of each linear solve
Definition: newtoninverseoperator.hh:122
virtual int maxLineSearchIterations() const
Definition: newtoninverseoperator.hh:135
virtual void setMaxLineSearchIterations(const int maxLineSearchIter)
Definition: newtoninverseoperator.hh:142
const SolverParam & solverParameter() const
Definition: newtoninverseoperator.hh:61
NewtonParameter(const SolverParam &baseParameter, const std::string keyPrefix="fem.solver.newton.")
Definition: newtoninverseoperator.hh:40
LineSearchMethod
Definition: newtoninverseoperator.hh:148
virtual void setMaxLinearIterations(const int maxLinearIter)
Definition: newtoninverseoperator.hh:129
virtual int maxIterations() const
Definition: newtoninverseoperator.hh:107
const ParameterReader & parameter() const
Definition: newtoninverseoperator.hh:60
virtual void setLineSearch(const LineSearchMethod method)
Definition: newtoninverseoperator.hh:159
std::shared_ptr< SolverParam > baseParam_
Definition: newtoninverseoperator.hh:33
virtual double tolerance() const
Definition: newtoninverseoperator.hh:75
NewtonParameter(const ParamReader &parameter, const std::string keyPrefix="fem.solver.newton.")
Definition: newtoninverseoperator.hh:54
virtual void setTolerance(const double tol)
Definition: newtoninverseoperator.hh:82
const std::string keyPrefix_
Definition: newtoninverseoperator.hh:35
virtual LineSearchMethod lineSearch() const
Definition: newtoninverseoperator.hh:153
NewtonParameter(const Parameter &solverParameter, const std::string keyPrefix="fem.solver.newton.")
Definition: newtoninverseoperator.hh:47
virtual void setMaxIterations(const int maxIter)
Definition: newtoninverseoperator.hh:114
ParameterReader parameter_
Definition: newtoninverseoperator.hh:37
virtual bool verbose() const
Definition: newtoninverseoperator.hh:88
inverse operator based on a newton scheme
Definition: newtoninverseoperator.hh:209
std::function< bool(const RangeFunctionType &w, const RangeFunctionType &dw, double residualNorm) > ErrorMeasureType
Definition: newtoninverseoperator.hh:230
NewtonInverseOperator(LinearInverseOperatorType jInv, const DomainFieldType &epsilon, const ParameterType &parameter)
Definition: newtoninverseoperator.hh:250
void unbind()
Definition: newtoninverseoperator.hh:305
virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const
Definition: newtoninverseoperator.hh:409
int linearIterations() const
Definition: newtoninverseoperator.hh:311
JacobianOperator JacobianOperatorType
type of operator's Jacobian
Definition: newtoninverseoperator.hh:215
bool converged() const
Definition: newtoninverseoperator.hh:333
void setMaxIterations(int maxIterations)
Definition: newtoninverseoperator.hh:310
BaseType::DomainFieldType DomainFieldType
Definition: newtoninverseoperator.hh:226
NewtonInverseOperator(const ParameterType &parameter=ParameterType(Parameter::container()))
Definition: newtoninverseoperator.hh:271
JacobianOperatorType & jacobian(Args &&... args) const
Definition: newtoninverseoperator.hh:379
int iterations() const
Definition: newtoninverseoperator.hh:309
BaseType::DomainFunctionType DomainFunctionType
Definition: newtoninverseoperator.hh:223
void bind(const OperatorType &op)
Definition: newtoninverseoperator.hh:303
NewtonInverseOperator(const DomainFieldType &epsilon, const ParameterReader &parameter=Parameter::container())
Definition: newtoninverseoperator.hh:287
BaseType::RangeFunctionType RangeFunctionType
Definition: newtoninverseoperator.hh:224
void setErrorMeasure(ErrorMeasureType finished)
Definition: newtoninverseoperator.hh:301
DifferentiableOperator< JacobianOperatorType > OperatorType
type of operator to invert
Definition: newtoninverseoperator.hh:218
void setMaxLinearIterations(int maxLinearIterations)
Definition: newtoninverseoperator.hh:312
bool verbose() const
Definition: newtoninverseoperator.hh:313
LInvOp LinearInverseOperatorType
type of linear inverse operator
Definition: newtoninverseoperator.hh:221
NewtonInverseOperator(const DomainFieldType &epsilon, const ParameterType &parameter)
Definition: newtoninverseoperator.hh:281
NewtonFailure failed() const
Definition: newtoninverseoperator.hh:315
virtual int lineSearch(RangeFunctionType &w, RangeFunctionType &dw, const DomainFunctionType &u, DomainFunctionType &residual) const
Definition: newtoninverseoperator.hh:335
NewtonParameter< typename LinearInverseOperatorType::SolverParameterType > ParameterType
Definition: newtoninverseoperator.hh:228