dune-fem 2.8.0
Loading...
Searching...
No Matches
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
17
18namespace 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
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,
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
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 );
345 if (! (test < std::numeric_limits< DomainFieldType >::max() &&
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 >
408 inline void NewtonInverseOperator< JacobianOperator, LInvOp >
409 ::operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const
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
STL namespace.
double sqrt(const Dune::Fem::Double &v)
Definition: double.hh:977
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
NewtonFailure
Definition: newtoninverseoperator.hh:180
Container for User Specified Parameters.
Definition: io/parameter.hh:191
static ParameterContainer & container()
Definition: io/parameter.hh:193
static void append(int &argc, char **argv)
add parameters from the command line RangeType gRight;
Definition: io/parameter.hh:208
Definition: io/parameter.hh:551
virtual ParamDefault * clone() const
Definition: io/parameter.hh:555
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
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
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
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 SolverParam & linear() const
Definition: newtoninverseoperator.hh:62
const ParameterReader & parameter() const
Definition: newtoninverseoperator.hh:60
const SolverParam & solverParameter() const
Definition: newtoninverseoperator.hh:61
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
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
JacobianOperatorType & jacobian(Args &&... args) const
Definition: newtoninverseoperator.hh:379
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