3 #ifndef DUNE_PDELAB_NEWTON_NEWTON_HH 4 #define DUNE_PDELAB_NEWTON_NEWTON_HH 11 #include <type_traits> 15 #include <dune/common/stdstreams.hh> 16 #include <dune/common/exceptions.hh> 17 #include <dune/common/ios_state.hh> 18 #include <dune/common/timer.hh> 19 #include <dune/common/parametertree.hh> 30 template<
typename T1,
typename =
void>
36 struct HasSetReuse<T, decltype(
std::declval<T>().setReuse(true), void())>
41 inline void setLinearSystemReuse(T& solver_backend,
bool reuse, std::true_type)
43 if (!solver_backend.getReuse() && reuse)
44 dwarn <<
"WARNING: Newton needed to override your choice to reuse the linear system in order to work!" << std::endl;
45 solver_backend.setReuse(reuse);
49 inline void setLinearSystemReuse(T&,
bool, std::false_type)
53 inline void setLinearSystemReuse(T& solver_backend,
bool reuse)
55 setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
67 template<
class RFType>
77 first_defect(0.0), defect(0.0), assembler_time(0.0), linear_solver_time(0.0),
78 linear_solver_iterations(0) {}
81 template<
class GOS,
class TrlV,
class TstV>
84 typedef GOS GridOperator;
85 typedef TrlV TrialVector;
86 typedef TstV TestVector;
88 typedef typename TestVector::ElementType RFType;
89 typedef typename GOS::Traits::Jacobian
Matrix;
98 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
101 verbosity_level_ = verbosity_level;
126 std::shared_ptr<TrialVector>
z_;
127 std::shared_ptr<TestVector>
r_;
128 std::shared_ptr<Matrix>
A_;
141 , verbosity_level_(1)
144 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
145 verbosity_level_ = 0;
151 , verbosity_level_(1)
154 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
155 verbosity_level_ = 0;
160 virtual bool terminate() = 0;
161 virtual void prepare_step(Matrix& A, TestVector& r) = 0;
162 virtual void line_search(TrialVector& z, TestVector& r) = 0;
163 virtual void defect(TestVector& r) = 0;
166 template<
class GOS,
class S,
class TrlV,
class TstV>
170 typedef GOS GridOperator;
171 typedef TrlV TrialVector;
172 typedef TstV TestVector;
174 typedef typename TestVector::ElementType RFType;
175 typedef typename GOS::Traits::Jacobian
Matrix;
183 , result_valid_(false)
189 , result_valid_(false)
194 void apply(TrialVector& u_);
200 "NewtonSolver::result() called before NewtonSolver::solve()");
205 virtual void defect(TestVector& r)
override 208 this->gridoperator_.residual(*this->u_, r);
209 this->res_.defect = this->solver_.norm(r);
210 if (!std::isfinite(this->res_.defect))
212 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
217 void linearSolve(Matrix& A, TrialVector& z, TestVector& r)
const 219 if (this->verbosity_level_ >= 4)
220 std::cout <<
" Solving linear system..." << std::endl;
223 Impl::setLinearSystemReuse(this->solver_, !this->reassembled_);
224 this->solver_.apply(A, z, r, this->linear_reduction_);
226 ios_base_all_saver restorer(std::cout);
228 if (!this->solver_.result().converged)
230 "NewtonSolver::linearSolve(): Linear solver did not converge " 231 "in " << this->solver_.result().iterations <<
" iterations");
232 if (this->verbosity_level_ >= 4)
233 std::cout <<
" linear solver iterations: " 234 << std::setw(12) << solver_.result().iterations << std::endl
235 <<
" linear defect reduction: " 236 << std::setw(12) << std::setprecision(4) << std::scientific
237 << solver_.result().reduction << std::endl;
244 template<
class GOS,
class S,
class TrlV,
class TstV>
251 template<
class GOS,
class S,
class TrlV,
class TstV>
254 this->res_.iterations = 0;
255 this->res_.converged =
false;
256 this->res_.reduction = 1.0;
257 this->res_.conv_rate = 1.0;
258 this->res_.elapsed = 0.0;
259 this->res_.assembler_time = 0.0;
260 this->res_.linear_solver_time = 0.0;
261 this->res_.linear_solver_iterations = 0;
262 result_valid_ =
true;
269 this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
276 this->defect(*this->r_);
277 this->res_.first_defect = this->res_.defect;
278 this->prev_defect_ = this->res_.defect;
280 if (this->verbosity_level_ >= 2)
283 ios_base_all_saver restorer(std::cout);
284 std::cout <<
" Initial defect: " 285 << std::setw(12) << std::setprecision(4) << std::scientific
286 << this->res_.defect << std::endl;
291 this->A_ = std::make_shared<Matrix>(this->gridoperator_);
295 this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
298 while (!this->terminate())
300 if (this->verbosity_level_ >= 3)
301 std::cout <<
" Newton iteration " << this->res_.iterations
302 <<
" --------------------------------" << std::endl;
304 Timer assembler_timer;
313 this->prepare_step(*this->A_,*this->r_);
317 this->res_.assembler_time += assembler_timer.elapsed();
320 double assembler_time = assembler_timer.elapsed();
321 this->res_.assembler_time += assembler_time;
322 if (this->verbosity_level_ >= 3)
323 std::cout <<
" matrix assembly time: " 324 << std::setw(12) << std::setprecision(4) << std::scientific
325 << assembler_time << std::endl;
327 Timer linear_solver_timer;
334 this->linearSolve(*this->A_, *this->z_, *this->r_);
338 this->res_.linear_solver_time += linear_solver_timer.elapsed();
339 this->res_.linear_solver_iterations += this->solver_.result().iterations;
342 double linear_solver_time = linear_solver_timer.elapsed();
343 this->res_.linear_solver_time += linear_solver_time;
344 this->res_.linear_solver_iterations += this->solver_.result().iterations;
350 this->line_search(*this->z_, *this->r_);
354 if (this->reassembled_)
356 if (this->verbosity_level_ >= 3)
357 std::cout <<
" line search failed - trying again with reassembled matrix" << std::endl;
361 this->res_.reduction = this->res_.defect/this->res_.first_defect;
362 this->res_.iterations++;
363 this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
366 ios_base_all_saver restorer(std::cout);
368 if (this->verbosity_level_ >= 3)
369 std::cout <<
" linear solver time: " 370 << std::setw(12) << std::setprecision(4) << std::scientific
371 << linear_solver_time << std::endl
372 <<
" defect reduction (this iteration):" 373 << std::setw(12) << std::setprecision(4) << std::scientific
374 << this->res_.defect/this->prev_defect_ << std::endl
375 <<
" defect reduction (total): " 376 << std::setw(12) << std::setprecision(4) << std::scientific
377 << this->res_.reduction << std::endl
379 << std::setw(12) << std::setprecision(4) << std::scientific
380 << this->res_.defect << std::endl;
381 if (this->verbosity_level_ == 2)
382 std::cout <<
" Newton iteration " << std::setw(2) << this->res_.iterations
384 << std::setw(12) << std::setprecision(4) << std::scientific
386 <<
". Reduction (this): " 387 << std::setw(12) << std::setprecision(4) << std::scientific
388 << this->res_.defect/this->prev_defect_
389 <<
". Reduction (total): " 390 << std::setw(12) << std::setprecision(4) << std::scientific
391 << this->res_.reduction << std::endl;
396 this->res_.elapsed = timer.elapsed();
399 this->res_.elapsed = timer.elapsed();
401 ios_base_all_saver restorer(std::cout);
403 if (this->verbosity_level_ == 1)
404 std::cout <<
" Newton converged after " << std::setw(2) << this->res_.iterations
405 <<
" iterations. Reduction: " 406 << std::setw(12) << std::setprecision(4) << std::scientific
407 << this->res_.reduction
408 <<
" (" << std::setprecision(4) << this->res_.elapsed <<
"s)" 411 if(!this->keep_matrix_)
415 template<
class GOS,
class TrlV,
class TstV>
418 typedef GOS GridOperator;
419 typedef TrlV TrialVector;
421 typedef typename TstV::ElementType RFType;
427 , force_iteration_(false)
429 this->reduction_ = 1
e-8;
430 this->abs_limit_ = 1
e-12;
436 , force_iteration_(false)
438 this->reduction_ = 1
e-8;
439 this->abs_limit_ = 1
e-12;
444 this->reduction_ = reduction;
454 force_iteration_ = force_iteration;
459 this->abs_limit_ = abs_limit_;
464 if (force_iteration_ && this->res_.iterations == 0)
466 this->res_.converged = this->res_.defect < this->abs_limit_
467 || this->res_.defect < this->res_.first_defect * this->reduction_;
468 if (this->res_.iterations >= maxit_ && !this->res_.converged)
470 "NewtonTerminate::terminate(): Maximum iteration count reached");
471 return this->res_.converged;
476 bool force_iteration_;
479 template<
class GOS,
class TrlV,
class TstV>
482 typedef GOS GridOperator;
483 typedef TrlV TrialVector;
485 typedef typename TstV::ElementType RFType;
486 typedef typename GOS::Traits::Jacobian
Matrix;
491 , min_linear_reduction_(1
e-3)
492 , fixed_linear_reduction_(0.0)
493 , reassemble_threshold_(0.0)
498 , min_linear_reduction_(1
e-3)
499 , fixed_linear_reduction_(0.0)
500 , reassemble_threshold_(0.0)
511 min_linear_reduction_ = min_linear_reduction;
520 fixed_linear_reduction_ = fixed_linear_reduction;
532 reassemble_threshold_ = reassemble_threshold;
537 this->reassembled_ =
false;
538 if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
540 if (this->verbosity_level_ >= 3)
541 std::cout <<
" Reassembling matrix..." << std::endl;
543 this->gridoperator_.jacobian(*this->u_, A);
544 this->reassembled_ =
true;
547 if (fixed_linear_reduction_ ==
true)
548 this->linear_reduction_ = min_linear_reduction_;
552 std::max(this->res_.first_defect * this->reduction_,
563 if ( stop_defect/(10*this->res_.defect) >
564 this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_) )
565 this->linear_reduction_ =
566 stop_defect/(10*this->res_.defect);
568 this->linear_reduction_ =
569 std::min(min_linear_reduction_,this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_));
572 this->prev_defect_ = this->res_.defect;
574 ios_base_all_saver restorer(std::cout);
576 if (this->verbosity_level_ >= 3)
577 std::cout <<
" requested linear reduction: " 578 << std::setw(12) << std::setprecision(4) << std::scientific
579 << this->linear_reduction_ << std::endl;
583 RFType min_linear_reduction_;
584 bool fixed_linear_reduction_;
585 RFType reassemble_threshold_;
588 template<
class GOS,
class TrlV,
class TstV>
591 typedef GOS GridOperator;
592 typedef TrlV TrialVector;
593 typedef TstV TestVector;
595 typedef typename TestVector::ElementType RFType;
606 hackbuschReuskenAcceptBest };
610 , strategy_(hackbuschReusken)
612 , damping_factor_(0.5)
617 , strategy_(hackbuschReusken)
619 , damping_factor_(0.5)
624 strategy_ = strategy;
629 strategy_ = strategyFromName(strategy);
639 damping_factor_ = damping_factor;
644 if (strategy_ == noLineSearch)
646 this->u_->axpy(-1.0, z);
651 if (this->verbosity_level_ >= 4)
652 std::cout <<
" Performing line search..." << std::endl;
654 RFType best_lambda = 0.0;
655 RFType best_defect = this->res_.defect;
656 TrialVector prev_u(*this->u_);
658 ios_base_all_saver restorer(std::cout);
662 if (this->verbosity_level_ >= 4)
663 std::cout <<
" trying line search damping factor: " 664 << std::setw(12) << std::setprecision(4) << std::scientific
668 this->u_->axpy(-lambda, z);
674 if (this->verbosity_level_ >= 4)
675 std::cout <<
" NaNs detected" << std::endl;
678 if (this->res_.defect <= (1.0 - lambda/4) * this->prev_defect_)
680 if (this->verbosity_level_ >= 4)
681 std::cout <<
" line search converged" << std::endl;
685 if (this->res_.defect < best_defect)
687 best_defect = this->res_.defect;
688 best_lambda = lambda;
693 if (this->verbosity_level_ >= 4)
694 std::cout <<
" max line search iterations exceeded" << std::endl;
697 case hackbuschReusken:
701 "NewtonLineSearch::line_search(): line search failed, " 702 "max iteration count reached, " 703 "defect did not improve enough");
704 case hackbuschReuskenAcceptBest:
705 if (best_lambda == 0.0)
710 "NewtonLineSearch::line_search(): line search failed, " 711 "max iteration count reached, " 712 "defect did not improve in any of the iterations");
714 if (best_lambda != lambda)
717 this->u_->axpy(-best_lambda, z);
727 lambda *= damping_factor_;
730 if (this->verbosity_level_ >= 4)
731 std::cout <<
" line search damping factor: " 732 << std::setw(12) << std::setprecision(4) << std::scientific
733 << lambda << std::endl;
739 if (s ==
"noLineSearch")
741 if (s ==
"hackbuschReusken")
742 return hackbuschReusken;
743 if (s ==
"hackbuschReuskenAcceptBest")
744 return hackbuschReuskenAcceptBest;
745 DUNE_THROW(
Exception,
"unknown line search strategy" << s);
751 RFType damping_factor_;
754 template<
class GOS,
class S,
class TrlV,
class TstV = TrlV>
760 typedef GOS GridOperator;
762 typedef TrlV TrialVector;
765 Newton(
const GridOperator& go, TrialVector& u_, Solver& solver_)
772 Newton(
const GridOperator& go, Solver& solver_)
803 typedef typename TstV::ElementType RFType;
804 if (param.hasKey(
"VerbosityLevel"))
805 this->setVerbosityLevel(
806 param.get<
unsigned int>(
"VerbosityLevel"));
807 if (param.hasKey(
"Reduction"))
809 param.get<RFType>(
"Reduction"));
810 if (param.hasKey(
"MaxIterations"))
811 this->setMaxIterations(
812 param.get<
unsigned int>(
"MaxIterations"));
813 if (param.hasKey(
"ForceIteration"))
814 this->setForceIteration(
815 param.get<
bool>(
"ForceIteration"));
816 if (param.hasKey(
"AbsoluteLimit"))
817 this->setAbsoluteLimit(
818 param.get<RFType>(
"AbsoluteLimit"));
819 if (param.hasKey(
"MinLinearReduction"))
820 this->setMinLinearReduction(
821 param.get<RFType>(
"MinLinearReduction"));
822 if (param.hasKey(
"FixedLinearReduction"))
823 this->setFixedLinearReduction(
824 param.get<
bool>(
"FixedLinearReduction"));
825 if (param.hasKey(
"ReassembleThreshold"))
826 this->setReassembleThreshold(
827 param.get<RFType>(
"ReassembleThreshold"));
828 if (param.hasKey(
"LineSearchStrategy"))
829 this->setLineSearchStrategy(
830 param.get<std::string>(
"LineSearchStrategy"));
831 if (param.hasKey(
"LineSearchMaxIterations"))
832 this->setLineSearchMaxIterations(
833 param.get<
unsigned int>(
"LineSearchMaxIterations"));
834 if (param.hasKey(
"LineSearchDampingFactor"))
835 this->setLineSearchDampingFactor(
836 param.get<RFType>(
"LineSearchDampingFactor"));
837 if (param.hasKey(
"KeepMatrix"))
839 param.get<
bool>(
"KeepMatrix"));
845 #endif // DUNE_PDELAB_NEWTON_NEWTON_HH void setForceIteration(bool force_iteration)
Definition: newton.hh:452
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: newton.hh:105
virtual bool terminate() override
Definition: newton.hh:462
const Entity & e
Definition: localfunctionspace.hh:120
virtual void prepare_step(Matrix &A, TstV &) override
Definition: newton.hh:535
NewtonResult()
Definition: newton.hh:76
NewtonLineSearch(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:608
NewtonPrepareStep(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:489
const GridOperator & gridoperator_
Definition: newton.hh:124
void setReassembleThreshold(RFType reassemble_threshold)
set a threshold, when the linear operator is reassembled
Definition: newton.hh:530
void setVerbosityLevel(unsigned int verbosity_level)
Definition: newton.hh:96
void setAbsoluteLimit(RFType abs_limit_)
Definition: newton.hh:457
Newton(const GridOperator &go, TrialVector &u_, Solver &solver_)
Definition: newton.hh:765
TrialVector * u_
Definition: newton.hh:125
NewtonSolver(const GridOperator &go, TrialVector &u_, Solver &solver)
Definition: newton.hh:180
NewtonTerminate(const GridOperator &go)
Definition: newton.hh:433
RFType linear_reduction_
Definition: newton.hh:132
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: newton.hh:111
unsigned int verbosity_level_
Definition: newton.hh:130
virtual void line_search(TrialVector &z, TestVector &r) override
Definition: newton.hh:642
NewtonPrepareStep(const GridOperator &go)
Definition: newton.hh:496
typename impl::BackendMatrixSelector< Backend, VU, VV, E >::Type Matrix
alias of the return type of BackendMatrixSelector
Definition: backend/interface.hh:127
void apply()
Definition: newton.hh:252
Definition: newton.hh:416
virtual ~NewtonBase()
Definition: newton.hh:158
RFType prev_defect_
Definition: newton.hh:131
RFType reduction_
Definition: newton.hh:134
void setParameters(const Dune::ParameterTree ¶m)
interpret a parameter tree as a set of options for the newton solver
Definition: newton.hh:801
don't do any line search or damping
Definition: newton.hh:600
NewtonTerminate(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:424
Strategy strategyFromName(const std::string &s)
Definition: newton.hh:738
Definition: newton.hh:589
NewtonSolver(const GridOperator &go, Solver &solver)
Definition: newton.hh:186
RFType abs_limit_
Definition: newton.hh:135
void setLineSearchMaxIterations(unsigned int maxit)
Definition: newton.hh:632
double assembler_time
Definition: newton.hh:72
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
void setLineSearchDampingFactor(RFType damping_factor)
Definition: newton.hh:637
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: newton.hh:117
NewtonResult< RFType > Result
Definition: newton.hh:94
bool keep_matrix_
Definition: newton.hh:136
Definition: newton.hh:167
void setMinLinearReduction(RFType min_linear_reduction)
set the minimal reduction in the linear solver
Definition: newton.hh:509
Definition: newton.hh:480
void setLineSearchStrategy(Strategy strategy)
Definition: newton.hh:622
const Result & result() const
Definition: newton.hh:196
NewtonLineSearch(const GridOperator &go)
Definition: newton.hh:615
double linear_solver_time
Definition: newton.hh:73
Newton(const GridOperator &go, Solver &solver_)
Definition: newton.hh:772
RFType defect
Definition: newton.hh:71
void setLineSearchStrategy(std::string strategy)
Definition: newton.hh:627
Definition: newton.hh:755
Result res_
Definition: newton.hh:129
NewtonBase(const GridOperator &go, TrialVector &u)
Definition: newton.hh:138
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
void setMaxIterations(unsigned int maxit)
Definition: newton.hh:447
int linear_solver_iterations
Definition: newton.hh:74
std::shared_ptr< Matrix > A_
Definition: newton.hh:128
Strategy
Definition: newton.hh:598
RFType first_defect
Definition: newton.hh:70
NewtonBase(const GridOperator &go)
Definition: newton.hh:148
std::shared_ptr< TestVector > r_
Definition: newton.hh:127
NewtonResult< RFType > Result
Definition: newton.hh:178
const std::string s
Definition: function.hh:830
std::shared_ptr< TrialVector > z_
Definition: newton.hh:126
void setReduction(RFType reduction)
Definition: newton.hh:442
perform a linear search for the optimal damping parameter with multiples of damping ...
Definition: newton.hh:604
bool reassembled_
Definition: newton.hh:133
virtual void defect(TestVector &r) override
Definition: newton.hh:205
void setFixedLinearReduction(bool fixed_linear_reduction)
set a fixed reduction in the linear solver (overwrites setMinLinearReduction)
Definition: newton.hh:518