3#ifndef DUNE_PDELAB_SOLVER_NEWTON_HH
4#define DUNE_PDELAB_SOLVER_NEWTON_HH
6#include <dune/common/exceptions.hh>
7#include <dune/common/ios_state.hh>
19 template<
typename T1,
typename =
void>
25 struct HasSetReuse<T, decltype(
std::declval<T>().setReuse(true), void())>
30 inline void setLinearSystemReuse(T& solver_backend,
bool reuse, std::true_type)
32 if (!solver_backend.getReuse() && reuse)
33 dwarn <<
"WARNING: Newton needed to override your choice to reuse the linear system in order to work!" << std::endl;
34 solver_backend.setReuse(reuse);
38 inline void setLinearSystemReuse(T&,
bool, std::false_type)
42 inline void setLinearSystemReuse(T& solver_backend,
bool reuse)
44 setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
61 template <
typename Gr
idOperator_,
typename LinearSolver_>
81 using Real =
typename Dune::FieldTraits<typename Domain::ElementType>::real_type;
90 DUNE_THROW(NewtonError,
"NewtonMethod::result() called before NewtonMethod::solve()");
97 if (_result.
defect/_previousDefect > _reassembleThreshold){
98 if (_hangingNodeModifications){
99 auto dirichletValues = solution;
107 _gridOperator.localAssembler().backtransform(solution);
110 std::cout <<
" Reassembling matrix..." << std::endl;
111 *_jacobian =
Real(0.0);
112 _gridOperator.jacobian(solution, *_jacobian);
116 _linearReduction = _minLinearReduction;
117 if (not _fixedLinearReduction){
121 auto stop_defect = max(_result.
first_defect*_reduction, _absoluteLimit);
127 if (stop_defect/(10*_result.
defect) > _result.
defect*_result.
defect/(_previousDefect*_previousDefect))
128 _linearReduction = stop_defect/(10*_result.
defect);
130 _linearReduction = min(_minLinearReduction, _result.
defect*_result.
defect/(_previousDefect*_previousDefect));
134 std::cout <<
" requested linear reduction: "
135 << std::setw(12) << std::setprecision(4) << std::scientific
136 << _linearReduction << std::endl;
142 std::cout <<
" Solving linear system..." << std::endl;
145 Impl::setLinearSystemReuse(_linearSolver, not _reassembled);
149 _linearSolver.apply(*_jacobian, _correction, _residual, _linearReduction);
151 if (not _linearSolver.result().converged)
152 DUNE_THROW(NewtonLinearSolverError,
153 "NewtonMethod::linearSolve(): Linear solver did not converge "
154 "in " << _linearSolver.result().iterations <<
" iterations");
156 std::cout <<
" linear solver iterations: "
157 << std::setw(12) << _linearSolver.result().iterations << std::endl
158 <<
" linear defect reduction: "
159 << std::setw(12) << std::setprecision(4) << std::scientific
160 << _linearSolver.result().reduction << std::endl;
171 ios_base_all_saver restorer(std::cout);
174 using Clock = std::chrono::steady_clock;
175 using Duration = Clock::duration;
176 auto assembler_time = Duration::zero();
177 auto linear_solver_time = Duration::zero();
178 auto line_search_time = Duration::zero();
180 [](Duration duration){
181 return std::chrono::duration<double>(duration).count();
183 auto start_solve = Clock::now();
190 _previousDefect = _result.
defect;
193 std::cout <<
" Initial defect: "
194 << std::setw(12) << std::setprecision(4) << std::scientific
195 << _result.
defect << std::endl;
201 _jacobian = std::make_shared<Jacobian>(_gridOperator);
206 while (not _terminate->terminate()){
208 std::cout <<
" Newton iteration " << _result.
iterations
209 <<
" --------------------------------" << std::endl;
214 auto start = Clock::now();
216 auto end = Clock::now();
217 assembler_time += end -start;
220 _previousDefect = _result.
defect;
225 start = Clock::now();
228 linear_solver_time += end -start;
233 start = Clock::now();
234 _lineSearch->lineSearch(solution, _correction);
237 line_search_time += end -start;
247 std::cout <<
" linear solver time: "
248 << std::setw(12) << std::setprecision(4) << std::scientific
249 << to_seconds(linear_solver_time) << std::endl
250 <<
" defect reduction (this iteration):"
251 << std::setw(12) << std::setprecision(4) << std::scientific
252 << _result.
defect/_previousDefect << std::endl
253 <<
" defect reduction (total): "
254 << std::setw(12) << std::setprecision(4) << std::scientific
257 << std::setw(12) << std::setprecision(4) << std::scientific
258 << _result.
defect << std::endl;
260 std::cout <<
" Newton iteration "
264 << std::setw(12) << std::setprecision(4) << std::scientific
266 <<
". Reduction (this): "
267 << std::setw(12) << std::setprecision(4) << std::scientific
268 << _result.
defect/_previousDefect
269 <<
". Reduction (total): "
270 << std::setw(12) << std::setprecision(4) << std::scientific
275 auto end_solve = Clock::now();
276 auto solve_time = end_solve - start_solve;
277 _result.
elapsed = to_seconds(solve_time);
280 std::cout <<
" Newton converged after "
283 <<
" iterations. Reduction: "
284 << std::setw(12) << std::setprecision(4) << std::scientific
286 <<
" (" << std::setprecision(4) << _result.
elapsed <<
"s)"
296 if (_hangingNodeModifications){
297 auto dirichletValues = solution;
305 _gridOperator.localAssembler().backtransform(solution);
309 _gridOperator.residual(solution, _residual);
316 _result.
defect = _gridOperator.testGridFunctionSpace().gridView().comm().max(rankMax);
319 _result.
defect = _linearSolver.norm(_residual);
325 if (_gridOperator.trialGridFunctionSpace().gridView().comm().rank()>0)
328 _verbosity = verbosity;
340 _reduction = reduction;
352 _absoluteLimit = absoluteLimit;
357 return _absoluteLimit;
375 _hangingNodeModifications = b;
401 _minLinearReduction = minLinearReduction;
411 _fixedLinearReduction = fixedLinearReduction;
422 _reassembleThreshold = reassembleThreshold;
459 if (parameterTree.hasKey(
"verbosity"))
461 if (parameterTree.hasKey(
"reduction"))
463 if (parameterTree.hasKey(
"absolute_limit"))
465 if (parameterTree.hasKey(
"keeep_matrix"))
467 if (parameterTree.hasKey(
"use_max_norm"))
469 if (parameterTree.hasKey(
"hanging_node_modifications"))
470 _hangingNodeModifications = parameterTree.get<
bool>(
"hanging_node_modifications");
472 if (parameterTree.hasKey(
"min_linear_reduction"))
474 if (parameterTree.hasKey(
"fixed_linear_reduction"))
476 if (parameterTree.hasKey(
"reassemble_threshold"))
479 _terminate->setParameters(parameterTree.sub(
"terminate"));
480 _lineSearch->setParameters(parameterTree.sub(
"line_search"));
486 _terminate = terminate;
495 _lineSearch = lineSearch;
502 const std::string& lineSearchStrategy=
"hackbusch_reusken")
503 : _gridOperator(gridOperator)
504 , _linearSolver(linearSolver)
505 , _residual(gridOperator.testGridFunctionSpace())
506 , _correction(gridOperator.trialGridFunctionSpace())
508 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
516 const ParameterTree& parameterTree,
517 const std::string& lineSearchStrategy=
"hackbusch_reusken")
518 : _gridOperator(gridOperator)
519 , _linearSolver(linearSolver)
520 , _residual(gridOperator.testGridFunctionSpace())
521 , _correction(gridOperator.trialGridFunctionSpace())
525 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
536 std::shared_ptr<Jacobian> _jacobian;
537 std::shared_ptr<Domain> _previousSolution;
539 std::shared_ptr<TerminateInterface> _terminate;
540 std::shared_ptr<LineSearchInterface<Domain>> _lineSearch;
544 bool _resultValid =
false;
545 Real _previousDefect = 0.0;
548 bool _reassembled =
true;
549 Real _linearReduction = 0.0;
552 unsigned int _verbosity = 0;
553 Real _reduction = 1e-8;
554 Real _absoluteLimit = 1e-12;
555 bool _keepMatrix =
true;
556 bool _useMaxNorm =
false;
559 bool _hangingNodeModifications =
false;
562 Real _minLinearReduction = 1e-3;
563 bool _fixedLinearReduction =
false;
564 Real _reassembleThreshold = 0.0;
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
Definition adaptivity.hh:29
std::shared_ptr< LineSearchInterface< typename Newton::Domain > > getLineSearch(Newton &newton, const std::string &name)
Get a pointer to a line search.
Definition newtonlinesearch.hh:225
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition backend/interface.hh:192
RFType conv_rate
Definition solver.hh:36
RFType reduction
Definition solver.hh:35
unsigned int iterations
Definition solver.hh:33
double elapsed
Definition solver.hh:34
Dune::PDELab::Backend::Matrix< MB, Domain, Range, JF > Jacobian
The type of the jacobian.
Definition gridoperatorutilities.hh:72
Dune::PDELab::Backend::Vector< GFSV, RF > Range
The type of the range (residual).
Definition gridoperatorutilities.hh:65
Dune::PDELab::Backend::Vector< GFSU, DF > Domain
The type of the domain (solution).
Definition gridoperatorutilities.hh:58
Standard grid operator implementation.
Definition gridoperator.hh:36
Newton solver for solving non-linear problems.
Definition solver/newton.hh:63
typename GridOperator::Traits::Jacobian Jacobian
Type of the Jacobian matrix.
Definition solver/newton.hh:78
const Result & result() const
Return results.
Definition solver/newton.hh:87
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver, const ParameterTree ¶meterTree, const std::string &lineSearchStrategy="hackbusch_reusken")
Construct Newton passing a parameter tree.
Definition solver/newton.hh:513
void setMinLinearReduction(Real minLinearReduction)
Set the minimal reduction in the linear solver.
Definition solver/newton.hh:399
void setTerminate(std::shared_ptr< TerminateInterface > terminate)
Set the termination criterion.
Definition solver/newton.hh:484
virtual void linearSolve()
Definition solver/newton.hh:139
Real getAbsoluteLimit() const
Definition solver/newton.hh:355
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition solver/newton.hh:380
unsigned int getVerbosityLevel() const
Get verbosity level.
Definition solver/newton.hh:332
typename GridOperator::Traits::Domain Domain
Type of the domain (solution)
Definition solver/newton.hh:72
void setLineSearch(std::shared_ptr< LineSearchInterface< Domain > > lineSearch)
Set the line search.
Definition solver/newton.hh:493
void setParameters(const ParameterTree ¶meterTree)
Interpret a parameter tree as a set of options for the newton solver.
Definition solver/newton.hh:458
virtual void prepareStep(Domain &solution)
Definition solver/newton.hh:94
void setHangingNodeModifications(bool b)
Does the problem have hanging nodes.
Definition solver/newton.hh:373
void setFixedLinearReduction(bool fixedLinearReduction)
Set wether to use a fixed reduction in the linear solver.
Definition solver/newton.hh:409
void discardMatrix()
Discard the stored Jacobian matrix.
Definition solver/newton.hh:386
GridOperator_ GridOperator
Type of the grid operator.
Definition solver/newton.hh:66
void setReduction(Real reduction)
Set reduction Newton needs to achieve.
Definition solver/newton.hh:338
typename Dune::FieldTraits< typename Domain::ElementType >::real_type Real
Number type.
Definition solver/newton.hh:81
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition solver/newton.hh:361
LinearSolver_ LinearSolver
Type of the linear solver.
Definition solver/newton.hh:69
void setUseMaxNorm(bool b)
Set whether to use the maximum norm for stopping criteria.
Definition solver/newton.hh:367
typename GridOperator::Traits::Range Range
Type of the range (residual)
Definition solver/newton.hh:75
void setAbsoluteLimit(Real absoluteLimit)
Set absolute convergence limit.
Definition solver/newton.hh:350
virtual void updateDefect(Domain &solution)
Update _residual and defect in _result.
Definition solver/newton.hh:294
virtual void apply(Domain &solution)
Solve the nonlinear problem using solution as initial guess and for storing the result.
Definition solver/newton.hh:164
PDESolverResult< Real > Result
Type of results.
Definition solver/newton.hh:84
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver, const std::string &lineSearchStrategy="hackbusch_reusken")
Construct Newton using default parameters.
Definition solver/newton.hh:499
void setReassembleThreshold(Real reassembleThreshold)
Set a threshold, when the linear operator is reassembled.
Definition solver/newton.hh:420
void setVerbosityLevel(unsigned int verbosity)
Set how much output you get.
Definition solver/newton.hh:323
Real getReduction() const
Get reduction.
Definition solver/newton.hh:344
Abstract base class describing the line search interface.
Definition newtonlinesearch.hh:13
RFType defect
Definition solver/utility.hh:11
RFType first_defect
Definition solver/utility.hh:10
void clear()
Definition solver/utility.hh:21