DUNE-FEM (unstable)

newtoninverseoperator.hh
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 <regex>
12#include <utility>
13
14#include <dune/common/timer.hh>
16
17#include <dune/fem/common/staticlistofint.hh>
18#include <dune/fem/solver/parameter.hh>
19#include <dune/fem/io/parameter.hh>
20#include <dune/fem/operator/common/operator.hh>
21#include <dune/fem/operator/common/differentiableoperator.hh>
22#include <dune/fem/solver/inverseoperatorinterface.hh>
23
24namespace Dune
25{
26
27 namespace Fem
28 {
29
38 {
39 protected:
40 const double etaMax_ = 0.99;
41 const double gamma_ = 0.1;
42 mutable double previousEta_ = -1.0;
43 mutable double previousResidual_ = -1.0;
44 mutable double newtonTolerance_;
45
46 public:
50 EisenstatWalkerStrategy(const double newtonTolerance) : newtonTolerance_(newtonTolerance) {}
51 double nextLinearTolerance(const double currentResidual) const
52 {
53 double eta = etaMax_;
54 // First call previousEta_ is negative
55 if (previousEta_ >= 0.0)
56 {
57 const double etaA = gamma_ * currentResidual * currentResidual / (previousResidual_ * previousResidual_);
58 const double indicator = gamma_ * previousEta_ * previousEta_;
59 const double etaC = indicator < 0.1 ? std::min(etaA, etaMax_) : std::min(etaMax_, std::max(etaA, indicator));
60 eta = std::min(etaMax_, std::max(etaC, 0.5 * newtonTolerance_ / currentResidual));
61 }
62 previousResidual_ = currentResidual;
63 previousEta_ = eta;
64 return eta;
65 }
66 void setTolerance(const double newtonTolerance) { newtonTolerance_ = newtonTolerance; }
67 };
68
69 // NewtonParameter
70 // ---------------
71
72 template <class SolverParam = SolverParameter>
73 struct NewtonParameter
74 : public Dune::Fem::LocalParameter< NewtonParameter<SolverParam>, NewtonParameter<SolverParam> >
75 {
76 protected:
77
78 std::shared_ptr<SolverParam> baseParam_;
79 // key prefix, default is fem.solver.nonlinear. (can be overloaded by user)
80 const std::string keyPrefix_;
81
82 ParameterReader parameter_;
83
84 void checkDeprecatedParameters() const
85 {
86 const std::string newton("newton.");
87 const std::size_t pos = keyPrefix_.find( newton );
88 if( pos != std::string::npos )
89 {
90 DUNE_THROW(InvalidStateException,"Keyprefix 'newton' is deprecated, replace with 'nonlinear'!");
91 }
92
93 const std::string params[]
94 = { "tolerance", "lineSearch", "maxterations", "linear", "maxlinesearchiterations" };
95 for( const auto& p : params )
96 {
97 std::string key( "fem.solver.newton." );
98 key += p;
99 if( parameter_.exists( key ) )
100 DUNE_THROW(InvalidStateException,"Keyprefix 'newton' is deprecated, replace with 'nonlinear'!");
101 }
102 }
103
104 std::string replaceNonLinearWithLinear( const std::string& keyPrefix ) const
105 {
106 if( keyPrefix.find( "nonlinear" ) != std::string::npos )
107 {
108 return std::regex_replace(keyPrefix, std::regex("nonlinear"), "linear" );
109 }
110 else
111 return keyPrefix;
112 }
113 public:
114 NewtonParameter( const SolverParam& baseParameter, const std::string keyPrefix = "fem.solver.nonlinear." )
115 : baseParam_( static_cast< SolverParam* > (baseParameter.clone()) ),
116 keyPrefix_( keyPrefix ),
117 parameter_( baseParameter.parameter() )
118 {
119 checkDeprecatedParameters();
120 checkForcingErrorMeasure();
121 }
122
123 template <class Parameter, std::enable_if_t<!std::is_base_of<SolverParam,Parameter>::value && !std::is_same<Parameter,ParameterReader>::value,int> i=0>
124 NewtonParameter( const Parameter& solverParameter, const std::string keyPrefix = "fem.solver.nonlinear." )
125 : baseParam_( new SolverParam(solverParameter) ),
126 keyPrefix_( keyPrefix ),
127 parameter_( solverParameter.parameter() )
128 {
129 checkDeprecatedParameters();
130 checkForcingErrorMeasure();
131 }
132
133 template <class ParamReader, std::enable_if_t<!std::is_same<ParamReader,SolverParam>::value && std::is_same<ParamReader,ParameterReader>::value,int> i=0>
134 NewtonParameter( const ParamReader &parameter, const std::string keyPrefix = "fem.solver.nonlinear." )
135 // pass keyprefix for linear solvers, which is the same as keyprefix with nonlinear replaced by linear
136 : baseParam_( std::make_shared<SolverParam>( replaceNonLinearWithLinear(keyPrefix), parameter) ),
137 keyPrefix_( keyPrefix ),
138 parameter_( parameter )
139 {
140 checkDeprecatedParameters();
141 checkForcingErrorMeasure();
142 }
143
144 void checkForcingErrorMeasure()
145 {
146 if (forcing() == Forcing::eisenstatwalker)
147 {
148 baseParam_->setDefaultErrorMeasure(2);
149 if (baseParam_->errorMeasure() != LinearSolver::ToleranceCriteria::residualReduction)
150 DUNE_THROW( InvalidStateException, "Parameter `linear.errormeasure` selecting the tolerance criteria in the linear solver must be `residualreduction` when using Eisenstat-Walker." );
151 }
152 }
153
154 const ParameterReader &parameter () const { return parameter_; }
155 const SolverParam& solverParameter () const { return *baseParam_; }
156 const SolverParam& linear () const { return *baseParam_; }
157
158 virtual void reset ()
159 {
160 baseParam_->reset();
161 tolerance_ = -1;
162 verbose_ = -1;
163 maxIterations_ = -1;
164 maxLinearIterations_ = -1;
165 maxLineSearchIterations_ = -1;
166 }
167
168 //These methods affect the nonlinear solver
169 virtual double tolerance () const
170 {
171 if(tolerance_ < 0)
172 tolerance_ = parameter_.getValue< double >( keyPrefix_ + "tolerance", 1e-6 );
173 return tolerance_;
174 }
175
176 virtual void setTolerance ( const double tol )
177 {
178 assert( tol > 0 );
179 tolerance_ = tol;
180 }
181
182 virtual bool verbose () const
183 {
184 if(verbose_ < 0)
185 {
186 // the following causes problems with different default values
187 // used if baseParam_->keyPrefix is not default but the default is
188 // also used in the program
189 // const bool v = baseParam_? baseParam_->verbose() : false;
190 const bool v = false;
191 verbose_ = parameter_.getValue< bool >(keyPrefix_ + "verbose", v ) ? 1 : 0 ;
192 }
193 return verbose_;
194 }
195
196 virtual void setVerbose( bool verb)
197 {
198 verbose_ = verb ? 1 : 0;
199 }
200
201 virtual int maxIterations () const
202 {
203 if(maxIterations_ < 0)
204 maxIterations_ = parameter_.getValue< int >( keyPrefix_ + "maxiterations", std::numeric_limits< int >::max() );
205 return maxIterations_;
206 }
207
208 virtual void setMaxIterations ( const int maxIter )
209 {
210 assert(maxIter >= 0);
211 maxIterations_ = maxIter;
212 }
213
214 //Maximum Linear Iterations in total
216 virtual int maxLinearIterations () const
217 {
218 if(maxLinearIterations_ < 0)
219 maxLinearIterations_ = linear().maxIterations();
220 return maxLinearIterations_;
221 }
222
223 virtual void setMaxLinearIterations ( const int maxLinearIter )
224 {
225 assert(maxLinearIter >=0);
226 maxLinearIterations_ = maxLinearIter;
227 }
228
229 virtual int maxLineSearchIterations () const
230 {
231 if(maxLineSearchIterations_ < 0)
232 maxLineSearchIterations_ = parameter_.getValue< int >( keyPrefix_ + "maxlinesearchiterations", std::numeric_limits< int >::max() );
233 return maxLineSearchIterations_;
234 }
235
236 virtual void setMaxLineSearchIterations ( const int maxLineSearchIter )
237 {
238 assert( maxLineSearchIter >= 0);
239 maxLineSearchIterations_ = maxLineSearchIter;
240 }
241
242 // LineSearchMethod: none, simple
243 LIST_OF_INT(LineSearchMethod,
244 none=0,
245 simple=1);
246
247 virtual int lineSearch () const
248 {
249 if( parameter_.exists( keyPrefix_ + "lineSearch" ) )
250 {
251 std::cout << "WARNING: using old parameter name '" << keyPrefix_ + "lineSearch" << "',\n"
252 << "please switch to '" << keyPrefix_ + "linesearch" << "' (all lower caps)!" <<std::endl;
253 return Forcing::to_id( parameter_.getEnum( keyPrefix_ + "lineSearch", LineSearchMethod::names(), LineSearchMethod::none ) );
254 }
255 return Forcing::to_id( parameter_.getEnum( keyPrefix_ + "linesearch", LineSearchMethod::names(), LineSearchMethod::none ) );
256 }
257
258 virtual void setLineSearch ( const int method )
259 {
260 Parameter::append( keyPrefix_ + "linesearch", LineSearchMethod::to_string(method), true );
261 }
262
263 // Forcing: none, eisenstatwalker
264 LIST_OF_INT(Forcing,
265 none = 0, // the provided linear solver tol is used in every iteration
266 eisenstatwalker=1); // Eistenstat-Walker criterion
267
268 virtual int forcing () const
269 {
270 if( parameter_.exists( keyPrefix_ + "linear.tolerance.strategy" ) )
271 {
272 std::string keypref( keyPrefix_ );
273 std::string femsolver("fem.solver.");
274 size_t pos = keypref.find( femsolver );
275 if (pos != std::string::npos)
276 {
277 // If found then erase it from string
278 keypref.erase(pos, femsolver.length());
279 }
280 std::cout << "WARNING: using old parameter name '" << keypref + "linear.tolerance.strategy" << "',\n"
281 << "please switch to '" << keypref + "forcing" << "'!" <<std::endl;
282 return Forcing::to_id( parameter_.getEnum( keyPrefix_ + "linear.tolerance.strategy", Forcing::names(), Forcing::none ) );
283 }
284 return Forcing::to_id( parameter_.getEnum( keyPrefix_ + "forcing", Forcing::names(), Forcing::none ) );
285 }
286
287 virtual void setForcing ( const int strategy )
288 {
289 Parameter::append( keyPrefix_ + "forcing", Forcing::to_string( strategy ), true );
290 }
291
293 virtual bool simplified () const
294 {
295 return parameter_.getValue< bool >( keyPrefix_ + "simplified", 0 );
296 }
297
298 // allow to override the automatic choice of nonlinear or linear solver to
299 // force nonlinear all the time
300 virtual bool forceNonLinear () const
301 {
302 bool v = false;
303 v = parameter_.getValue< bool >(keyPrefix_ + "forcenonlinear", v );
304 return v;
305 }
306
307 private:
308 mutable double tolerance_ = -1;
309 mutable int verbose_ = -1;
310 mutable int maxIterations_ = -1;
311 mutable int maxLinearIterations_ = -1;
312 mutable int maxLineSearchIterations_ = -1;
313 };
314
315
316
317 // NewtonFailure
318 // -------------
319
320 enum class NewtonFailure
321 // : public int
322 {
323 Success = 0,
324 InvalidResidual = 1,
325 IterationsExceeded = 2,
326 LinearIterationsExceeded = 3,
327 LineSearchFailed = 4,
328 TooManyIterations = 5,
329 TooManyLinearIterations = 6,
330 LinearSolverFailed = 7
331 };
332
333
334
335 // NewtonInverseOperator
336 // ---------------------
337
354 template< class JacobianOperator, class LInvOp >
356 : public Operator< typename JacobianOperator::RangeFunctionType, typename JacobianOperator::DomainFunctionType >
357 {
360
362 template <class Obj, bool defaultValue = false >
363 struct SelectPreconditioning
364 {
365 private:
366 template <typename T, typename = bool>
367 struct CheckMember : public std::false_type { };
368
369 template <typename T>
370 struct CheckMember<T, decltype((void) T::preconditioningAvailable, true)> : public std::true_type { };
371
372 template <class T, bool>
373 struct SelectValue
374 {
375 static const bool value = defaultValue;
376 typedef BaseType type;
377 };
378
379 template <class T>
380 struct SelectValue< T, true >
381 {
382 static const bool value = T::preconditioningAvailable;
383 typedef typename T::PreconditionerType type;
384 };
385 public:
386 static constexpr bool value = SelectValue< Obj, CheckMember< Obj >::value >::value;
387 typedef typename SelectValue< Obj, CheckMember< Obj >::value >::type type;
388 };
389
390 public:
392 typedef JacobianOperator JacobianOperatorType;
393
396
399
401 static constexpr bool preconditioningAvailable = SelectPreconditioning< LinearInverseOperatorType > :: value;
402 typedef typename SelectPreconditioning< LinearInverseOperatorType > :: type PreconditionerType;
403
404 typedef typename BaseType::DomainFunctionType DomainFunctionType;
405 typedef typename BaseType::RangeFunctionType RangeFunctionType;
406
407 typedef typename BaseType::DomainFieldType DomainFieldType;
408
409 typedef NewtonParameter<typename LinearInverseOperatorType::SolverParameterType> ParameterType;
410
411 typedef std::function< bool ( const RangeFunctionType &w, const RangeFunctionType &dw, double residualNorm ) > ErrorMeasureType;
412
414 typedef Impl::SolverInfo SolverInfoType;
415
433 // main constructor
434 NewtonInverseOperator ( LinearInverseOperatorType jInv, const DomainFieldType &epsilon, const ParameterType &parameter )
435 : verbose_( parameter.verbose() ),
436 maxLineSearchIterations_( parameter.maxLineSearchIterations() ),
437 jInv_( std::move( jInv ) ),
438 parameter_(parameter),
439 lsMethod_( parameter.lineSearch() ),
440 finished_( [ epsilon ] ( const RangeFunctionType &w, const RangeFunctionType &dw, double res ) { return res < epsilon; } ),
441 forcing_ ( parameter.forcing() ),
442 eisenstatWalker_ ( epsilon ),
443 timing_(3, 0.0)
444 {
445 }
446
447
453 explicit NewtonInverseOperator ( const ParameterType &parameter = ParameterType( Parameter::container() ) )
454 : NewtonInverseOperator( parameter.tolerance(), parameter )
455 {
456 // std::cout << "in Newton inv op should use:" << parameter.linear().solverMethod({SolverParameter::gmres,SolverParameter::bicgstab,SolverParameter::minres}) << std::endl;
457 }
458
464 NewtonInverseOperator ( const DomainFieldType &epsilon, const ParameterType &parameter )
466 LinearInverseOperatorType( parameter.linear() ),
467 epsilon, parameter )
468 {}
469
475 NewtonInverseOperator ( const DomainFieldType &epsilon,
476 const ParameterReader &parameter = Parameter::container() )
477 : NewtonInverseOperator( epsilon, ParameterType( parameter ) )
478 {}
479
480
488 void setErrorMeasure ( ErrorMeasureType finished ) { finished_ = std::move( finished ); }
489
490 EisenstatWalkerStrategy& eisenstatWalker () { return eisenstatWalker_; }
491
492 void bind ( const OperatorType &op ) { op_ = &op; }
493
494 void bind ( const OperatorType &op, const PreconditionerType& preconditioner )
495 {
496 bind( op );
498 preconditioner_ = &preconditioner;
499 else
500 std::cerr << "WARNING: linear inverse operator does not support external preconditioning!" << std::endl;
501 }
502
503 void unbind () { op_ = nullptr; preconditioner_ = nullptr; }
504
505 virtual void operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const;
506
507 int iterations () const { return iterations_; }
508 void setMaxIterations ( int maxIterations ) { parameter_.setMaxIterations( maxIterations ); }
509 int linearIterations () const { return linearIterations_; }
510 void setMaxLinearIterations ( int maxLinearIterations ) { parameter_.setMaxLinearIterations( maxLinearIterations ); }
511 void updateLinearTolerance () const {
512 if (forcing_ == ParameterType::Forcing::eisenstatwalker) {
513 double newTol = eisenstatWalker_.nextLinearTolerance( delta_ );
514 jInv_.parameter().setTolerance( newTol );
515 }
516 }
517 bool verbose() const { return verbose_ && Parameter::verbose( Parameter::solverStatistics ); }
518 double residual () const { return delta_; }
519
520 NewtonFailure failed () const
521 {
522 // check for finite |residual| - this also works for -ffinite-math-only (gcc)
523 // nan test not working with optimization flags...
524 if( !(delta_ < std::numeric_limits< DomainFieldType >::max()) || std::isnan( delta_ ) )
525 return NewtonFailure::InvalidResidual;
526 else if( iterations_ >= parameter_.maxIterations() )
527 return NewtonFailure::TooManyIterations;
528 else if( linearIterations_ >= parameter_.maxLinearIterations() )
529 return NewtonFailure::TooManyLinearIterations;
530 else if( linearIterations_ < 0)
531 return NewtonFailure::LinearSolverFailed;
532 else if( !stepCompleted_ )
533 return NewtonFailure::LineSearchFailed;
534 else
535 return NewtonFailure::Success;
536 }
537
538 bool converged () const { return failed() == NewtonFailure::Success; }
539
540 virtual int lineSearch(RangeFunctionType &w, RangeFunctionType &dw,
541 const DomainFunctionType &u, DomainFunctionType &residual) const
542 {
543 double deltaOld = delta_;
544 delta_ = std::sqrt( residual.scalarProductDofs( residual ) );
546 return 0;
547 if (failed() == NewtonFailure::InvalidResidual)
548 {
549 double test = dw.scalarProductDofs( dw );
551 !std::isnan(test)) )
552 delta_ = 2.*deltaOld; // enter line search
553 }
554 double factor = 1.0;
555 int noLineSearch = (delta_ < deltaOld)?1:0;
556 int lineSearchIteration = 0;
557 const bool lsVerbose = verbose() && Parameter::verbose( Parameter::extendedStatistics );
558 while (delta_ >= deltaOld)
559 {
560 double deltaPrev = delta_;
561 factor *= 0.5;
562 if( lsVerbose )
563 std::cout << " line search:" << delta_ << ">" << deltaOld << std::endl;
564 if (std::abs(delta_-deltaOld) < 1e-5*delta_) // || !converged()) // line search not working
565 return -1; // failed
566 w.axpy(factor,dw);
567 (*op_)( w, residual );
568 residual -= u;
569 delta_ = std::sqrt( residual.scalarProductDofs( residual ) );
570 if (std::abs(delta_-deltaPrev) < 1e-15)
571 return -1;
572 if (failed() == NewtonFailure::InvalidResidual)
573 delta_ = 2.*deltaOld; // remain in line search
574
575 ++lineSearchIteration;
576 if( lineSearchIteration >= maxLineSearchIterations_ )
577 return -1; // failed
578 }
579 return noLineSearch;
580 }
581
583 const std::vector<double>& timing() const { return timing_; }
584
586 SolverInfoType info() const { return SolverInfoType( converged(), linearIterations(), iterations(), timing() ); }
587
588 protected:
589 void bindOperatorAndPreconditioner( JacobianOperatorType& jOp ) const
590 {
591 // if preconditioner was given pass it on to linear solver
592 if constexpr ( preconditioningAvailable )
593 {
594 if( preconditioner_ )
595 {
596 jInv_.bind( jOp, *preconditioner_ );
597 return ;
598 }
599 }
600
601 // if preconditioning not enabled or set, then only set jOp
602 jInv_.bind( jOp );
603 }
604
605 // hold pointer to jacobian operator, if memory reallocation is needed, the operator should know how to handle this.
606 template< class ... Args>
607 JacobianOperatorType& jacobian ( Args && ... args ) const
608 {
609 if( !jOp_ )
610 jOp_.reset( new JacobianOperatorType( std::forward< Args >( args ) ...) ); //, parameter_.parameter() ) );
611 return *jOp_;
612 }
613
614 protected:
615 const OperatorType *op_ = nullptr;
616 const PreconditionerType* preconditioner_ = nullptr;
617
618 const bool verbose_;
619 const int maxLineSearchIterations_;
620
621 mutable DomainFieldType delta_;
622 mutable int iterations_;
623 mutable int linearIterations_;
624 mutable LinearInverseOperatorType jInv_;
625 mutable std::unique_ptr< JacobianOperatorType > jOp_;
626 ParameterType parameter_;
627 mutable int stepCompleted_;
628 const int lsMethod_;
629 ErrorMeasureType finished_;
630 const int forcing_;
631 EisenstatWalkerStrategy eisenstatWalker_;
632
633 mutable std::vector<double> timing_;
634 };
635
636
637 // Implementation of NewtonInverseOperator
638 // ---------------------------------------
639
640 template< class JacobianOperator, class LInvOp >
641 inline void NewtonInverseOperator< JacobianOperator, LInvOp >
642 ::operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const
643 {
644 assert( op_ );
645 std::fill(timing_.begin(), timing_.end(), 0.0 );
646
647 // obtain information about operator to invert
648 const bool nonlinear = op_->nonlinear() || parameter_.forceNonLinear();
649
650 Dune::Timer allTimer;
651 DomainFunctionType residual( u );
652 RangeFunctionType dw( w );
653
654 Dune::Timer jacTimer;
655 double jacTime = 0.0;
656 JacobianOperatorType& jOp = jacobian( "jacobianOperator", dw.space(), u.space(), parameter_.solverParameter() );
657 jacTime += jacTimer.elapsed();
658
659 stepCompleted_ = true;
660 iterations_ = 0;
661 linearIterations_ = 0;
662 // compute initial residual
663 (*op_)( w, residual ); // r=S[w], r=w-g on bnd
664 residual -= u; // r=S[w]-u, r=w-g-u on bnd (note: we should start with w=g+u on bnd so r=0)
665 delta_ = std::sqrt( residual.scalarProductDofs( residual ) );
666 updateLinearTolerance();
667
668 // this is true for Newton, and false simplified Newton after first iteration
669 bool computeJacobian = true;
670 const bool simplifiedNewton = parameter_.simplified();
671
672 const bool newtonVerbose = verbose() && nonlinear;
673 if( newtonVerbose )
674 {
675 std::cout << "Start Newton: tol = " << parameter_.tolerance() << " (forcing = " << ParameterType::Forcing::to_string(forcing_) << " | linear tol = " << parameter_.linear().tolerance() << ")"<<std::endl;
676 std::cout << "Newton iteration " << iterations_ << ": |residual| = " << delta_;
677 }
678 while( true )
679 {
680 if( newtonVerbose )
681 std::cout << std::endl;
682
683 if( computeJacobian )
684 {
685 // evaluate operator's Jacobian
686 jacTimer.reset();
687 (*op_).jacobian( w, jOp );
688
689 // if simplified Newton is activated do not compute Jacobian again
690 computeJacobian = ! simplifiedNewton;
691 }
692
693 // bind jOp to jInv including preconditioner if enabled and set
694 bindOperatorAndPreconditioner( jOp );
695 jacTime += jacTimer.elapsed();
696
697 if ( parameter_.maxLinearIterations() - linearIterations_ <= 0 )
698 break;
699 jInv_.setMaxIterations( parameter_.maxLinearIterations() - linearIterations_ );
700
701 dw.clear();
702 jInv_( residual, dw ); // dw = DS[w]^{-1}(S[w]-u)
703 // dw = w-g-u on bnd
704 if (jInv_.iterations() < 0) // iterations are negative if solver didn't converge
705 {
706 linearIterations_ = jInv_.iterations();
707 break;
708 }
709 linearIterations_ += jInv_.iterations();
710 w -= dw; // w = w - DS[w]^{-1}(S[w]-u)
711 // w = g+u
712
713 // the following only for nonlinear problems
714 if( nonlinear )
715 {
716 // compute new residual
717 (*op_)( w, residual ); // res = S[w]
718 // res = w-g = g+u-g = u
719 residual -= u; // res = S[w] - u
720 // res = 0
721
722 // line search if enabled
723 int ls = lineSearch(w,dw,u,residual);
724 stepCompleted_ = ls >= 0;
725 updateLinearTolerance();
726 ++iterations_;
727
728 if( newtonVerbose )
729 std::cout << "Newton iteration " << iterations_ << ": |residual| = " << delta_ << std::flush;
730 // if ( (ls==1 && finished_(w, dw, delta_)) || !converged())
731 if ( (finished_(w, dw, delta_)) || !converged())
732 {
733 if( newtonVerbose )
734 {
735 std::cout << std::endl;
736 std::cout << "Linear iterations: " << linearIterations_ << std::endl;
737 }
738 break;
739 }
740 }
741 else // in linear case do not continue
742 break ;
743 } // end Newton loop
744
745 if( newtonVerbose )
746 std::cout << std::endl;
747
748 jInv_.unbind();
749
750 // store time measurements
751 timing_[0] = allTimer.elapsed();
752 timing_[1] = jacTime;
753 timing_[2] = timing_[0] - jacTime;
754 }
755
756 } // namespace Fem
757
758} // namespace Dune
759
760#endif // #ifndef DUNE_FEM_NEWTONINVERSEOPERATOR_HH
abstract differentiable operator
Definition: differentiableoperator.hh:29
Adaptive tolerance selection for linear solver.
Definition: newtoninverseoperator.hh:38
EisenstatWalkerStrategy(const double newtonTolerance)
Definition: newtoninverseoperator.hh:50
inverse operator based on a newton scheme
Definition: newtoninverseoperator.hh:357
NewtonInverseOperator(LinearInverseOperatorType jInv, const DomainFieldType &epsilon, const ParameterType &parameter)
Definition: newtoninverseoperator.hh:434
Impl::SolverInfo SolverInfoType
performance info about last solver call
Definition: newtoninverseoperator.hh:414
static constexpr bool preconditioningAvailable
type of preconditioner for linear solver
Definition: newtoninverseoperator.hh:401
JacobianOperator JacobianOperatorType
type of operator's Jacobian
Definition: newtoninverseoperator.hh:392
NewtonInverseOperator(const ParameterType &parameter=ParameterType(Parameter::container()))
Definition: newtoninverseoperator.hh:453
const std::vector< double > & timing() const
returns [overall, jacobian, solve] timings in seconds for last operator () call.
Definition: newtoninverseoperator.hh:583
SolverInfoType info() const
return performance information about last solver run *‍/
Definition: newtoninverseoperator.hh:586
NewtonInverseOperator(const DomainFieldType &epsilon, const ParameterReader &parameter=Parameter::container())
Definition: newtoninverseoperator.hh:475
void setErrorMeasure(ErrorMeasureType finished)
Definition: newtoninverseoperator.hh:488
DifferentiableOperator< JacobianOperatorType > OperatorType
type of operator to invert
Definition: newtoninverseoperator.hh:395
LInvOp LinearInverseOperatorType
type of linear inverse operator
Definition: newtoninverseoperator.hh:398
NewtonInverseOperator(const DomainFieldType &epsilon, const ParameterType &parameter)
Definition: newtoninverseoperator.hh:464
static bool verbose()
obtain the cached value for fem.verbose with default verbosity level 2
Definition: parameter.hh:466
static void append(int &argc, char **argv)
add parameters from the command line RangeType gRight;
Definition: parameter.hh:219
A simple stop watch.
Definition: timer.hh:43
void reset() noexcept
Reset timer while keeping the running/stopped state.
Definition: timer.hh:57
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:77
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType none(unsigned int dim)
Returns a GeometryType representing a singular of dimension dim.
Definition: type.hh:471
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Dune namespace.
Definition: alignedallocator.hh:13
STL namespace.
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
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)