1#ifndef DUNE_FEM_NEWTONINVERSEOPERATOR_HH
2#define DUNE_FEM_NEWTONINVERSEOPERATOR_HH
27 template <
class SolverParam = SolverParameter>
40 NewtonParameter(
const SolverParam& baseParameter,
const std::string keyPrefix =
"fem.solver.newton." )
46 template <class Parameter, std::enable_if_t<!std::is_base_of<SolverParam,Parameter>::value && !std::is_same<Parameter,ParameterReader>::value,
int> i=0>
53 template <class ParamReader, std::enable_if_t<!std::is_same<ParamReader,SolverParam>::value && std::is_same<ParamReader,ParameterReader>::value,
int> i=0>
70 maxLinearIterations_ = -1;
71 maxLineSearchIterations_ = -1;
104 verbose_ = verb ? 1 : 0;
109 if(maxIterations_ < 0)
111 return maxIterations_;
116 assert(maxIter >= 0);
117 maxIterations_ = maxIter;
124 if(maxLinearIterations_ < 0)
125 maxLinearIterations_ =
linear().maxIterations();
126 return maxLinearIterations_;
131 assert(maxLinearIter >=0);
132 maxLinearIterations_ = maxLinearIter;
137 if(maxLineSearchIterations_ < 0)
139 return maxLineSearchIterations_;
144 assert( maxLineSearchIter >= 0);
145 maxLineSearchIterations_ = maxLineSearchIter;
155 const std::string lineSearchMethods[] = {
"none",
"simple" };
161 const std::string lineSearchMethods[] = {
"none",
"simple" };
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;
206 template<
class JacobianOperator,
class LInvOp >
208 :
public Operator< typename JacobianOperator::RangeFunctionType, typename JacobianOperator::DomainFunctionType >
252 maxLineSearchIterations_( parameter.maxLineSearchIterations() ),
253 jInv_(
std::move( jInv ) ),
254 parameter_(parameter),
319 if( !(delta_ < std::numeric_limits< DomainFieldType >::max()) || std::isnan( delta_ ) )
325 else if( linearIterations_ < 0)
327 else if( !stepCompleted_ )
338 double deltaOld = delta_;
339 delta_ =
std::sqrt( residual.scalarProductDofs( residual ) );
340 if (lsMethod_ == ParameterType::LineSearchMethod::none)
344 double test = dw.scalarProductDofs( dw );
345 if (! (test < std::numeric_limits< DomainFieldType >::max() &&
347 delta_ = 2.*deltaOld;
350 int noLineSearch = (delta_ < deltaOld)?1:0;
351 int lineSearchIteration = 0;
352 while (delta_ >= deltaOld)
354 double deltaPrev = delta_;
357 std::cerr <<
" line search:" << delta_ <<
">" << deltaOld << std::endl;
358 if (
std::abs(delta_-deltaOld) < 1e-5*delta_)
361 (*op_)( w, residual );
363 delta_ =
std::sqrt( residual.scalarProductDofs( residual ) );
364 if (
std::abs(delta_-deltaPrev) < 1e-15)
367 delta_ = 2.*deltaOld;
369 ++lineSearchIteration;
370 if( lineSearchIteration >= maxLineSearchIterations_ )
378 template<
class ... Args>
390 const int maxLineSearchIterations_;
393 mutable int iterations_;
394 mutable int linearIterations_;
396 mutable std::unique_ptr< JacobianOperatorType > jOp_;
398 mutable int stepCompleted_;
407 template<
class JacobianOperator,
class LInvOp >
408 inline void NewtonInverseOperator< JacobianOperator, LInvOp >
417 stepCompleted_ =
true;
419 linearIterations_ = 0;
421 (*op_)( w, residual );
423 delta_ =
std::sqrt( residual.scalarProductDofs( residual ) );
426 std::cerr <<
"Newton iteration " << iterations_ <<
": |residual| = " << delta_;
430 std::cerr << std::endl;
432 (*op_).jacobian( w, jOp );
438 if ( parameter_.maxLinearIterations() - linearIterations_ <= 0 )
440 jInv_.setMaxIterations( parameter_.maxLinearIterations() - linearIterations_ );
443 jInv_( residual, dw );
444 if (jInv_.iterations() < 0)
446 linearIterations_ = jInv_.iterations();
449 linearIterations_ += jInv_.iterations();
452 (*op_)( w, residual );
454 int ls = lineSearch(w,dw,u,residual);
455 stepCompleted_ = ls >= 0;
458 std::cerr <<
"Newton iteration " << iterations_ <<
": |residual| = " << delta_ << std::flush;
460 if ( (finished_(w, dw, delta_)) || !converged())
463 std::cerr << std::endl;
468 std::cerr << std::endl;
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
@ LinearIterationsExceeded
@ TooManyLinearIterations
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 ¶meter, 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 ¶meter)
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 ¶meter=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 ¶meter=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 ¶meter)
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