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);
113 if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
115 std::cout <<
" Reassembling matrix..." << std::endl;
116 *_jacobian =
Real(0.0);
117 _gridOperator.jacobian(solution, *_jacobian);
122 _linearReduction = _minLinearReduction;
124 if (not _fixedLinearReduction && !(_result.iterations==0 && _result.first_defect<_absoluteLimit)){
128 auto stop_defect =
max(_result.first_defect*_reduction, _absoluteLimit);
135 max( stop_defect/(10*_result.defect),
136 min(_minLinearReduction, _result.defect*_result.defect/(_previousDefect*_previousDefect)) );
140 std::cout <<
" requested linear reduction: "
141 << std::setw(12) << std::setprecision(4) << std::scientific
142 << _linearReduction << std::endl;
145 virtual void linearSolve()
148 std::cout <<
" Solving linear system..." << std::endl;
150 if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
152 Impl::setLinearSystemReuse(_linearSolver, not _reassembled);
157 if constexpr (linearSolverIsMatrixFree<LinearSolver>()){
158 _linearSolver.apply(_correction, _residual, _linearReduction);
161 _linearSolver.apply(*_jacobian, _correction, _residual, _linearReduction);
164 if (not _linearSolver.result().converged)
166 "NewtonMethod::linearSolve(): Linear solver did not converge "
167 "in " << _linearSolver.result().iterations <<
" iterations");
169 std::cout <<
" linear solver iterations: "
170 << std::setw(12) << _linearSolver.result().iterations << std::endl
171 <<
" linear defect reduction: "
172 << std::setw(12) << std::setprecision(4) << std::scientific
173 << _linearSolver.result().reduction << std::endl;
187 using Clock = std::chrono::steady_clock;
188 using Duration = Clock::duration;
189 auto assembler_time = Duration::zero();
190 auto linear_solver_time = Duration::zero();
191 auto line_search_time = Duration::zero();
193 [](Duration duration){
194 return std::chrono::duration<double>(duration).count();
196 auto start_solve = Clock::now();
202 _result.first_defect = _result.defect;
203 _previousDefect = _result.defect;
206 std::cout <<
" Initial defect: "
207 << std::setw(12) << std::setprecision(4) << std::scientific
208 << _result.defect << std::endl;
213 if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
215 _jacobian = std::make_shared<Jacobian>(_gridOperator);
221 while (not _terminate->terminate()){
223 std::cout <<
" Newton iteration " << _result.iterations
224 <<
" --------------------------------" << std::endl;
229 auto start = Clock::now();
231 prepareStep(solution);
238 auto end = Clock::now();
239 assembler_time += end-start;
240 _result.assembler_time = to_seconds(assembler_time);
243 auto end = Clock::now();
244 assembler_time += end -start;
245 _result.assembler_time = to_seconds(assembler_time);
248 _previousDefect = _result.defect;
253 start = Clock::now();
256 if constexpr (linearSolverIsMatrixFree<LinearSolver>()){
257 _linearSolver.setLinearizationPoint(solution);
266 linear_solver_time += end-start;
267 _result.linear_solver_time = to_seconds(linear_solver_time);
268 _result.linear_solver_iterations += _linearSolver.result().iterations;
272 linear_solver_time += end -start;
273 _result.linear_solver_time = to_seconds(linear_solver_time);
274 _result.linear_solver_iterations += _linearSolver.result().iterations;
279 start = Clock::now();
280 _lineSearch->lineSearch(solution, _correction);
283 line_search_time += end -start;
288 _result.reduction = _result.defect/_result.first_defect;
289 _result.iterations++;
290 _result.conv_rate = std::pow(_result.reduction, 1.0/_result.iterations);
293 std::cout <<
" linear solver time: "
294 << std::setw(12) << std::setprecision(4) << std::scientific
295 << to_seconds(linear_solver_time) << std::endl
296 <<
" defect reduction (this iteration):"
297 << std::setw(12) << std::setprecision(4) << std::scientific
298 << _result.defect/_previousDefect << std::endl
299 <<
" defect reduction (total): "
300 << std::setw(12) << std::setprecision(4) << std::scientific
301 << _result.reduction << std::endl
303 << std::setw(12) << std::setprecision(4) << std::scientific
304 << _result.defect << std::endl;
306 std::cout <<
" Newton iteration "
308 << _result.iterations
310 << std::setw(12) << std::setprecision(4) << std::scientific
312 <<
". Reduction (this): "
313 << std::setw(12) << std::setprecision(4) << std::scientific
314 << _result.defect/_previousDefect
315 <<
". Reduction (total): "
316 << std::setw(12) << std::setprecision(4) << std::scientific
321 auto end_solve = Clock::now();
322 auto solve_time = end_solve - start_solve;
323 _result.elapsed = to_seconds(solve_time);
326 std::cout <<
" Newton converged after "
328 << _result.iterations
329 <<
" iterations. Reduction: "
330 << std::setw(12) << std::setprecision(4) << std::scientific
332 <<
" (" << std::setprecision(4) << _result.elapsed <<
"s)"
335 if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
344 if (_hangingNodeModifications){
345 auto dirichletValues = solution;
353 _gridOperator.localAssembler().backtransform(solution);
357 _gridOperator.residual(solution, _residual);
363 auto rankMax = Backend::native(_residual).infinity_norm();
364 _result.defect = _gridOperator.testGridFunctionSpace().gridView().comm().max(rankMax);
367 _result.defect = _linearSolver.norm(_residual);
373 if (_gridOperator.trialGridFunctionSpace().gridView().comm().rank()>0)
376 _verbosity = verbosity;
388 _reduction = reduction;
400 _absoluteLimit = absoluteLimit;
403 Real getAbsoluteLimit()
const
405 return _absoluteLimit;
423 _hangingNodeModifications = b;
449 _minLinearReduction = minLinearReduction;
459 _fixedLinearReduction = fixedLinearReduction;
470 _reassembleThreshold = reassembleThreshold;
503 _verbosity = parameterTree.
get(
"VerbosityLevel", _verbosity);
504 _reduction = parameterTree.
get(
"Reduction", _reduction);
505 _absoluteLimit = parameterTree.
get(
"AbsoluteLimit", _absoluteLimit);
506 _keepMatrix = parameterTree.
get(
"KeepMatrix", _keepMatrix);
507 _useMaxNorm = parameterTree.
get(
"UseMaxNorm", _useMaxNorm);
508 _hangingNodeModifications = parameterTree.
get(
"HangingNodeModifications", _hangingNodeModifications);
509 _minLinearReduction = parameterTree.
get(
"MinLinearReduction", _minLinearReduction);
510 _fixedLinearReduction = parameterTree.
get(
"FixedLinearReduction", _fixedLinearReduction);
511 _reassembleThreshold = parameterTree.
get(
"ReassembleThreshold", _reassembleThreshold);
514 std::string lineSearchStrategy = parameterTree.
get(
"LineSearchStrategy",
"hackbuschReusken");
515 auto strategy = lineSearchStrategyFromString(lineSearchStrategy);
516 _lineSearch = createLineSearch(*
this, strategy);
519 if (parameterTree.
hasSub(
"Terminate")){
520 _terminate->setParameters(parameterTree.
sub(
"Terminate"));
524 terminateTree[
"MaxIterations"] = std::to_string(parameterTree.
get(
"MaxIterations", 40));
525 terminateTree[
"ForceIteration"] = std::to_string(parameterTree.
get(
"ForceIteration",
false));
526 _terminate->setParameters(terminateTree);
528 if (parameterTree.
hasSub(
"LineSearch")){
529 _lineSearch->setParameters(parameterTree.
sub(
"LineSearch"));
533 lineSearchTree[
"MaxIterations"] = std::to_string(parameterTree.
get(
"LineSearchMaxIterations", 10));
534 lineSearchTree[
"DampingFactor"] = std::to_string(parameterTree.
get(
"LineSearchDampingFactor", 0.5));
535 lineSearchTree[
"AcceptBest"] = std::to_string(parameterTree.
get(
"LineSearchAcceptBest",
false));
536 _lineSearch->setParameters(lineSearchTree);
543 _terminate = terminate;
558 _lineSearch = lineSearch;
576 auto ioflags = std::cout.flags();
577 std::cout.setf(std::ios_base::boolalpha);
578 std::cout << _name <<
" parameters:\n";
579 std::cout <<
"Verbosity............... " << _verbosity << std::endl;
580 std::cout <<
"Reduction............... " << _reduction << std::endl;
581 std::cout <<
"AbsoluteLimit........... " << _absoluteLimit << std::endl;
582 std::cout <<
"KeepMatrix.............. " << _keepMatrix << std::endl;
583 std::cout <<
"UseMaxNorm.............. " << _useMaxNorm << std::endl;
584 std::cout <<
"MinLinearReduction...... " << _minLinearReduction << std::endl;
585 std::cout <<
"FixedLinearReduction.... " << _fixedLinearReduction << std::endl;
586 std::cout <<
"ReassembleThreshold..... " << _reassembleThreshold << std::endl;
587 std::cout <<
"HangingNodeModifications " << _hangingNodeModifications << std::endl;
590 _terminate->printParameters();
591 _lineSearch->printParameters();
596 std::cout.flags(ioflags);
599 std::cout.flags(ioflags);
609 : _gridOperator(gridOperator)
610 , _linearSolver(linearSolver)
611 , _residual(gridOperator.testGridFunctionSpace())
612 , _correction(gridOperator.trialGridFunctionSpace())
614 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
615 _lineSearch = createLineSearch(*
this, LineSearchStrategy::hackbuschReusken);
623 : _gridOperator(gridOperator)
624 , _linearSolver(linearSolver)
625 , _residual(gridOperator.testGridFunctionSpace())
626 , _correction(gridOperator.trialGridFunctionSpace())
629 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
642 std::shared_ptr<Jacobian> _jacobian;
644 std::shared_ptr<TerminateInterface> _terminate;
645 std::shared_ptr<LineSearch> _lineSearch;
649 bool _resultValid =
false;
650 Real _previousDefect = 0.0;
653 bool _reassembled =
true;
654 Real _linearReduction = 0.0;
657 unsigned int _verbosity = 0;
658 Real _reduction = 1e-8;
659 Real _absoluteLimit = 1e-12;
660 bool _keepMatrix =
true;
661 bool _useMaxNorm =
false;
664 bool _hangingNodeModifications =
false;
667 Real _minLinearReduction = 1e-3;
668 bool _fixedLinearReduction =
false;
669 Real _reassembleThreshold = 0.0;
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
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver, const ParameterTree ¶meterTree)
Construct Newton passing a parameter tree.
Definition: newton.hh:619
void setMinLinearReduction(Real minLinearReduction)
Set the minimal reduction in the linear solver.
Definition: newton.hh:447
void setTerminate(std::shared_ptr< TerminateInterface > terminate)
Set the termination criterion.
Definition: newton.hh:541
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: newton.hh:428
unsigned int getVerbosityLevel() const
Get verbosity level.
Definition: newton.hh:380
std::shared_ptr< TerminateInterface > getTerminate() const
Return a pointer to the stored termination criterion.
Definition: newton.hh:547
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:502
void setHangingNodeModifications(bool b)
Does the problem have hanging nodes.
Definition: newton.hh:421
void setFixedLinearReduction(bool fixedLinearReduction)
Set wether to use a fixed reduction in the linear solver.
Definition: newton.hh:457
std::shared_ptr< LineSearch > getLineSearch() const
Return a pointer to the stored line search.
Definition: newton.hh:562
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: newton.hh:434
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:556
void setReduction(Real reduction)
Set reduction Newton needs to achieve.
Definition: newton.hh:386
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:409
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:415
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:398
virtual void updateDefect(Domain &solution)
Update _residual and defect in _result.
Definition: newton.hh:342
void printParameters(const std::string &_name="NewtonMethod") const
Output NewtonMethod parameters.
Definition: newton.hh:572
virtual void apply(Domain &solution)
Solve the nonlinear problem using solution as initial guess and for storing the result.
Definition: newton.hh:177
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:606
Result & result()
Return results.
Definition: newton.hh:91
void setReassembleThreshold(Real reassembleThreshold)
Set a threshold, when the linear operator is reassembled.
Definition: newton.hh:468
void setVerbosityLevel(unsigned int verbosity)
Set how much output you get.
Definition: newton.hh:371
Real getReduction() const
Get reduction.
Definition: newton.hh:392
Hierarchical structure of string parameters.
Definition: parametertree.hh:37
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:188
ParameterTree & sub(const std::string &sub)
get substructure by name
Definition: parametertree.cc:106
bool hasSub(const std::string &sub) const
test for substructure
Definition: parametertree.cc:79
Utility class for storing and resetting stream attributes.
Definition: ios_state.hh:34
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:162
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