8#ifndef DUNE_PDELAB_NEWTON_NEWTON_HH
9#define DUNE_PDELAB_NEWTON_NEWTON_HH
20#include <dune/common/stdstreams.hh>
21#include <dune/common/exceptions.hh>
22#include <dune/common/ios_state.hh>
23#include <dune/common/timer.hh>
24#include <dune/common/parametertree.hh>
40 template<
class RFType>
41 struct NewtonResult : LinearSolverResult<RFType>
45 double assembler_time;
46 double linear_solver_time;
47 int linear_solver_iterations;
50 first_defect(0.0), defect(0.0), assembler_time(0.0), linear_solver_time(0.0),
51 linear_solver_iterations(0) {}
54 template<
class GOS,
class TrlV,
class TstV>
57 typedef GOS GridOperator;
58 typedef TrlV TrialVector;
59 typedef TstV TestVector;
61 typedef typename TestVector::ElementType RFType;
62 typedef typename GOS::Traits::Jacobian
Matrix;
67 typedef NewtonResult<RFType> Result;
69 void setVerbosityLevel(
unsigned int verbosity_level)
71 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
74 verbosity_level_ = verbosity_level;
78 void setKeepMatrix(
bool b)
84 bool keepMatrix()
const
97 const GridOperator& gridoperator_;
99 std::shared_ptr<TrialVector> z_;
100 std::shared_ptr<TestVector> r_;
101 std::shared_ptr<Matrix> A_;
103 unsigned int verbosity_level_;
105 RFType linear_reduction_;
111 NewtonBase(
const GridOperator& go, TrialVector& u)
114 , verbosity_level_(1)
117 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
118 verbosity_level_ = 0;
121 NewtonBase(
const GridOperator& go)
124 , verbosity_level_(1)
127 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
128 verbosity_level_ = 0;
131 virtual ~NewtonBase() { }
133 virtual bool terminate() = 0;
134 virtual void prepare_step(Matrix& A, TestVector& r) = 0;
135 virtual void line_search(TrialVector& z, TestVector& r) = 0;
136 virtual void defect(TestVector& r) = 0;
139 template<
class GOS,
class S,
class TrlV,
class TstV>
140 class NewtonSolver :
public virtual NewtonBase<GOS,TrlV,TstV>
143 typedef GOS GridOperator;
144 typedef TrlV TrialVector;
145 typedef TstV TestVector;
147 typedef typename TestVector::ElementType RFType;
148 typedef typename GOS::Traits::Jacobian
Matrix;
151 typedef NewtonResult<RFType> Result;
153 NewtonSolver(
const GridOperator& go, TrialVector& u_, Solver& solver)
154 : NewtonBase<GOS,TrlV,TstV>(go,u_)
156 , result_valid_(false)
157 , use_maxnorm_(false)
160 NewtonSolver(
const GridOperator& go, Solver& solver)
161 : NewtonBase<GOS,TrlV,TstV>(go)
163 , result_valid_(false)
164 , use_maxnorm_(false)
165 , hanging_node_modifications_(false)
170 void apply(TrialVector& u_);
172 const Result& result()
const
175 DUNE_THROW(NewtonError,
176 "NewtonSolver::result() called before NewtonSolver::solve()");
181 void setUseMaxNorm(
bool b)
186 void setHangingNodeModifications(
bool b)
188 hanging_node_modifications_ = b;
192 virtual void defect(TestVector& r)
override
194 if (hanging_node_modifications_){
195 auto dirichletValues = *this->u_;
203 this->gridoperator_.localAssembler().backtransform(*this->u_);
207 this->gridoperator_.residual(*this->u_, r);
213 auto rank_max = Backend::native(r).infinity_norm();
214 this->res_.defect = this->gridoperator_.testGridFunctionSpace().gridView().comm().max(rank_max);
218 this->res_.defect = this->solver_.norm(r);
220 if (!std::isfinite(this->res_.defect))
221 DUNE_THROW(NewtonDefectError,
222 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
227 void linearSolve(Matrix& A, TrialVector& z, TestVector& r)
const
229 if (this->verbosity_level_ >= 4)
230 std::cout <<
" Solving linear system..." << std::endl;
233 Impl::setLinearSystemReuse(this->solver_, !this->reassembled_);
234 this->solver_.apply(A, z, r, this->linear_reduction_);
236 ios_base_all_saver restorer(std::cout);
238 if (!this->solver_.result().converged)
239 DUNE_THROW(NewtonLinearSolverError,
240 "NewtonSolver::linearSolve(): Linear solver did not converge "
241 "in " << this->solver_.result().iterations <<
" iterations");
242 if (this->verbosity_level_ >= 4)
243 std::cout <<
" linear solver iterations: "
244 << std::setw(12) << solver_.result().iterations << std::endl
245 <<
" linear defect reduction: "
246 << std::setw(12) << std::setprecision(4) << std::scientific
247 << solver_.result().reduction << std::endl;
253 bool hanging_node_modifications_;
256 template<
class GOS,
class S,
class TrlV,
class TstV>
257 void NewtonSolver<GOS,S,TrlV,TstV>::apply(TrialVector& u)
263 template<
class GOS,
class S,
class TrlV,
class TstV>
264 void NewtonSolver<GOS,S,TrlV,TstV>::apply()
266 this->res_.iterations = 0;
267 this->res_.converged =
false;
268 this->res_.reduction = 1.0;
269 this->res_.conv_rate = 1.0;
270 this->res_.elapsed = 0.0;
271 this->res_.assembler_time = 0.0;
272 this->res_.linear_solver_time = 0.0;
273 this->res_.linear_solver_iterations = 0;
274 result_valid_ =
true;
281 this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
288 this->defect(*this->r_);
289 this->res_.first_defect = this->res_.defect;
290 this->prev_defect_ = this->res_.defect;
292 if (this->verbosity_level_ >= 2)
295 ios_base_all_saver restorer(std::cout);
296 std::cout <<
" Initial defect: "
297 << std::setw(12) << std::setprecision(4) << std::scientific
298 << this->res_.defect << std::endl;
303 this->A_ = std::make_shared<Matrix>(this->gridoperator_);
307 this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
310 while (!this->terminate())
312 if (this->verbosity_level_ >= 3)
313 std::cout <<
" Newton iteration " << this->res_.iterations
314 <<
" --------------------------------" << std::endl;
316 Timer assembler_timer;
325 this->prepare_step(*this->A_,*this->r_);
329 this->res_.assembler_time += assembler_timer.elapsed();
332 double assembler_time = assembler_timer.elapsed();
333 this->res_.assembler_time += assembler_time;
334 if (this->verbosity_level_ >= 3)
335 std::cout <<
" matrix assembly time: "
336 << std::setw(12) << std::setprecision(4) << std::scientific
337 << assembler_time << std::endl;
339 Timer linear_solver_timer;
346 this->linearSolve(*this->A_, *this->z_, *this->r_);
350 this->res_.linear_solver_time += linear_solver_timer.elapsed();
351 this->res_.linear_solver_iterations += this->solver_.result().iterations;
354 double linear_solver_time = linear_solver_timer.elapsed();
355 this->res_.linear_solver_time += linear_solver_time;
356 this->res_.linear_solver_iterations += this->solver_.result().iterations;
362 this->line_search(*this->z_, *this->r_);
364 catch (NewtonLineSearchError&)
366 if (this->reassembled_)
368 if (this->verbosity_level_ >= 3)
369 std::cout <<
" line search failed - trying again with reassembled matrix" << std::endl;
373 this->res_.reduction = this->res_.defect/this->res_.first_defect;
374 this->res_.iterations++;
375 this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
378 ios_base_all_saver restorer(std::cout);
380 if (this->verbosity_level_ >= 3)
381 std::cout <<
" linear solver time: "
382 << std::setw(12) << std::setprecision(4) << std::scientific
383 << linear_solver_time << std::endl
384 <<
" defect reduction (this iteration):"
385 << std::setw(12) << std::setprecision(4) << std::scientific
386 << this->res_.defect/this->prev_defect_ << std::endl
387 <<
" defect reduction (total): "
388 << std::setw(12) << std::setprecision(4) << std::scientific
389 << this->res_.reduction << std::endl
391 << std::setw(12) << std::setprecision(4) << std::scientific
392 << this->res_.defect << std::endl;
393 if (this->verbosity_level_ == 2)
394 std::cout <<
" Newton iteration " << std::setw(2) << this->res_.iterations
396 << std::setw(12) << std::setprecision(4) << std::scientific
398 <<
". Reduction (this): "
399 << std::setw(12) << std::setprecision(4) << std::scientific
400 << this->res_.defect/this->prev_defect_
401 <<
". Reduction (total): "
402 << std::setw(12) << std::setprecision(4) << std::scientific
403 << this->res_.reduction << std::endl;
408 this->res_.elapsed = timer.elapsed();
411 this->res_.elapsed = timer.elapsed();
413 ios_base_all_saver restorer(std::cout);
415 if (this->verbosity_level_ == 1)
416 std::cout <<
" Newton converged after " << std::setw(2) << this->res_.iterations
417 <<
" iterations. Reduction: "
418 << std::setw(12) << std::setprecision(4) << std::scientific
419 << this->res_.reduction
420 <<
" (" << std::setprecision(4) << this->res_.elapsed <<
"s)"
423 if(!this->keep_matrix_)
427 template<
class GOS,
class TrlV,
class TstV>
428 class NewtonTerminate :
public virtual NewtonBase<GOS,TrlV,TstV>
430 typedef GOS GridOperator;
431 typedef TrlV TrialVector;
433 typedef typename TstV::ElementType RFType;
436 NewtonTerminate(
const GridOperator& go, TrialVector& u_)
437 : NewtonBase<GOS,TrlV,TstV>(go,u_)
439 , force_iteration_(false)
441 this->reduction_ = 1e-8;
442 this->abs_limit_ = 1e-12;
445 NewtonTerminate(
const GridOperator& go)
446 : NewtonBase<GOS,TrlV,TstV>(go)
448 , force_iteration_(false)
450 this->reduction_ = 1e-8;
451 this->abs_limit_ = 1e-12;
454 void setReduction(RFType reduction)
456 this->reduction_ = reduction;
459 void setMaxIterations(
unsigned int maxit)
464 void setForceIteration(
bool force_iteration)
466 force_iteration_ = force_iteration;
469 void setAbsoluteLimit(RFType abs_limit_)
471 this->abs_limit_ = abs_limit_;
474 virtual bool terminate()
override
476 if (force_iteration_ && this->res_.iterations == 0)
478 this->res_.converged = this->res_.defect < this->abs_limit_
479 || this->res_.defect < this->res_.first_defect * this->reduction_;
480 if (this->res_.iterations >= maxit_ && !this->res_.converged)
481 DUNE_THROW(NewtonNotConverged,
482 "NewtonTerminate::terminate(): Maximum iteration count reached");
483 return this->res_.converged;
488 bool force_iteration_;
491 template<
class GOS,
class TrlV,
class TstV>
492 class NewtonPrepareStep :
public virtual NewtonBase<GOS,TrlV,TstV>
494 typedef GOS GridOperator;
495 typedef TrlV TrialVector;
497 typedef typename TstV::ElementType RFType;
498 typedef typename GOS::Traits::Jacobian
Matrix;
501 NewtonPrepareStep(
const GridOperator& go, TrialVector& u_)
502 : NewtonBase<GOS,TrlV,TstV>(go,u_)
503 , min_linear_reduction_(1e-3)
504 , fixed_linear_reduction_(0.0)
505 , reassemble_threshold_(0.0)
508 NewtonPrepareStep(
const GridOperator& go)
509 : NewtonBase<GOS,TrlV,TstV>(go)
510 , min_linear_reduction_(1e-3)
511 , fixed_linear_reduction_(0.0)
512 , reassemble_threshold_(0.0)
513 , hanging_node_modifications_(false)
522 void setMinLinearReduction(RFType min_linear_reduction)
524 min_linear_reduction_ = min_linear_reduction;
531 void setFixedLinearReduction(
bool fixed_linear_reduction)
533 fixed_linear_reduction_ = fixed_linear_reduction;
543 void setReassembleThreshold(RFType reassemble_threshold)
545 reassemble_threshold_ = reassemble_threshold;
548 void setHangingNodeModifications(
bool b)
550 hanging_node_modifications_ = b;
553 virtual void prepare_step(Matrix& A, TstV& )
override
555 this->reassembled_ =
false;
556 if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
558 if (hanging_node_modifications_){
559 auto dirichletValues = *this->u_;
567 this->gridoperator_.localAssembler().backtransform(*this->u_);
570 if (this->verbosity_level_ >= 3)
571 std::cout <<
" Reassembling matrix..." << std::endl;
573 this->gridoperator_.jacobian(*this->u_, A);
574 this->reassembled_ =
true;
577 if (fixed_linear_reduction_ ==
true)
578 this->linear_reduction_ = min_linear_reduction_;
582 std::max(this->res_.first_defect * this->reduction_,
593 if ( stop_defect/(10*this->res_.defect) >
594 this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_) )
595 this->linear_reduction_ =
596 stop_defect/(10*this->res_.defect);
598 this->linear_reduction_ =
599 std::min(min_linear_reduction_,this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_));
602 this->prev_defect_ = this->res_.defect;
604 ios_base_all_saver restorer(std::cout);
606 if (this->verbosity_level_ >= 3)
607 std::cout <<
" requested linear reduction: "
608 << std::setw(12) << std::setprecision(4) << std::scientific
609 << this->linear_reduction_ << std::endl;
613 RFType min_linear_reduction_;
614 bool fixed_linear_reduction_;
615 RFType reassemble_threshold_;
616 bool hanging_node_modifications_;
619 template<
class GOS,
class TrlV,
class TstV>
620 class NewtonLineSearch :
public virtual NewtonBase<GOS,TrlV,TstV>
622 typedef GOS GridOperator;
623 typedef TrlV TrialVector;
624 typedef TstV TestVector;
626 typedef typename TestVector::ElementType RFType;
637 hackbuschReuskenAcceptBest };
639 NewtonLineSearch(
const GridOperator& go, TrialVector& u_)
640 : NewtonBase<GOS,TrlV,TstV>(go,u_)
643 , damping_factor_(0.5)
646 NewtonLineSearch(
const GridOperator& go)
647 : NewtonBase<GOS,TrlV,TstV>(go)
650 , damping_factor_(0.5)
653 void setLineSearchStrategy(Strategy strategy)
655 strategy_ = strategy;
658 void setLineSearchStrategy(std::string strategy)
660 strategy_ = strategyFromName(strategy);
663 void setLineSearchMaxIterations(
unsigned int maxit)
668 void setLineSearchDampingFactor(RFType damping_factor)
670 damping_factor_ = damping_factor;
673 virtual void line_search(TrialVector& z, TestVector& r)
override
675 if ((strategy_ == noLineSearch) || (this->res_.defect < this->abs_limit_))
677 this->u_->axpy(-1.0, z);
682 if (this->verbosity_level_ >= 4)
683 std::cout <<
" Performing line search..." << std::endl;
685 RFType best_lambda = 0.0;
686 RFType best_defect = this->res_.defect;
687 TrialVector prev_u(*this->u_);
689 ios_base_all_saver restorer(std::cout);
693 if (this->verbosity_level_ >= 4)
694 std::cout <<
" trying line search damping factor: "
695 << std::setw(12) << std::setprecision(4) << std::scientific
699 this->u_->axpy(-lambda, z);
703 catch (NewtonDefectError&)
705 if (this->verbosity_level_ >= 4)
706 std::cout <<
" NaNs detected" << std::endl;
709 if (this->res_.defect <= (1.0 - lambda/4) * this->prev_defect_)
711 if (this->verbosity_level_ >= 4)
712 std::cout <<
" line search converged" << std::endl;
716 if (this->res_.defect < best_defect)
718 best_defect = this->res_.defect;
719 best_lambda = lambda;
724 if (this->verbosity_level_ >= 4)
725 std::cout <<
" max line search iterations exceeded" << std::endl;
731 DUNE_THROW(NewtonLineSearchError,
732 "NewtonLineSearch::line_search(): line search failed, "
733 "max iteration count reached, "
734 "defect did not improve enough");
735 case hackbuschReuskenAcceptBest:
736 if (best_lambda == 0.0)
740 DUNE_THROW(NewtonLineSearchError,
741 "NewtonLineSearch::line_search(): line search failed, "
742 "max iteration count reached, "
743 "defect did not improve in any of the iterations");
745 if (best_lambda != lambda)
748 this->u_->axpy(-best_lambda, z);
758 lambda *= damping_factor_;
761 if (this->verbosity_level_ >= 4)
762 std::cout <<
" line search damping factor: "
763 << std::setw(12) << std::setprecision(4) << std::scientific
764 << lambda << std::endl;
769 Strategy strategyFromName(
const std::string &
s) {
770 if (
s ==
"noLineSearch")
772 if (
s ==
"hackbuschReusken")
774 if (
s ==
"hackbuschReuskenAcceptBest")
775 return hackbuschReuskenAcceptBest;
776 DUNE_THROW(Exception,
"unknown line search strategy" <<
s);
782 RFType damping_factor_;
785 template<
class GOS,
class S,
class TrlV,
class TstV = TrlV>
786 class Newton :
public NewtonSolver<GOS,S,TrlV,TstV>
787 ,
public NewtonTerminate<GOS,TrlV,TstV>
788 ,
public NewtonLineSearch<GOS,TrlV,TstV>
789 ,
public NewtonPrepareStep<GOS,TrlV,TstV>
791 typedef GOS GridOperator;
793 typedef TrlV TrialVector;
796 Newton(
const GridOperator& go, TrialVector& u_, Solver& solver_)
797 DUNE_DEPRECATED_MSG(
"This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")
798 : NewtonBase<GOS,TrlV,TstV>(go,u_)
799 , NewtonSolver<GOS,S,TrlV,TstV>(go,u_,solver_)
800 , NewtonTerminate<GOS,TrlV,TstV>(go,u_)
801 , NewtonLineSearch<GOS,TrlV,TstV>(go,u_)
802 , NewtonPrepareStep<GOS,TrlV,TstV>(go,u_)
804 Newton(
const GridOperator& go, Solver& solver_)
805 DUNE_DEPRECATED_MSG(
"This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")
806 : NewtonBase<GOS,TrlV,TstV>(go)
807 , NewtonSolver<GOS,S,TrlV,TstV>(go,solver_)
808 , NewtonTerminate<GOS,TrlV,TstV>(go)
809 , NewtonLineSearch<GOS,TrlV,TstV>(go)
810 , NewtonPrepareStep<GOS,TrlV,TstV>(go)
834 void setParameters(
const Dune::ParameterTree & param)
836 typedef typename TstV::ElementType RFType;
837 if (param.hasKey(
"VerbosityLevel"))
838 this->setVerbosityLevel(
839 param.get<
unsigned int>(
"VerbosityLevel"));
840 if (param.hasKey(
"Reduction"))
842 param.get<RFType>(
"Reduction"));
843 if (param.hasKey(
"MaxIterations"))
844 this->setMaxIterations(
845 param.get<
unsigned int>(
"MaxIterations"));
846 if (param.hasKey(
"ForceIteration"))
847 this->setForceIteration(
848 param.get<
bool>(
"ForceIteration"));
849 if (param.hasKey(
"AbsoluteLimit"))
850 this->setAbsoluteLimit(
851 param.get<RFType>(
"AbsoluteLimit"));
852 if (param.hasKey(
"MinLinearReduction"))
853 this->setMinLinearReduction(
854 param.get<RFType>(
"MinLinearReduction"));
855 if (param.hasKey(
"FixedLinearReduction"))
856 this->setFixedLinearReduction(
857 param.get<
bool>(
"FixedLinearReduction"));
858 if (param.hasKey(
"ReassembleThreshold"))
859 this->setReassembleThreshold(
860 param.get<RFType>(
"ReassembleThreshold"));
861 if (param.hasKey(
"LineSearchStrategy"))
862 this->setLineSearchStrategy(
863 param.get<std::string>(
"LineSearchStrategy"));
864 if (param.hasKey(
"LineSearchMaxIterations"))
865 this->setLineSearchMaxIterations(
866 param.get<
unsigned int>(
"LineSearchMaxIterations"));
867 if (param.hasKey(
"LineSearchDampingFactor"))
868 this->setLineSearchDampingFactor(
869 param.get<RFType>(
"LineSearchDampingFactor"));
870 if (param.hasKey(
"KeepMatrix"))
872 param.get<
bool>(
"KeepMatrix"));
873 if (param.hasKey(
"UseMaxNorm"))
875 param.get<
bool>(
"UseMaxNorm"));
const std::string s
Definition function.hh:843
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition constraints.hh:796
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition constraints.hh:1014
void copy_constrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition constraints.hh:936
For backward compatibility – Do not use this!
Definition adaptivity.hh:28
@ noLineSearch
Definition newtonlinesearch.hh:192
@ hackbuschReusken
Definition newtonlinesearch.hh:193
typename impl::BackendMatrixSelector< Backend, VU, VV, E >::Type Matrix
alias of the return type of BackendMatrixSelector
Definition backend/interface.hh:127