3#ifndef DUNE_PDELAB_SOLVER_NEWTON_HH
4#define DUNE_PDELAB_SOLVER_NEWTON_HH
9#include <dune/pdelab/constraints/common/constraints.hh>
10#include <dune/pdelab/solver/newtonerrors.hh>
11#include <dune/pdelab/solver/linesearch.hh>
12#include <dune/pdelab/solver/terminate.hh>
13#include <dune/pdelab/solver/utility.hh>
20 template<
typename T1,
typename =
void>
26 struct HasSetReuse<T, decltype(
std::declval<T>().setReuse(true), void())>
31 inline void setLinearSystemReuse(T& solver_backend,
bool reuse, std::true_type)
33 if (!solver_backend.getReuse() && reuse)
34 dwarn <<
"WARNING: Newton needed to override your choice to reuse the linear system in order to work!" << std::endl;
35 solver_backend.setReuse(reuse);
39 inline void setLinearSystemReuse(T&,
bool, std::false_type)
43 inline void setLinearSystemReuse(T& solver_backend,
bool reuse)
45 setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
62 template <
typename Gr
idOperator_,
typename LinearSolver_>
82 using Real =
typename Dune::FieldTraits<typename Domain::ElementType>::real_type;
85 using Result = PDESolverResult<Real>;
94 DUNE_THROW(NewtonError,
"NewtonMethod::result() called before NewtonMethod::solve()");
98 virtual void prepareStep(
Domain& solution)
100 _reassembled =
false;
101 if (_result.defect/_previousDefect > _reassembleThreshold){
102 if (_hangingNodeModifications){
103 auto dirichletValues = solution;
111 _gridOperator.localAssembler().backtransform(solution);
114 std::cout <<
" Reassembling matrix..." << std::endl;
115 *_jacobian =
Real(0.0);
116 _gridOperator.jacobian(solution, *_jacobian);
120 _linearReduction = _minLinearReduction;
121 if (not _fixedLinearReduction){
125 auto stop_defect =
max(_result.first_defect*_reduction, _absoluteLimit);
131 if (stop_defect/(10*_result.defect) > _result.defect*_result.defect/(_previousDefect*_previousDefect))
132 _linearReduction = stop_defect/(10*_result.defect);
134 _linearReduction =
min(_minLinearReduction, _result.defect*_result.defect/(_previousDefect*_previousDefect));
138 std::cout <<
" requested linear reduction: "
139 << std::setw(12) << std::setprecision(4) << std::scientific
140 << _linearReduction << std::endl;
143 virtual void linearSolve()
146 std::cout <<
" Solving linear system..." << std::endl;
149 Impl::setLinearSystemReuse(_linearSolver, not _reassembled);
153 _linearSolver.apply(*_jacobian, _correction, _residual, _linearReduction);
155 if (not _linearSolver.result().converged)
157 "NewtonMethod::linearSolve(): Linear solver did not converge "
158 "in " << _linearSolver.result().iterations <<
" iterations");
160 std::cout <<
" linear solver iterations: "
161 << std::setw(12) << _linearSolver.result().iterations << std::endl
162 <<
" linear defect reduction: "
163 << std::setw(12) << std::setprecision(4) << std::scientific
164 << _linearSolver.result().reduction << std::endl;
178 using Clock = std::chrono::steady_clock;
179 using Duration = Clock::duration;
180 auto assembler_time = Duration::zero();
181 auto linear_solver_time = Duration::zero();
182 auto line_search_time = Duration::zero();
184 [](Duration duration){
185 return std::chrono::duration<double>(duration).count();
187 auto start_solve = Clock::now();
193 _result.first_defect = _result.defect;
194 _previousDefect = _result.defect;
197 std::cout <<
" Initial defect: "
198 << std::setw(12) << std::setprecision(4) << std::scientific
199 << _result.defect << std::endl;
205 _jacobian = std::make_shared<Jacobian>(_gridOperator);
210 while (not _terminate->terminate()){
212 std::cout <<
" Newton iteration " << _result.iterations
213 <<
" --------------------------------" << std::endl;
218 auto start = Clock::now();
220 prepareStep(solution);
227 auto end = Clock::now();
228 assembler_time += end-start;
229 _result.assembler_time = to_seconds(assembler_time);
232 auto end = Clock::now();
233 assembler_time += end -start;
234 _result.assembler_time = to_seconds(assembler_time);
237 _previousDefect = _result.defect;
242 start = Clock::now();
250 linear_solver_time += end-start;
251 _result.linear_solver_time = to_seconds(linear_solver_time);
252 _result.linear_solver_iterations = _linearSolver.result().iterations;
256 linear_solver_time += end -start;
257 _result.linear_solver_time = to_seconds(linear_solver_time);
258 _result.linear_solver_iterations = _linearSolver.result().iterations;
263 start = Clock::now();
264 _lineSearch->lineSearch(solution, _correction);
267 line_search_time += end -start;
272 _result.reduction = _result.defect/_result.first_defect;
273 _result.iterations++;
274 _result.conv_rate = std::pow(_result.reduction, 1.0/_result.iterations);
277 std::cout <<
" linear solver time: "
278 << std::setw(12) << std::setprecision(4) << std::scientific
279 << to_seconds(linear_solver_time) << std::endl
280 <<
" defect reduction (this iteration):"
281 << std::setw(12) << std::setprecision(4) << std::scientific
282 << _result.defect/_previousDefect << std::endl
283 <<
" defect reduction (total): "
284 << std::setw(12) << std::setprecision(4) << std::scientific
285 << _result.reduction << std::endl
287 << std::setw(12) << std::setprecision(4) << std::scientific
288 << _result.defect << std::endl;
290 std::cout <<
" Newton iteration "
292 << _result.iterations
294 << std::setw(12) << std::setprecision(4) << std::scientific
296 <<
". Reduction (this): "
297 << std::setw(12) << std::setprecision(4) << std::scientific
298 << _result.defect/_previousDefect
299 <<
". Reduction (total): "
300 << std::setw(12) << std::setprecision(4) << std::scientific
305 auto end_solve = Clock::now();
306 auto solve_time = end_solve - start_solve;
307 _result.elapsed = to_seconds(solve_time);
310 std::cout <<
" Newton converged after "
312 << _result.iterations
313 <<
" iterations. Reduction: "
314 << std::setw(12) << std::setprecision(4) << std::scientific
316 <<
" (" << std::setprecision(4) << _result.elapsed <<
"s)"
326 if (_hangingNodeModifications){
327 auto dirichletValues = solution;
335 _gridOperator.localAssembler().backtransform(solution);
339 _gridOperator.residual(solution, _residual);
345 auto rankMax = Backend::native(_residual).infinity_norm();
346 _result.defect = _gridOperator.testGridFunctionSpace().gridView().comm().max(rankMax);
349 _result.defect = _linearSolver.norm(_residual);
355 if (_gridOperator.trialGridFunctionSpace().gridView().comm().rank()>0)
358 _verbosity = verbosity;
370 _reduction = reduction;
382 _absoluteLimit = absoluteLimit;
385 Real getAbsoluteLimit()
const
387 return _absoluteLimit;
405 _hangingNodeModifications = b;
431 _minLinearReduction = minLinearReduction;
441 _fixedLinearReduction = fixedLinearReduction;
452 _reassembleThreshold = reassembleThreshold;
485 _verbosity = parameterTree.
get(
"VerbosityLevel", _verbosity);
486 _reduction = parameterTree.
get(
"Reduction", _reduction);
487 _absoluteLimit = parameterTree.
get(
"AbsoluteLimit", _absoluteLimit);
488 _keepMatrix = parameterTree.
get(
"KeepMatrix", _keepMatrix);
489 _useMaxNorm = parameterTree.
get(
"UseMaxNorm", _useMaxNorm);
490 _hangingNodeModifications = parameterTree.
get(
"HangingNodeModifications", _hangingNodeModifications);
491 _minLinearReduction = parameterTree.
get(
"MinLinearReduction", _minLinearReduction);
492 _fixedLinearReduction = parameterTree.
get(
"FixedLinearReduction", _fixedLinearReduction);
493 _reassembleThreshold = parameterTree.
get(
"ReassembleThreshold", _reassembleThreshold);
496 std::string lineSearchStrategy = parameterTree.
get(
"LineSearchStrategy",
"hackbuschReusken");
497 auto strategy = lineSearchStrategyFromString(lineSearchStrategy);
498 _lineSearch = createLineSearch(*
this, strategy);
501 if (parameterTree.
hasSub(
"Terminate")){
502 _terminate->setParameters(parameterTree.
sub(
"Terminate"));
506 terminateTree[
"MaxIterations"] = std::to_string(parameterTree.
get(
"MaxIterations", 40));
507 terminateTree[
"ForceIteration"] = std::to_string(parameterTree.
get(
"ForceIteration",
false));
508 _terminate->setParameters(terminateTree);
510 if (parameterTree.
hasSub(
"LineSearch")){
511 _lineSearch->setParameters(parameterTree.
sub(
"LineSearch"));
515 lineSearchTree[
"MaxIterations"] = std::to_string(parameterTree.
get(
"LineSearchMaxIterations", 10));
516 lineSearchTree[
"DampingFactor"] = std::to_string(parameterTree.
get(
"LineSearchDampingFactor", 0.5));
517 lineSearchTree[
"AcceptBest"] = std::to_string(parameterTree.
get(
"LineSearchAcceptBest",
false));
518 _lineSearch->setParameters(lineSearchTree);
525 _terminate = terminate;
540 _lineSearch = lineSearch;
556 : _gridOperator(gridOperator)
557 , _linearSolver(linearSolver)
558 , _residual(gridOperator.testGridFunctionSpace())
559 , _correction(gridOperator.trialGridFunctionSpace())
561 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
562 _lineSearch = createLineSearch(*
this, LineSearchStrategy::hackbuschReusken);
570 : _gridOperator(gridOperator)
571 , _linearSolver(linearSolver)
572 , _residual(gridOperator.testGridFunctionSpace())
573 , _correction(gridOperator.trialGridFunctionSpace())
576 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
587 std::shared_ptr<Jacobian> _jacobian;
588 std::shared_ptr<Domain> _previousSolution;
590 std::shared_ptr<TerminateInterface> _terminate;
591 std::shared_ptr<LineSearch> _lineSearch;
595 bool _resultValid =
false;
596 Real _previousDefect = 0.0;
599 bool _reassembled =
true;
600 Real _linearReduction = 0.0;
603 unsigned int _verbosity = 0;
604 Real _reduction = 1e-8;
605 Real _absoluteLimit = 1e-12;
606 bool _keepMatrix =
true;
607 bool _useMaxNorm =
false;
610 bool _hangingNodeModifications =
false;
613 Real _minLinearReduction = 1e-3;
614 bool _fixedLinearReduction =
false;
615 Real _reassembleThreshold = 0.0;
Standard grid operator implementation.
Definition: gridoperator.hh:36
Abstract base class describing the line search interface.
Definition: linesearch.hh:15
Newton solver for solving non-linear problems.
Definition: newton.hh:64
typename GridOperator::Traits::Jacobian Jacobian
Type of the Jacobian matrix.
Definition: newton.hh:79
const Result & result() const
Return results.
Definition: newton.hh:91
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver, const ParameterTree ¶meterTree)
Construct Newton passing a parameter tree.
Definition: newton.hh:566
void setMinLinearReduction(Real minLinearReduction)
Set the minimal reduction in the linear solver.
Definition: newton.hh:429
void setTerminate(std::shared_ptr< TerminateInterface > terminate)
Set the termination criterion.
Definition: newton.hh:523
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: newton.hh:410
unsigned int getVerbosityLevel() const
Get verbosity level.
Definition: newton.hh:362
std::shared_ptr< TerminateInterface > getTerminate() const
Return a pointer to the stored termination criterion.
Definition: newton.hh:529
typename GridOperator::Traits::Domain Domain
Type of the domain (solution)
Definition: newton.hh:73
void setParameters(const ParameterTree ¶meterTree)
Interpret a parameter tree as a set of options for the newton solver.
Definition: newton.hh:484
void setHangingNodeModifications(bool b)
Does the problem have hanging nodes.
Definition: newton.hh:403
void setFixedLinearReduction(bool fixedLinearReduction)
Set wether to use a fixed reduction in the linear solver.
Definition: newton.hh:439
std::shared_ptr< LineSearch > getLineSearch() const
Return a pointer to the stored line search.
Definition: newton.hh:544
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: newton.hh:416
GridOperator_ GridOperator
Type of the grid operator.
Definition: newton.hh:67
void setLineSearch(std::shared_ptr< LineSearch > lineSearch)
Set the line search.
Definition: newton.hh:538
void setReduction(Real reduction)
Set reduction Newton needs to achieve.
Definition: newton.hh:368
typename Dune::FieldTraits< typename Domain::ElementType >::real_type Real
Number type.
Definition: newton.hh:82
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: newton.hh:391
LinearSolver_ LinearSolver
Type of the linear solver.
Definition: newton.hh:70
void setUseMaxNorm(bool b)
Set whether to use the maximum norm for stopping criteria.
Definition: newton.hh:397
typename GridOperator::Traits::Range Range
Type of the range (residual)
Definition: newton.hh:76
void setAbsoluteLimit(Real absoluteLimit)
Set absolute convergence limit.
Definition: newton.hh:380
virtual void updateDefect(Domain &solution)
Update _residual and defect in _result.
Definition: newton.hh:324
virtual void apply(Domain &solution)
Solve the nonlinear problem using solution as initial guess and for storing the result.
Definition: newton.hh:168
PDESolverResult< Real > Result
Type of results.
Definition: newton.hh:85
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver)
Construct Newton using default parameters with default parameters.
Definition: newton.hh:553
void setReassembleThreshold(Real reassembleThreshold)
Set a threshold, when the linear operator is reassembled.
Definition: newton.hh:450
void setVerbosityLevel(unsigned int verbosity)
Set how much output you get.
Definition: newton.hh:353
Real getReduction() const
Get reduction.
Definition: newton.hh:374
Hierarchical structure of string parameters.
Definition: parametertree.hh:35
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:183
ParameterTree & sub(const std::string &sub)
get substructure by name
Definition: parametertree.cc:101
bool hasSub(const std::string &sub) const
test for substructure
Definition: parametertree.cc:74
Utility class for storing and resetting stream attributes.
Definition: ios_state.hh:32
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
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
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:159
Utility class for storing and resetting stream attributes.
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