8#ifndef DUNE_PDELAB_NEWTON_NEWTON_HH
9#define DUNE_PDELAB_NEWTON_NEWTON_HH
26#include <dune/pdelab/backend/solver.hh>
27#include <dune/pdelab/constraints/common/constraints.hh>
28#include <dune/pdelab/solver/newton.hh>
41 class NewtonLineSearchError :
public NewtonError {};
42 class NewtonNotConverged :
public NewtonError {};
45 template<
class RFType>
46 struct NewtonResult : LinearSolverResult<RFType>
50 double assembler_time;
51 double linear_solver_time;
52 int linear_solver_iterations;
55 first_defect(0.0), defect(0.0), assembler_time(0.0), linear_solver_time(0.0),
56 linear_solver_iterations(0) {}
59 template<
class GOS,
class TrlV,
class TstV>
62 typedef GOS GridOperator;
63 typedef TrlV TrialVector;
64 typedef TstV TestVector;
66 typedef typename TestVector::ElementType RFType;
67 typedef typename GOS::Traits::Jacobian Matrix;
72 typedef NewtonResult<RFType> Result;
74 void setVerbosityLevel(
unsigned int verbosity_level)
76 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
79 verbosity_level_ = verbosity_level;
83 void setKeepMatrix(
bool b)
89 bool keepMatrix()
const
102 const GridOperator& gridoperator_;
104 std::shared_ptr<TrialVector> z_;
105 std::shared_ptr<TestVector> r_;
106 std::shared_ptr<Matrix> A_;
108 unsigned int verbosity_level_;
110 RFType linear_reduction_;
116 NewtonBase(
const GridOperator& go, TrialVector& u)
119 , verbosity_level_(1)
122 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
123 verbosity_level_ = 0;
126 NewtonBase(
const GridOperator& go)
129 , verbosity_level_(1)
132 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
133 verbosity_level_ = 0;
136 virtual ~NewtonBase() { }
138 virtual bool terminate() = 0;
139 virtual void prepare_step(Matrix& A, TestVector& r) = 0;
140 virtual void line_search(TrialVector& z, TestVector& r) = 0;
141 virtual void defect(TestVector& r) = 0;
144 template<
class GOS,
class S,
class TrlV,
class TstV>
145 class NewtonSolver :
public virtual NewtonBase<GOS,TrlV,TstV>
148 typedef GOS GridOperator;
149 typedef TrlV TrialVector;
150 typedef TstV TestVector;
152 typedef typename TestVector::ElementType RFType;
153 typedef typename GOS::Traits::Jacobian Matrix;
156 typedef NewtonResult<RFType> Result;
158 NewtonSolver(
const GridOperator& go, TrialVector& u_, Solver& solver)
159 : NewtonBase<GOS,TrlV,TstV>(go,u_)
161 , result_valid_(false)
162 , use_maxnorm_(false)
165 NewtonSolver(
const GridOperator& go, Solver& solver)
166 : NewtonBase<GOS,TrlV,TstV>(go)
168 , result_valid_(false)
169 , use_maxnorm_(false)
170 , hanging_node_modifications_(false)
175 void apply(TrialVector& u_);
177 const Result& result()
const
181 "NewtonSolver::result() called before NewtonSolver::solve()");
186 void setUseMaxNorm(
bool b)
191 void setHangingNodeModifications(
bool b)
193 hanging_node_modifications_ = b;
197 virtual void defect(TestVector& r)
override
199 if (hanging_node_modifications_){
200 auto dirichletValues = *this->u_;
208 this->gridoperator_.localAssembler().backtransform(*this->u_);
212 this->gridoperator_.residual(*this->u_, r);
218 auto rank_max = Backend::native(r).infinity_norm();
219 this->res_.defect = this->gridoperator_.testGridFunctionSpace().gridView().comm().max(rank_max);
223 this->res_.defect = this->solver_.norm(r);
225 if (!std::isfinite(this->res_.defect))
227 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
232 void linearSolve(Matrix& A, TrialVector& z, TestVector& r)
const
234 if (this->verbosity_level_ >= 4)
235 std::cout <<
" Solving linear system..." << std::endl;
238 Impl::setLinearSystemReuse(this->solver_, !this->reassembled_);
239 this->solver_.apply(A, z, r, this->linear_reduction_);
241 ios_base_all_saver restorer(std::cout);
243 if (!this->solver_.result().converged)
245 "NewtonSolver::linearSolve(): Linear solver did not converge "
246 "in " << this->solver_.result().iterations <<
" iterations");
247 if (this->verbosity_level_ >= 4)
248 std::cout <<
" linear solver iterations: "
249 << std::setw(12) << solver_.result().iterations << std::endl
250 <<
" linear defect reduction: "
251 << std::setw(12) << std::setprecision(4) << std::scientific
252 << solver_.result().reduction << std::endl;
258 bool hanging_node_modifications_;
261 template<
class GOS,
class S,
class TrlV,
class TstV>
262 void NewtonSolver<GOS,S,TrlV,TstV>::apply(TrialVector& u)
268 template<
class GOS,
class S,
class TrlV,
class TstV>
269 void NewtonSolver<GOS,S,TrlV,TstV>::apply()
271 this->res_.iterations = 0;
272 this->res_.converged =
false;
273 this->res_.reduction = 1.0;
274 this->res_.conv_rate = 1.0;
275 this->res_.elapsed = 0.0;
276 this->res_.assembler_time = 0.0;
277 this->res_.linear_solver_time = 0.0;
278 this->res_.linear_solver_iterations = 0;
279 result_valid_ =
true;
286 this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
293 this->defect(*this->r_);
294 this->res_.first_defect = this->res_.defect;
295 this->prev_defect_ = this->res_.defect;
297 if (this->verbosity_level_ >= 2)
300 ios_base_all_saver restorer(std::cout);
301 std::cout <<
" Initial defect: "
302 << std::setw(12) << std::setprecision(4) << std::scientific
303 << this->res_.defect << std::endl;
308 this->A_ = std::make_shared<Matrix>(this->gridoperator_);
312 this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
315 while (!this->terminate())
317 if (this->verbosity_level_ >= 3)
318 std::cout <<
" Newton iteration " << this->res_.iterations
319 <<
" --------------------------------" << std::endl;
321 Timer assembler_timer;
330 this->prepare_step(*this->A_,*this->r_);
334 this->res_.assembler_time += assembler_timer.elapsed();
337 double assembler_time = assembler_timer.elapsed();
338 this->res_.assembler_time += assembler_time;
339 if (this->verbosity_level_ >= 3)
340 std::cout <<
" matrix assembly time: "
341 << std::setw(12) << std::setprecision(4) << std::scientific
342 << assembler_time << std::endl;
344 Timer linear_solver_timer;
351 this->linearSolve(*this->A_, *this->z_, *this->r_);
355 this->res_.linear_solver_time += linear_solver_timer.elapsed();
356 this->res_.linear_solver_iterations += this->solver_.result().iterations;
359 double linear_solver_time = linear_solver_timer.elapsed();
360 this->res_.linear_solver_time += linear_solver_time;
361 this->res_.linear_solver_iterations += this->solver_.result().iterations;
367 this->line_search(*this->z_, *this->r_);
369 catch (NewtonLineSearchError&)
371 if (this->reassembled_)
373 if (this->verbosity_level_ >= 3)
374 std::cout <<
" line search failed - trying again with reassembled matrix" << std::endl;
378 this->res_.reduction = this->res_.defect/this->res_.first_defect;
379 this->res_.iterations++;
380 this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
383 ios_base_all_saver restorer(std::cout);
385 if (this->verbosity_level_ >= 3)
386 std::cout <<
" linear solver time: "
387 << std::setw(12) << std::setprecision(4) << std::scientific
388 << linear_solver_time << std::endl
389 <<
" defect reduction (this iteration):"
390 << std::setw(12) << std::setprecision(4) << std::scientific
391 << this->res_.defect/this->prev_defect_ << std::endl
392 <<
" defect reduction (total): "
393 << std::setw(12) << std::setprecision(4) << std::scientific
394 << this->res_.reduction << std::endl
396 << std::setw(12) << std::setprecision(4) << std::scientific
397 << this->res_.defect << std::endl;
398 if (this->verbosity_level_ == 2)
399 std::cout <<
" Newton iteration " << std::setw(2) << this->res_.iterations
401 << std::setw(12) << std::setprecision(4) << std::scientific
403 <<
". Reduction (this): "
404 << std::setw(12) << std::setprecision(4) << std::scientific
405 << this->res_.defect/this->prev_defect_
406 <<
". Reduction (total): "
407 << std::setw(12) << std::setprecision(4) << std::scientific
408 << this->res_.reduction << std::endl;
413 this->res_.elapsed = timer.elapsed();
416 this->res_.elapsed = timer.elapsed();
418 ios_base_all_saver restorer(std::cout);
420 if (this->verbosity_level_ == 1)
421 std::cout <<
" Newton converged after " << std::setw(2) << this->res_.iterations
422 <<
" iterations. Reduction: "
423 << std::setw(12) << std::setprecision(4) << std::scientific
424 << this->res_.reduction
425 <<
" (" << std::setprecision(4) << this->res_.elapsed <<
"s)"
428 if(!this->keep_matrix_)
432 template<
class GOS,
class TrlV,
class TstV>
433 class NewtonTerminate :
public virtual NewtonBase<GOS,TrlV,TstV>
435 typedef GOS GridOperator;
436 typedef TrlV TrialVector;
438 typedef typename TstV::ElementType RFType;
441 NewtonTerminate(
const GridOperator& go, TrialVector& u_)
442 : NewtonBase<GOS,TrlV,TstV>(go,u_)
444 , force_iteration_(false)
446 this->reduction_ = 1e-8;
447 this->abs_limit_ = 1e-12;
450 NewtonTerminate(
const GridOperator& go)
451 : NewtonBase<GOS,TrlV,TstV>(go)
453 , force_iteration_(false)
455 this->reduction_ = 1e-8;
456 this->abs_limit_ = 1e-12;
459 void setReduction(RFType reduction)
461 this->reduction_ = reduction;
464 void setMaxIterations(
unsigned int maxit)
469 void setForceIteration(
bool force_iteration)
471 force_iteration_ = force_iteration;
474 void setAbsoluteLimit(RFType abs_limit_)
476 this->abs_limit_ = abs_limit_;
479 virtual bool terminate()
override
481 if (force_iteration_ && this->res_.iterations == 0)
483 this->res_.converged = this->res_.defect < this->abs_limit_
484 || this->res_.defect < this->res_.first_defect * this->reduction_;
485 if (this->res_.iterations >= maxit_ && !this->res_.converged)
487 "NewtonTerminate::terminate(): Maximum iteration count reached");
488 return this->res_.converged;
493 bool force_iteration_;
496 template<
class GOS,
class TrlV,
class TstV>
497 class NewtonPrepareStep :
public virtual NewtonBase<GOS,TrlV,TstV>
499 typedef GOS GridOperator;
500 typedef TrlV TrialVector;
502 typedef typename TstV::ElementType RFType;
503 typedef typename GOS::Traits::Jacobian Matrix;
506 NewtonPrepareStep(
const GridOperator& go, TrialVector& u_)
507 : NewtonBase<GOS,TrlV,TstV>(go,u_)
508 , min_linear_reduction_(1e-3)
509 , fixed_linear_reduction_(0.0)
510 , reassemble_threshold_(0.0)
513 NewtonPrepareStep(
const GridOperator& go)
514 : NewtonBase<GOS,TrlV,TstV>(go)
515 , min_linear_reduction_(1e-3)
516 , fixed_linear_reduction_(0.0)
517 , reassemble_threshold_(0.0)
518 , hanging_node_modifications_(false)
527 void setMinLinearReduction(RFType min_linear_reduction)
529 min_linear_reduction_ = min_linear_reduction;
536 void setFixedLinearReduction(
bool fixed_linear_reduction)
538 fixed_linear_reduction_ = fixed_linear_reduction;
548 void setReassembleThreshold(RFType reassemble_threshold)
550 reassemble_threshold_ = reassemble_threshold;
553 void setHangingNodeModifications(
bool b)
555 hanging_node_modifications_ = b;
558 virtual void prepare_step(Matrix& A, TstV& )
override
560 this->reassembled_ =
false;
561 if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
563 if (hanging_node_modifications_){
564 auto dirichletValues = *this->u_;
572 this->gridoperator_.localAssembler().backtransform(*this->u_);
575 if (this->verbosity_level_ >= 3)
576 std::cout <<
" Reassembling matrix..." << std::endl;
578 this->gridoperator_.jacobian(*this->u_, A);
579 this->reassembled_ =
true;
583 if (fixed_linear_reduction_ ==
true || (this->res_.iterations==0 && this->res_.first_defect<this->abs_limit_))
584 this->linear_reduction_ = min_linear_reduction_;
588 std::max(this->res_.first_defect * this->reduction_,
599 this->linear_reduction_ =
600 std::max( stop_defect/(10*this->res_.defect),
601 std::min(min_linear_reduction_,this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_)) );
604 this->prev_defect_ = this->res_.defect;
606 ios_base_all_saver restorer(std::cout);
608 if (this->verbosity_level_ >= 3)
609 std::cout <<
" requested linear reduction: "
610 << std::setw(12) << std::setprecision(4) << std::scientific
611 << this->linear_reduction_ << std::endl;
615 RFType min_linear_reduction_;
616 bool fixed_linear_reduction_;
617 RFType reassemble_threshold_;
618 bool hanging_node_modifications_;
621 template<
class GOS,
class TrlV,
class TstV>
622 class NewtonLineSearch :
public virtual NewtonBase<GOS,TrlV,TstV>
624 typedef GOS GridOperator;
625 typedef TrlV TrialVector;
626 typedef TstV TestVector;
628 typedef typename TestVector::ElementType RFType;
639 hackbuschReuskenAcceptBest };
641 NewtonLineSearch(
const GridOperator& go, TrialVector& u_)
642 : NewtonBase<GOS,TrlV,TstV>(go,u_)
643 , strategy_(hackbuschReusken)
645 , damping_factor_(0.5)
648 NewtonLineSearch(
const GridOperator& go)
649 : NewtonBase<GOS,TrlV,TstV>(go)
650 , strategy_(hackbuschReusken)
652 , damping_factor_(0.5)
655 void setLineSearchStrategy(Strategy strategy)
657 strategy_ = strategy;
660 void setLineSearchStrategy(std::string strategy)
662 strategy_ = strategyFromName(strategy);
665 void setLineSearchMaxIterations(
unsigned int maxit)
670 void setLineSearchDampingFactor(RFType damping_factor)
672 damping_factor_ = damping_factor;
675 virtual void line_search(TrialVector& z, TestVector& r)
override
677 if ((strategy_ == noLineSearch) || (this->res_.defect < this->abs_limit_))
679 this->u_->axpy(-1.0, z);
684 if (this->verbosity_level_ >= 4)
685 std::cout <<
" Performing line search..." << std::endl;
687 RFType best_lambda = 0.0;
688 RFType best_defect = this->res_.defect;
689 TrialVector prev_u(*this->u_);
691 ios_base_all_saver restorer(std::cout);
695 if (this->verbosity_level_ >= 4)
696 std::cout <<
" trying line search damping factor: "
697 << std::setw(12) << std::setprecision(4) << std::scientific
701 this->u_->axpy(-lambda, z);
705 catch (NewtonDefectError&)
707 if (this->verbosity_level_ >= 4)
708 std::cout <<
" NaNs detected" << std::endl;
711 if (this->res_.defect <= (1.0 - lambda/4) * this->prev_defect_)
713 if (this->verbosity_level_ >= 4)
714 std::cout <<
" line search converged" << std::endl;
718 if (this->res_.defect < best_defect)
720 best_defect = this->res_.defect;
721 best_lambda = lambda;
726 if (this->verbosity_level_ >= 4)
727 std::cout <<
" max line search iterations exceeded" << std::endl;
730 case hackbuschReusken:
734 "NewtonLineSearch::line_search(): line search failed, "
735 "max iteration count reached, "
736 "defect did not improve enough");
737 case hackbuschReuskenAcceptBest:
738 if (best_lambda == 0.0)
743 "NewtonLineSearch::line_search(): line search failed, "
744 "max iteration count reached, "
745 "defect did not improve in any of the iterations");
747 if (best_lambda != lambda)
750 this->u_->axpy(-best_lambda, z);
760 lambda *= damping_factor_;
763 if (this->verbosity_level_ >= 4)
764 std::cout <<
" line search damping factor: "
765 << std::setw(12) << std::setprecision(4) << std::scientific
766 << lambda << std::endl;
771 Strategy strategyFromName(
const std::string & s) {
772 if (s ==
"noLineSearch")
774 if (s ==
"hackbuschReusken")
775 return hackbuschReusken;
776 if (s ==
"hackbuschReuskenAcceptBest")
777 return hackbuschReuskenAcceptBest;
778 DUNE_THROW(Exception,
"unknown line search strategy" << s);
784 RFType damping_factor_;
787 template<
class GOS,
class S,
class TrlV,
class TstV = TrlV>
788 class Newton :
public NewtonSolver<GOS,S,TrlV,TstV>
789 ,
public NewtonTerminate<GOS,TrlV,TstV>
790 ,
public NewtonLineSearch<GOS,TrlV,TstV>
791 ,
public NewtonPrepareStep<GOS,TrlV,TstV>
793 typedef GOS GridOperator;
795 typedef TrlV TrialVector;
798 [[deprecated(
"This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")]]
799 Newton(
const GridOperator& go, TrialVector& u_, Solver& solver_)
800 : NewtonBase<GOS,TrlV,TstV>(go,u_)
801 , NewtonSolver<GOS,S,TrlV,TstV>(go,u_,solver_)
802 , NewtonTerminate<GOS,TrlV,TstV>(go,u_)
803 , NewtonLineSearch<GOS,TrlV,TstV>(go,u_)
804 , NewtonPrepareStep<GOS,TrlV,TstV>(go,u_)
806 [[deprecated(
"This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")]]
807 Newton(
const GridOperator& go, Solver& solver_)
808 : NewtonBase<GOS,TrlV,TstV>(go)
809 , NewtonSolver<GOS,S,TrlV,TstV>(go,solver_)
810 , NewtonTerminate<GOS,TrlV,TstV>(go)
811 , NewtonLineSearch<GOS,TrlV,TstV>(go)
812 , NewtonPrepareStep<GOS,TrlV,TstV>(go)
838 typedef typename TstV::ElementType RFType;
839 if (param.
hasKey(
"VerbosityLevel"))
840 this->setVerbosityLevel(
841 param.
get<
unsigned int>(
"VerbosityLevel"));
842 if (param.
hasKey(
"Reduction"))
844 param.
get<RFType>(
"Reduction"));
845 if (param.
hasKey(
"MaxIterations"))
846 this->setMaxIterations(
847 param.
get<
unsigned int>(
"MaxIterations"));
848 if (param.
hasKey(
"ForceIteration"))
849 this->setForceIteration(
850 param.
get<
bool>(
"ForceIteration"));
851 if (param.
hasKey(
"AbsoluteLimit"))
852 this->setAbsoluteLimit(
853 param.
get<RFType>(
"AbsoluteLimit"));
854 if (param.
hasKey(
"MinLinearReduction"))
855 this->setMinLinearReduction(
856 param.
get<RFType>(
"MinLinearReduction"));
857 if (param.
hasKey(
"FixedLinearReduction"))
858 this->setFixedLinearReduction(
859 param.
get<
bool>(
"FixedLinearReduction"));
860 if (param.
hasKey(
"ReassembleThreshold"))
861 this->setReassembleThreshold(
862 param.
get<RFType>(
"ReassembleThreshold"));
863 if (param.
hasKey(
"LineSearchStrategy"))
864 this->setLineSearchStrategy(
865 param.
get<std::string>(
"LineSearchStrategy"));
866 if (param.
hasKey(
"LineSearchMaxIterations"))
867 this->setLineSearchMaxIterations(
868 param.
get<
unsigned int>(
"LineSearchMaxIterations"));
869 if (param.
hasKey(
"LineSearchDampingFactor"))
870 this->setLineSearchDampingFactor(
871 param.
get<RFType>(
"LineSearchDampingFactor"));
872 if (param.
hasKey(
"KeepMatrix"))
874 param.
get<
bool>(
"KeepMatrix"));
875 if (param.
hasKey(
"UseMaxNorm"))
877 param.
get<
bool>(
"UseMaxNorm"));
Hierarchical structure of string parameters.
Definition: parametertree.hh:37
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:188
bool hasKey(const std::string &key) const
test for key
Definition: parametertree.cc:51
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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
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
Utility class for storing and resetting stream attributes.
Dune namespace.
Definition: alignedallocator.hh:13
A hierarchical structure of string parameters.
Standard Dune debug streams.