3 #ifndef DUNE_PDELAB_NEWTON_NEWTON_HH 4 #define DUNE_PDELAB_NEWTON_NEWTON_HH 13 #include <dune/common/exceptions.hh> 14 #include <dune/common/ios_state.hh> 15 #include <dune/common/timer.hh> 16 #include <dune/common/parametertree.hh> 32 template<
class RFType>
42 first_defect(0.0), defect(0.0), assembler_time(0.0), linear_solver_time(0.0),
43 linear_solver_iterations(0) {}
46 template<
class GOS,
class TrlV,
class TstV>
49 typedef GOS GridOperator;
50 typedef TrlV TrialVector;
51 typedef TstV TestVector;
53 typedef typename TestVector::ElementType RFType;
54 typedef typename GOS::Traits::Jacobian
Matrix;
63 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
66 verbosity_level_ = verbosity_level;
91 std::shared_ptr<TrialVector>
z_;
92 std::shared_ptr<TestVector>
r_;
93 std::shared_ptr<Matrix>
A_;
106 , verbosity_level_(1)
109 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
110 verbosity_level_ = 0;
116 , verbosity_level_(1)
119 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
120 verbosity_level_ = 0;
125 virtual bool terminate() = 0;
126 virtual void prepare_step(Matrix& A, TestVector& r) = 0;
127 virtual void line_search(TrialVector& z, TestVector& r) = 0;
128 virtual void defect(TestVector& r) = 0;
131 template<
class GOS,
class S,
class TrlV,
class TstV>
135 typedef GOS GridOperator;
136 typedef TrlV TrialVector;
137 typedef TstV TestVector;
139 typedef typename TestVector::ElementType RFType;
140 typedef typename GOS::Traits::Jacobian
Matrix;
148 , result_valid_(false)
154 , result_valid_(false)
159 void apply(TrialVector& u_);
165 "NewtonSolver::result() called before NewtonSolver::solve()");
173 this->gridoperator_.residual(*this->u_, r);
174 this->res_.defect = this->solver_.norm(r);
175 if (!std::isfinite(this->res_.defect))
177 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
182 void linearSolve(Matrix& A, TrialVector& z, TestVector& r)
const 184 if (this->verbosity_level_ >= 4)
185 std::cout <<
" Solving linear system..." << std::endl;
187 this->solver_.apply(A, z, r, this->linear_reduction_);
189 ios_base_all_saver restorer(std::cout);
191 if (!this->solver_.result().converged)
193 "NewtonSolver::linearSolve(): Linear solver did not converge " 194 "in " << this->solver_.result().iterations <<
" iterations");
195 if (this->verbosity_level_ >= 4)
196 std::cout <<
" linear solver iterations: " 197 << std::setw(12) << solver_.result().iterations << std::endl
198 <<
" linear defect reduction: " 199 << std::setw(12) << std::setprecision(4) << std::scientific
200 << solver_.result().reduction << std::endl;
207 template<
class GOS,
class S,
class TrlV,
class TstV>
214 template<
class GOS,
class S,
class TrlV,
class TstV>
217 this->res_.iterations = 0;
218 this->res_.converged =
false;
219 this->res_.reduction = 1.0;
220 this->res_.conv_rate = 1.0;
221 this->res_.elapsed = 0.0;
222 this->res_.assembler_time = 0.0;
223 this->res_.linear_solver_time = 0.0;
224 this->res_.linear_solver_iterations = 0;
225 result_valid_ =
true;
232 this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
239 this->defect(*this->r_);
240 this->res_.first_defect = this->res_.defect;
241 this->prev_defect_ = this->res_.defect;
243 if (this->verbosity_level_ >= 2)
246 ios_base_all_saver restorer(std::cout);
247 std::cout <<
" Initial defect: " 248 << std::setw(12) << std::setprecision(4) << std::scientific
249 << this->res_.defect << std::endl;
254 this->A_ = std::make_shared<Matrix>(this->gridoperator_);
258 this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
261 while (!this->terminate())
263 if (this->verbosity_level_ >= 3)
264 std::cout <<
" Newton iteration " << this->res_.iterations
265 <<
" --------------------------------" << std::endl;
267 Timer assembler_timer;
276 this->prepare_step(*this->A_,*this->r_);
280 this->res_.assembler_time += assembler_timer.elapsed();
283 double assembler_time = assembler_timer.elapsed();
284 this->res_.assembler_time += assembler_time;
285 if (this->verbosity_level_ >= 3)
286 std::cout <<
" matrix assembly time: " 287 << std::setw(12) << std::setprecision(4) << std::scientific
288 << assembler_time << std::endl;
290 Timer linear_solver_timer;
297 this->linearSolve(*this->A_, *this->z_, *this->r_);
301 this->res_.linear_solver_time += linear_solver_timer.elapsed();
302 this->res_.linear_solver_iterations += this->solver_.result().iterations;
305 double linear_solver_time = linear_solver_timer.elapsed();
306 this->res_.linear_solver_time += linear_solver_time;
307 this->res_.linear_solver_iterations += this->solver_.result().iterations;
313 this->line_search(*this->z_, *this->r_);
317 if (this->reassembled_)
319 if (this->verbosity_level_ >= 3)
320 std::cout <<
" line search failed - trying again with reassembled matrix" << std::endl;
324 this->res_.reduction = this->res_.defect/this->res_.first_defect;
325 this->res_.iterations++;
326 this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
329 ios_base_all_saver restorer(std::cout);
331 if (this->verbosity_level_ >= 3)
332 std::cout <<
" linear solver time: " 333 << std::setw(12) << std::setprecision(4) << std::scientific
334 << linear_solver_time << std::endl
335 <<
" defect reduction (this iteration):" 336 << std::setw(12) << std::setprecision(4) << std::scientific
337 << this->res_.defect/this->prev_defect_ << std::endl
338 <<
" defect reduction (total): " 339 << std::setw(12) << std::setprecision(4) << std::scientific
340 << this->res_.reduction << std::endl
342 << std::setw(12) << std::setprecision(4) << std::scientific
343 << this->res_.defect << std::endl;
344 if (this->verbosity_level_ == 2)
345 std::cout <<
" Newton iteration " << std::setw(2) << this->res_.iterations
347 << std::setw(12) << std::setprecision(4) << std::scientific
349 <<
". Reduction (this): " 350 << std::setw(12) << std::setprecision(4) << std::scientific
351 << this->res_.defect/this->prev_defect_
352 <<
". Reduction (total): " 353 << std::setw(12) << std::setprecision(4) << std::scientific
354 << this->res_.reduction << std::endl;
359 this->res_.elapsed = timer.elapsed();
362 this->res_.elapsed = timer.elapsed();
364 ios_base_all_saver restorer(std::cout);
366 if (this->verbosity_level_ == 1)
367 std::cout <<
" Newton converged after " << std::setw(2) << this->res_.iterations
368 <<
" iterations. Reduction: " 369 << std::setw(12) << std::setprecision(4) << std::scientific
370 << this->res_.reduction
371 <<
" (" << std::setprecision(4) << this->res_.elapsed <<
"s)" 374 if(!this->keep_matrix_)
378 template<
class GOS,
class TrlV,
class TstV>
381 typedef GOS GridOperator;
382 typedef TrlV TrialVector;
384 typedef typename TstV::ElementType RFType;
390 , force_iteration_(false)
392 this->reduction_ = 1
e-8;
393 this->abs_limit_ = 1
e-12;
399 , force_iteration_(false)
401 this->reduction_ = 1
e-8;
402 this->abs_limit_ = 1
e-12;
407 this->reduction_ = reduction;
417 force_iteration_ = force_iteration;
422 this->abs_limit_ = abs_limit_;
427 if (force_iteration_ && this->res_.iterations == 0)
429 this->res_.converged = this->res_.defect < this->abs_limit_
430 || this->res_.defect < this->res_.first_defect * this->reduction_;
431 if (this->res_.iterations >= maxit_ && !this->res_.converged)
433 "NewtonTerminate::terminate(): Maximum iteration count reached");
434 return this->res_.converged;
439 bool force_iteration_;
442 template<
class GOS,
class TrlV,
class TstV>
445 typedef GOS GridOperator;
446 typedef TrlV TrialVector;
448 typedef typename TstV::ElementType RFType;
449 typedef typename GOS::Traits::Jacobian
Matrix;
454 , min_linear_reduction_(1
e-3)
455 , fixed_linear_reduction_(0.0)
456 , reassemble_threshold_(0.0)
461 , min_linear_reduction_(1
e-3)
462 , fixed_linear_reduction_(0.0)
463 , reassemble_threshold_(0.0)
474 min_linear_reduction_ = min_linear_reduction;
483 fixed_linear_reduction_ = fixed_linear_reduction;
495 reassemble_threshold_ = reassemble_threshold;
500 this->reassembled_ =
false;
501 if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
503 if (this->verbosity_level_ >= 3)
504 std::cout <<
" Reassembling matrix..." << std::endl;
506 this->gridoperator_.jacobian(*this->u_, A);
507 this->reassembled_ =
true;
510 if (fixed_linear_reduction_ ==
true)
511 this->linear_reduction_ = min_linear_reduction_;
515 std::max(this->res_.first_defect * this->reduction_,
526 if ( stop_defect/(10*this->res_.defect) >
527 this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_) )
528 this->linear_reduction_ =
529 stop_defect/(10*this->res_.defect);
531 this->linear_reduction_ =
532 std::min(min_linear_reduction_,this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_));
535 this->prev_defect_ = this->res_.defect;
537 ios_base_all_saver restorer(std::cout);
539 if (this->verbosity_level_ >= 3)
540 std::cout <<
" requested linear reduction: " 541 << std::setw(12) << std::setprecision(4) << std::scientific
542 << this->linear_reduction_ << std::endl;
546 RFType min_linear_reduction_;
547 bool fixed_linear_reduction_;
548 RFType reassemble_threshold_;
551 template<
class GOS,
class TrlV,
class TstV>
554 typedef GOS GridOperator;
555 typedef TrlV TrialVector;
556 typedef TstV TestVector;
558 typedef typename TestVector::ElementType RFType;
569 hackbuschReuskenAcceptBest };
573 , strategy_(hackbuschReusken)
575 , damping_factor_(0.5)
580 , strategy_(hackbuschReusken)
582 , damping_factor_(0.5)
587 strategy_ = strategy;
592 strategy_ = strategyFromName(strategy);
602 damping_factor_ = damping_factor;
607 if (strategy_ == noLineSearch)
609 this->u_->axpy(-1.0, z);
614 if (this->verbosity_level_ >= 4)
615 std::cout <<
" Performing line search..." << std::endl;
617 RFType best_lambda = 0.0;
618 RFType best_defect = this->res_.defect;
619 TrialVector prev_u(*this->u_);
621 ios_base_all_saver restorer(std::cout);
625 if (this->verbosity_level_ >= 4)
626 std::cout <<
" trying line search damping factor: " 627 << std::setw(12) << std::setprecision(4) << std::scientific
631 this->u_->axpy(-lambda, z);
637 if (this->verbosity_level_ >= 4)
638 std::cout <<
" NaNs detected" << std::endl;
641 if (this->res_.defect <= (1.0 - lambda/4) * this->prev_defect_)
643 if (this->verbosity_level_ >= 4)
644 std::cout <<
" line search converged" << std::endl;
648 if (this->res_.defect < best_defect)
650 best_defect = this->res_.defect;
651 best_lambda = lambda;
656 if (this->verbosity_level_ >= 4)
657 std::cout <<
" max line search iterations exceeded" << std::endl;
660 case hackbuschReusken:
664 "NewtonLineSearch::line_search(): line search failed, " 665 "max iteration count reached, " 666 "defect did not improve enough");
667 case hackbuschReuskenAcceptBest:
668 if (best_lambda == 0.0)
673 "NewtonLineSearch::line_search(): line search failed, " 674 "max iteration count reached, " 675 "defect did not improve in any of the iterations");
677 if (best_lambda != lambda)
680 this->u_->axpy(-best_lambda, z);
690 lambda *= damping_factor_;
693 if (this->verbosity_level_ >= 4)
694 std::cout <<
" line search damping factor: " 695 << std::setw(12) << std::setprecision(4) << std::scientific
696 << lambda << std::endl;
702 if (s ==
"noLineSearch")
704 if (s ==
"hackbuschReusken")
705 return hackbuschReusken;
706 if (s ==
"hackbuschReuskenAcceptBest")
707 return hackbuschReuskenAcceptBest;
708 DUNE_THROW(
Exception,
"unknown line search strategy" << s);
714 RFType damping_factor_;
717 template<
class GOS,
class S,
class TrlV,
class TstV = TrlV>
723 typedef GOS GridOperator;
725 typedef TrlV TrialVector;
728 Newton(
const GridOperator& go, TrialVector& u_, Solver& solver_)
735 Newton(
const GridOperator& go, Solver& solver_)
766 typedef typename TstV::ElementType RFType;
767 if (param.hasKey(
"VerbosityLevel"))
768 this->setVerbosityLevel(
769 param.get<
unsigned int>(
"VerbosityLevel"));
770 if (param.hasKey(
"Reduction"))
772 param.get<RFType>(
"Reduction"));
773 if (param.hasKey(
"MaxIterations"))
774 this->setMaxIterations(
775 param.get<
unsigned int>(
"MaxIterations"));
776 if (param.hasKey(
"ForceIteration"))
777 this->setForceIteration(
778 param.get<
bool>(
"ForceIteration"));
779 if (param.hasKey(
"AbsoluteLimit"))
780 this->setAbsoluteLimit(
781 param.get<RFType>(
"AbsoluteLimit"));
782 if (param.hasKey(
"MinLinearReduction"))
783 this->setMinLinearReduction(
784 param.get<RFType>(
"MinLinearReduction"));
785 if (param.hasKey(
"FixedLinearReduction"))
786 this->setFixedLinearReduction(
787 param.get<
bool>(
"FixedLinearReduction"));
788 if (param.hasKey(
"ReassembleThreshold"))
789 this->setReassembleThreshold(
790 param.get<RFType>(
"ReassembleThreshold"));
791 if (param.hasKey(
"LineSearchStrategy"))
792 this->setLineSearchStrategy(
793 param.get<std::string>(
"LineSearchStrategy"));
794 if (param.hasKey(
"LineSearchMaxIterations"))
795 this->setLineSearchMaxIterations(
796 param.get<
unsigned int>(
"LineSearchMaxIterations"));
797 if (param.hasKey(
"LineSearchDampingFactor"))
798 this->setLineSearchDampingFactor(
799 param.get<RFType>(
"LineSearchDampingFactor"));
800 if (param.hasKey(
"KeepMatrix"))
802 param.get<
bool>(
"KeepMatrix"));
808 #endif // DUNE_PDELAB_NEWTON_NEWTON_HH perform a linear search for the optimal damping parameter with multiples of damping ...
Definition: newton.hh:567
Definition: newton.hh:379
RFType abs_limit_
Definition: newton.hh:100
virtual void defect(TestVector &r)
Definition: newton.hh:170
NewtonSolver(const GridOperator &go, TrialVector &u_, Solver &solver)
Definition: newton.hh:145
void setLineSearchMaxIterations(unsigned int maxit)
Definition: newton.hh:595
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
const GridOperator & gridoperator_
Definition: newton.hh:89
const std::string s
Definition: function.hh:1101
NewtonTerminate(const GridOperator &go)
Definition: newton.hh:396
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: newton.hh:82
std::shared_ptr< Matrix > A_
Definition: newton.hh:93
RFType defect
Definition: newton.hh:36
bool keep_matrix_
Definition: newton.hh:101
std::shared_ptr< TestVector > r_
Definition: newton.hh:92
NewtonPrepareStep(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:452
RFType linear_reduction_
Definition: newton.hh:97
Result res_
Definition: newton.hh:94
void apply()
Definition: newton.hh:215
void setReduction(RFType reduction)
Definition: newton.hh:405
Definition: newton.hh:552
Definition: newton.hh:132
void setVerbosityLevel(unsigned int verbosity_level)
Definition: newton.hh:61
NewtonLineSearch(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:571
void setLineSearchStrategy(std::string strategy)
Definition: newton.hh:590
RFType reduction_
Definition: newton.hh:99
typename impl::BackendMatrixSelector< Backend, VU, VV, E >::Type Matrix
alias of the return type of BackendMatrixSelector
Definition: backend/interface.hh:127
void setMinLinearReduction(RFType min_linear_reduction)
set the minimal reduction in the linear solver
Definition: newton.hh:472
don't do any line search or damping
Definition: newton.hh:563
NewtonTerminate(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:387
Definition: newton.hh:443
std::shared_ptr< TrialVector > z_
Definition: newton.hh:91
virtual bool terminate()
Definition: newton.hh:425
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: newton.hh:70
void setLineSearchStrategy(Strategy strategy)
Definition: newton.hh:585
unsigned int verbosity_level_
Definition: newton.hh:95
virtual void line_search(TrialVector &z, TestVector &r)
Definition: newton.hh:605
virtual ~NewtonBase()
Definition: newton.hh:123
int linear_solver_iterations
Definition: newton.hh:39
void setMaxIterations(unsigned int maxit)
Definition: newton.hh:410
NewtonLineSearch(const GridOperator &go)
Definition: newton.hh:578
NewtonResult< RFType > Result
Definition: newton.hh:143
void setParameters(const Dune::ParameterTree ¶m)
interpret a parameter tree as a set of options for the newton solver
Definition: newton.hh:764
const Result & result() const
Definition: newton.hh:161
void setLineSearchDampingFactor(RFType damping_factor)
Definition: newton.hh:600
double assembler_time
Definition: newton.hh:37
NewtonResult< RFType > Result
Definition: newton.hh:59
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
RFType prev_defect_
Definition: newton.hh:96
const Entity & e
Definition: localfunctionspace.hh:111
void setReassembleThreshold(RFType reassemble_threshold)
set a threshold, when the linear operator is reassembled
Definition: newton.hh:493
TrialVector * u_
Definition: newton.hh:90
NewtonBase(const GridOperator &go)
Definition: newton.hh:113
virtual void prepare_step(Matrix &A, TstV &)
Definition: newton.hh:498
NewtonPrepareStep(const GridOperator &go)
Definition: newton.hh:459
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: newton.hh:76
Newton(const GridOperator &go, Solver &solver_)
Definition: newton.hh:735
bool reassembled_
Definition: newton.hh:98
double linear_solver_time
Definition: newton.hh:38
void setAbsoluteLimit(RFType abs_limit_)
Definition: newton.hh:420
NewtonBase(const GridOperator &go, TrialVector &u)
Definition: newton.hh:103
NewtonResult()
Definition: newton.hh:41
NewtonSolver(const GridOperator &go, Solver &solver)
Definition: newton.hh:151
Strategy strategyFromName(const std::string &s)
Definition: newton.hh:701
void setFixedLinearReduction(bool fixed_linear_reduction)
set a fixed reduction in the linear solver (overwrites setMinLinearReduction)
Definition: newton.hh:481
void setForceIteration(bool force_iteration)
Definition: newton.hh:415
Newton(const GridOperator &go, TrialVector &u_, Solver &solver_)
Definition: newton.hh:728
RFType first_defect
Definition: newton.hh:35
Definition: newton.hh:718
Strategy
Definition: newton.hh:561