3 #ifndef DUNE_PDELAB_NEWTON_HH
4 #define DUNE_PDELAB_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)
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;
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)
165 "NewtonSolver::result() called before NewtonSolver::solve()");
177 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
182 void linearSolve(Matrix& A, TrialVector& z, TestVector& r)
const
185 std::cout <<
" Solving linear system..." << std::endl;
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");
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)
399 , force_iteration_(false)
417 force_iteration_ = force_iteration;
433 "NewtonTerminate::terminate(): Maximum iteration count reached");
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;
501 if (this->
res_.
defect/this->prev_defect_ > reassemble_threshold_)
504 std::cout <<
" Reassembling matrix..." << std::endl;
510 if (fixed_linear_reduction_ ==
true)
527 this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_) )
532 std::min(min_linear_reduction_,this->
res_.
defect*this->res_.defect/(this->prev_defect_*this->prev_defect_));
537 ios_base_all_saver restorer(std::cout);
540 std::cout <<
" requested linear reduction: "
541 << std::setw(12) << std::setprecision(4) << std::scientific
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;
575 , damping_factor_(0.5)
582 , damping_factor_(0.5)
587 strategy_ = strategy;
602 damping_factor_ = damping_factor;
609 this->
u_->axpy(-1.0, z);
615 std::cout <<
" Performing line search..." << std::endl;
617 RFType best_lambda = 0.0;
619 TrialVector prev_u(*this->
u_);
621 ios_base_all_saver restorer(std::cout);
626 std::cout <<
" trying line search damping factor: "
627 << std::setw(12) << std::setprecision(4) << std::scientific
631 this->
u_->axpy(-lambda, z);
638 std::cout <<
" Nans detected" << std::endl;
641 if (this->
res_.
defect <= (1.0 - lambda/4) * this->prev_defect_)
644 std::cout <<
" line search converged" << std::endl;
651 best_lambda = lambda;
657 std::cout <<
" max line search iterations exceeded" << std::endl;
664 "NewtonLineSearch::line_search(): line search failed, "
665 "max iteration count reached, "
666 "defect did not improve enough");
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_;
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")
706 if (s ==
"hackbuschReuskenAcceptBest")
708 DUNE_THROW(
Exception,
"unknown linesearch 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"))
769 param.get<
unsigned int>(
"VerbosityLevel"));
770 if (param.hasKey(
"Reduction"))
772 param.get<RFType>(
"Reduction"));
773 if (param.hasKey(
"MaxIterations"))
775 param.get<
unsigned int>(
"MaxIterations"));
776 if (param.hasKey(
"ForceIteration"))
778 param.get<
bool>(
"ForceIteration"));
779 if (param.hasKey(
"AbsoluteLimit"))
781 param.get<RFType>(
"AbsoluteLimit"));
782 if (param.hasKey(
"MinLinearReduction"))
784 param.get<RFType>(
"MinLinearReduction"));
785 if (param.hasKey(
"FixedLinearReduction"))
787 param.get<
bool>(
"FixedLinearReduction"));
788 if (param.hasKey(
"ReassembleThreshold"))
790 param.get<RFType>(
"ReassembleThreshold"));
791 if (param.hasKey(
"LineSearchStrategy"))
793 param.get<std::string>(
"LineSearchStrategy"));
794 if (param.hasKey(
"LineSearchMaxIterations"))
796 param.get<
unsigned int>(
"LineSearchMaxIterations"));
797 if (param.hasKey(
"LineSearchDampingFactor"))
799 param.get<RFType>(
"LineSearchDampingFactor"));
800 if (param.hasKey(
"KeepMatrix"))
802 param.get<
bool>(
"KeepMatrix"));
808 #endif // DUNE_PDELAB_NEWTON_HH
Definition: newton.hh:552
void setParameters(Dune::ParameterTree ¶m)
interpret a parameter tree as a set of options for the newton solver
Definition: newton.hh:764
double linear_solver_time
Definition: newton.hh:38
NewtonSolver(const GridOperator &go, Solver &solver)
Definition: newton.hh:151
NewtonSolver(const GridOperator &go, TrialVector &u_, Solver &solver)
Definition: newton.hh:145
void setFixedLinearReduction(bool fixed_linear_reduction)
set a fixed reduction in the linear solver (overwrites setMinLinearReduction)
Definition: newton.hh:481
virtual bool terminate()
Definition: newton.hh:425
Definition: newton.hh:443
Strategy
Definition: newton.hh:561
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: newton.hh:70
NewtonBase(const GridOperator &go)
Definition: newton.hh:113
bool reassembled_
Definition: newton.hh:98
virtual void prepare_step(Matrix &A, TestVector &r)=0
Definition: newton.hh:379
unsigned int verbosity_level_
Definition: newton.hh:95
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: newton.hh:76
void setLineSearchMaxIterations(unsigned int maxit)
Definition: newton.hh:595
void setAbsoluteLimit(RFType abs_limit_)
Definition: newton.hh:420
virtual bool terminate()=0
Newton(const GridOperator &go, TrialVector &u_, Solver &solver_)
Definition: newton.hh:728
NewtonLineSearch(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:571
const E & e
Definition: interpolate.hh:172
TrialVector * u_
Definition: newton.hh:90
RFType abs_limit_
Definition: newton.hh:100
NewtonPrepareStep(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:452
RFType prev_defect_
Definition: newton.hh:96
NewtonResult< RFType > Result
Definition: newton.hh:143
void setVerbosityLevel(unsigned int verbosity_level)
Definition: newton.hh:61
void setMinLinearReduction(RFType min_linear_reduction)
set the minimal reduction in the linear solver
Definition: newton.hh:472
virtual void line_search(TrialVector &z, TestVector &r)
Definition: newton.hh:605
virtual void defect(TestVector &r)
Definition: newton.hh:170
virtual void line_search(TrialVector &z, TestVector &r)=0
unsigned int iterations
Definition: solver.hh:33
void apply()
Definition: newton.hh:215
NewtonTerminate(const GridOperator &go)
Definition: newton.hh:396
virtual ~NewtonBase()
Definition: newton.hh:123
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: newton.hh:82
Newton(const GridOperator &go, Solver &solver_)
Definition: newton.hh:735
const Result & result() const
Definition: newton.hh:161
bool converged
Definition: solver.hh:32
Strategy strategyFromName(const std::string &s)
Definition: newton.hh:701
NewtonResult()
Definition: newton.hh:41
Definition: newton.hh:132
void setForceIteration(bool force_iteration)
Definition: newton.hh:415
Result res_
Definition: newton.hh:94
const GridOperator & gridoperator_
Definition: newton.hh:89
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
std::shared_ptr< Matrix > A_
Definition: newton.hh:93
Definition: adaptivity.hh:27
void setMaxIterations(unsigned int maxit)
Definition: newton.hh:410
same as hackbuschReusken, but doesn't fail if the best update is still not good enough ...
Definition: newton.hh:569
NewtonResult< RFType > Result
Definition: newton.hh:59
bool keep_matrix_
Definition: newton.hh:101
virtual void defect(TestVector &r)=0
perform a linear search for the optimal damping parameter with multiples of damping ...
Definition: newton.hh:567
RFType reduction_
Definition: newton.hh:99
NewtonLineSearch(const GridOperator &go)
Definition: newton.hh:578
void setReassembleThreshold(RFType reassemble_threshold)
set a threshold, when the linear operator is reassembled
Definition: newton.hh:493
const std::string s
Definition: function.hh:1102
void setLineSearchStrategy(Strategy strategy)
Definition: newton.hh:585
void setLineSearchStrategy(std::string strategy)
Definition: newton.hh:590
std::shared_ptr< TrialVector > z_
Definition: newton.hh:91
NewtonPrepareStep(const GridOperator &go)
Definition: newton.hh:459
NewtonBase(const GridOperator &go, TrialVector &u)
Definition: newton.hh:103
NewtonTerminate(const GridOperator &go, TrialVector &u_)
Definition: newton.hh:387
RFType linear_reduction_
Definition: newton.hh:97
RFType defect
Definition: newton.hh:36
int linear_solver_iterations
Definition: newton.hh:39
virtual void prepare_step(Matrix &A, TstV &)
Definition: newton.hh:498
void setReduction(RFType reduction)
Definition: newton.hh:405
double assembler_time
Definition: newton.hh:37
RFType first_defect
Definition: newton.hh:35
don't do any linesearch or damping
Definition: newton.hh:563
std::shared_ptr< TestVector > r_
Definition: newton.hh:92
void setLineSearchDampingFactor(RFType damping_factor)
Definition: newton.hh:600
Definition: newton.hh:718