3 #ifndef DUNE_PDELAB_NEWTON_HH
4 #define DUNE_PDELAB_NEWTON_HH
12 #include <dune/common/exceptions.hh>
13 #include <dune/common/ios_state.hh>
14 #include <dune/common/timer.hh>
15 #include <dune/common/parametertree.hh>
31 template<
class RFType>
45 template<
class GOS,
class TrlV,
class TstV>
48 typedef GOS GridOperator;
49 typedef TrlV TrialVector;
50 typedef TstV TestVector;
52 typedef typename TestVector::ElementType RFType;
53 typedef typename GOS::Traits::Jacobian Matrix;
62 if (
gridoperator.trialGridFunctionSpace().gridView().comm().rank()>0)
84 if (
gridoperator.trialGridFunctionSpace().gridView().comm().rank()>0)
93 if (
gridoperator.trialGridFunctionSpace().gridView().comm().rank()>0)
100 virtual void prepare_step(Matrix& A, TestVector& r) = 0;
101 virtual void line_search(TrialVector& z, TestVector& r) = 0;
102 virtual void defect(TestVector& r) = 0;
105 template<
class GOS,
class S,
class TrlV,
class TstV>
109 typedef GOS GridOperator;
110 typedef TrlV TrialVector;
111 typedef TstV TestVector;
113 typedef typename TestVector::ElementType RFType;
114 typedef typename GOS::Traits::Jacobian Matrix;
122 , result_valid(false)
128 , result_valid(false)
133 void apply(TrialVector& u_);
139 "NewtonSolver::result() called before NewtonSolver::solve()");
151 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
156 void linearSolve(Matrix& A, TrialVector& z, TestVector& r)
const
159 std::cout <<
" Solving linear system..." << std::endl;
163 ios_base_all_saver restorer(std::cout);
165 if (!this->solver.result().converged)
167 "NewtonSolver::linearSolve(): Linear solver did not converge "
168 "in " << this->solver.result().iterations <<
" iterations");
170 std::cout <<
" linear solver iterations: "
171 << std::setw(12) << solver.result().iterations << std::endl
172 <<
" linear defect reduction: "
173 << std::setw(12) << std::setprecision(4) << std::scientific
174 << solver.result().reduction << std::endl;
181 template<
class GOS,
class S,
class TrlV,
class TstV>
188 template<
class GOS,
class S,
class TrlV,
class TstV>
191 this->res.iterations = 0;
192 this->res.converged =
false;
193 this->res.reduction = 1.0;
194 this->res.conv_rate = 1.0;
195 this->res.elapsed = 0.0;
196 this->res.assembler_time = 0.0;
197 this->res.linear_solver_time = 0.0;
198 this->res.linear_solver_iterations = 0;
204 TestVector r(this->gridoperator.testGridFunctionSpace());
206 this->res.first_defect = this->res.defect;
207 this->prev_defect = this->res.defect;
209 if (this->verbosity_level >= 2)
212 ios_base_all_saver restorer(std::cout);
213 std::cout <<
" Initial defect: "
214 << std::setw(12) << std::setprecision(4) << std::scientific
215 << this->res.defect << std::endl;
218 Matrix A(this->gridoperator);
219 TrialVector z(this->gridoperator.trialGridFunctionSpace());
221 while (!this->terminate())
223 if (this->verbosity_level >= 3)
224 std::cout <<
" Newton iteration " << this->res.iterations
225 <<
" --------------------------------" << std::endl;
227 Timer assembler_timer;
230 this->prepare_step(A,r);
234 this->res.assembler_time += assembler_timer.elapsed();
237 double assembler_time = assembler_timer.elapsed();
238 this->res.assembler_time += assembler_time;
239 if (this->verbosity_level >= 3)
240 std::cout <<
" matrix assembly time: "
241 << std::setw(12) << std::setprecision(4) << std::scientific
242 << assembler_time << std::endl;
244 Timer linear_solver_timer;
247 this->linearSolve(A, z, r);
251 this->res.linear_solver_time += linear_solver_timer.elapsed();
252 this->res.linear_solver_iterations += this->solver.result().iterations;
255 double linear_solver_time = linear_solver_timer.elapsed();
256 this->res.linear_solver_time += linear_solver_time;
257 this->res.linear_solver_iterations += this->solver.result().iterations;
261 this->line_search(z, r);
265 if (this->reassembled)
267 if (this->verbosity_level >= 3)
268 std::cout <<
" line search failed - trying again with reassembled matrix" << std::endl;
272 this->res.reduction = this->res.defect/this->res.first_defect;
273 this->res.iterations++;
274 this->res.conv_rate = std::pow(this->res.reduction, 1.0/this->res.iterations);
277 ios_base_all_saver restorer(std::cout);
279 if (this->verbosity_level >= 3)
280 std::cout <<
" linear solver time: "
281 << std::setw(12) << std::setprecision(4) << std::scientific
282 << linear_solver_time << std::endl
283 <<
" defect reduction (this iteration):"
284 << std::setw(12) << std::setprecision(4) << std::scientific
285 << this->res.defect/this->prev_defect << std::endl
286 <<
" defect reduction (total): "
287 << std::setw(12) << std::setprecision(4) << std::scientific
288 << this->res.reduction << std::endl
290 << std::setw(12) << std::setprecision(4) << std::scientific
291 << this->res.defect << std::endl;
292 if (this->verbosity_level == 2)
293 std::cout <<
" Newton iteration " << std::setw(2) << this->res.iterations
295 << std::setw(12) << std::setprecision(4) << std::scientific
297 <<
". Reduction (this): "
298 << std::setw(12) << std::setprecision(4) << std::scientific
299 << this->res.defect/this->prev_defect
300 <<
". Reduction (total): "
301 << std::setw(12) << std::setprecision(4) << std::scientific
302 << this->res.reduction << std::endl;
307 this->res.elapsed = timer.elapsed();
310 this->res.elapsed = timer.elapsed();
312 ios_base_all_saver restorer(std::cout);
314 if (this->verbosity_level == 1)
315 std::cout <<
" Newton converged after " << std::setw(2) << this->res.iterations
316 <<
" iterations. Reduction: "
317 << std::setw(12) << std::setprecision(4) << std::scientific
318 << this->res.reduction
319 <<
" (" << std::setprecision(4) << this->res.elapsed <<
"s)"
323 template<
class GOS,
class TrlV,
class TstV>
326 typedef GOS GridOperator;
327 typedef TrlV TrialVector;
329 typedef typename TstV::ElementType RFType;
335 , force_iteration(false)
344 , force_iteration(false)
362 force_iteration = force_iteration_;
378 "NewtonTerminate::terminate(): Maximum iteration count reached");
384 bool force_iteration;
387 template<
class GOS,
class TrlV,
class TstV>
390 typedef GOS GridOperator;
391 typedef TrlV TrialVector;
393 typedef typename TstV::ElementType RFType;
394 typedef typename GOS::Traits::Jacobian Matrix;
399 , min_linear_reduction(1
e-3)
400 , fixed_linear_reduction(0.0)
401 , reassemble_threshold(0.0)
406 , min_linear_reduction(1
e-3)
407 , fixed_linear_reduction(0.0)
408 , reassemble_threshold(0.0)
417 min_linear_reduction = min_linear_reduction_;
424 fixed_linear_reduction = fixed_linear_reduction_;
429 reassemble_threshold = reassemble_threshold_;
435 if (this->
res.
defect/this->prev_defect > reassemble_threshold)
438 std::cout <<
" Reassembling matrix..." << std::endl;
444 if (fixed_linear_reduction ==
true)
461 this->res.defect*this->res.defect/(this->prev_defect*this->prev_defect) )
466 std::min(min_linear_reduction,this->
res.
defect*this->res.defect/(this->prev_defect*this->prev_defect));
471 ios_base_all_saver restorer(std::cout);
474 std::cout <<
" requested linear reduction: "
475 << std::setw(12) << std::setprecision(4) << std::scientific
480 RFType min_linear_reduction;
481 bool fixed_linear_reduction;
482 RFType reassemble_threshold;
485 template<
class GOS,
class TrlV,
class TstV>
488 typedef GOS GridOperator;
489 typedef TrlV TrialVector;
490 typedef TstV TestVector;
492 typedef typename TestVector::ElementType RFType;
503 , damping_factor(0.5)
510 , damping_factor(0.5)
515 strategy = strategy_;
530 damping_factor = damping_factor_;
537 this->
u->axpy(-1.0, z);
543 std::cout <<
" Performing line search..." << std::endl;
545 RFType best_lambda = 0.0;
547 TrialVector prev_u(*this->
u);
549 ios_base_all_saver restorer(std::cout);
554 std::cout <<
" trying line search damping factor: "
555 << std::setw(12) << std::setprecision(4) << std::scientific
559 this->
u->axpy(-lambda, z);
566 std::cout <<
" Nans detected" << std::endl;
569 if (this->
res.
defect <= (1.0 - lambda/4) * this->prev_defect)
572 std::cout <<
" line search converged" << std::endl;
579 best_lambda = lambda;
585 std::cout <<
" max line search iterations exceeded" << std::endl;
592 "NewtonLineSearch::line_search(): line search failed, "
593 "max iteration count reached, "
594 "defect did not improve enough");
596 if (best_lambda == 0.0)
601 "NewtonLineSearch::line_search(): line search failed, "
602 "max iteration count reached, "
603 "defect did not improve in any of the iterations");
605 if (best_lambda != lambda)
608 this->
u->axpy(-best_lambda, z);
618 lambda *= damping_factor;
622 std::cout <<
" line search damping factor: "
623 << std::setw(12) << std::setprecision(4) << std::scientific
624 << lambda << std::endl;
630 if (s ==
"noLineSearch")
632 if (s ==
"hackbuschReusken")
634 if (s ==
"hackbuschReuskenAcceptBest")
641 RFType damping_factor;
644 template<
class GOS,
class S,
class TrlV,
class TstV = TrlV>
650 typedef GOS GridOperator;
652 typedef TrlV TrialVector;
655 Newton(GridOperator& go, TrialVector& u_, Solver& solver_)
662 Newton(GridOperator& go, Solver& solver_)
693 typedef typename TstV::ElementType RFType;
694 if (param.hasKey(
"VerbosityLevel"))
696 param.get<
unsigned int>(
"VerbosityLevel"));
697 if (param.hasKey(
"VerbosityLevel"))
699 param.get<
unsigned int>(
"VerbosityLevel"));
700 if (param.hasKey(
"Reduction"))
702 param.get<RFType>(
"Reduction"));
703 if (param.hasKey(
"MaxIterations"))
705 param.get<
unsigned int>(
"MaxIterations"));
706 if (param.hasKey(
"ForceIteration"))
708 param.get<
bool>(
"ForceIteration"));
709 if (param.hasKey(
"AbsoluteLimit"))
711 param.get<RFType>(
"AbsoluteLimit"));
712 if (param.hasKey(
"MinLinearReduction"))
714 param.get<RFType>(
"MinLinearReduction"));
715 if (param.hasKey(
"FixedLinearReduction"))
717 param.get<
bool>(
"FixedLinearReduction"));
718 if (param.hasKey(
"ReassembleThreshold"))
720 param.get<RFType>(
"ReassembleThreshold"));
721 if (param.hasKey(
"LineSearchStrategy"))
723 param.get<std::string>(
"LineSearchStrategy"));
724 if (param.hasKey(
"LineSearchMaxIterations"))
726 param.get<
unsigned int>(
"LineSearchMaxIterations"));
727 if (param.hasKey(
"LineSearchDampingFactor"))
729 param.get<RFType>(
"LineSearchDampingFactor"));
735 #endif // DUNE_PDELAB_NEWTON_HH
Definition: newton.hh:106
RFType defect
Definition: newton.hh:35
void setMinLinearReduction(RFType min_linear_reduction_)
Definition: newton.hh:415
void setAbsoluteLimit(RFType abs_limit_)
Definition: newton.hh:365
virtual void prepare_step(Matrix &A, TestVector &r)=0
Strategy strategyFromName(const std::string &s)
Definition: newton.hh:629
Definition: newton.hh:497
NewtonSolver(GridOperator &go, Solver &solver_)
Definition: newton.hh:125
bool converged
Definition: solver.hh:32
RFType first_defect
Definition: newton.hh:34
NewtonResult()
Definition: newton.hh:40
void setLineSearchMaxIterations(unsigned int maxit_)
Definition: newton.hh:523
Definition: newton.hh:496
virtual ~NewtonBase()
Definition: newton.hh:97
virtual void line_search(TrialVector &z, TestVector &r)=0
void setLineSearchStrategy(std::string strategy_)
Definition: newton.hh:518
Newton(GridOperator &go, Solver &solver_)
Definition: newton.hh:662
RFType reduction
Definition: newton.hh:76
void setReassembleThreshold(RFType reassemble_threshold_)
Definition: newton.hh:427
void setLineSearchStrategy(Strategy strategy_)
Definition: newton.hh:513
virtual bool terminate()=0
NewtonTerminate(GridOperator &go, TrialVector &u_)
Definition: newton.hh:332
NewtonLineSearch(GridOperator &go, TrialVector &u_)
Definition: newton.hh:499
NewtonBase(GridOperator &go)
Definition: newton.hh:88
TrialVector * u
Definition: newton.hh:70
RFType abs_limit
Definition: newton.hh:77
virtual void defect(TestVector &r)
Definition: newton.hh:144
NewtonResult< RFType > Result
Definition: newton.hh:58
void setVerbosityLevel(unsigned int verbosity_level_)
Definition: newton.hh:60
bool reassembled
Definition: newton.hh:75
Newton(GridOperator &go, TrialVector &u_, Solver &solver_)
Definition: newton.hh:655
NewtonPrepareStep(GridOperator &go)
Definition: newton.hh:404
virtual void prepare_step(Matrix &A, TstV &)
Definition: newton.hh:432
Definition: newton.hh:645
void setReduction(RFType reduction_)
Definition: newton.hh:350
Definition: newton.hh:495
virtual void defect(TestVector &r)=0
NewtonResult< RFType > Result
Definition: newton.hh:117
NewtonLineSearch(GridOperator &go)
Definition: newton.hh:506
double assembler_time
Definition: newton.hh:36
Definition: newton.hh:388
int linear_solver_iterations
Definition: newton.hh:38
double linear_solver_time
Definition: newton.hh:37
GridOperator & gridoperator
Definition: newton.hh:69
Definition: newton.hh:324
Result res
Definition: newton.hh:71
Strategy
Definition: newton.hh:495
NewtonPrepareStep(GridOperator &go, TrialVector &u_)
Definition: newton.hh:397
NewtonTerminate(GridOperator &go)
Definition: newton.hh:341
virtual bool terminate()
Definition: newton.hh:370
void setParameters(Dune::ParameterTree ¶m)
interpret a parameter tree as a set of options for the newton solver
Definition: newton.hh:691
Definition: newton.hh:486
NewtonBase(GridOperator &go, TrialVector &u_)
Definition: newton.hh:79
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
void setLineSearchDampingFactor(RFType damping_factor_)
Definition: newton.hh:528
const E & e
Definition: interpolate.hh:172
void setForceIteration(bool force_iteration_)
Definition: newton.hh:360
void setFixedLinearReduction(bool fixed_linear_reduction_)
Definition: newton.hh:422
virtual void line_search(TrialVector &z, TestVector &r)
Definition: newton.hh:533
NewtonSolver(GridOperator &go, TrialVector &u_, Solver &solver_)
Definition: newton.hh:119
unsigned int verbosity_level
Definition: newton.hh:72
void apply()
Definition: newton.hh:189
void setMaxIterations(unsigned int maxit_)
Definition: newton.hh:355
RFType prev_defect
Definition: newton.hh:73
const std::string s
Definition: function.hh:1103
const Result & result() const
Definition: newton.hh:135
RFType linear_reduction
Definition: newton.hh:74
unsigned int iterations
Definition: solver.hh:33