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>
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;
94 DUNE_THROW(NewtonError,
"NewtonMethod::result() called before NewtonMethod::solve()");
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;
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;
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)
165 DUNE_THROW(NewtonLinearSolverError,
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;
184 ios_base_all_saver restorer(std::cout);
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();
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();
238 auto end = Clock::now();
239 assembler_time += end-start;
243 auto end = Clock::now();
244 assembler_time += end -start;
248 _previousDefect = _result.
defect;
253 start = Clock::now();
256 if constexpr (linearSolverIsMatrixFree<LinearSolver>()){
257 _linearSolver.setLinearizationPoint(solution);
266 linear_solver_time += end-start;
272 linear_solver_time += end -start;
279 start = Clock::now();
280 _lineSearch->lineSearch(solution, _correction);
283 line_search_time += end -start;
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
303 << std::setw(12) << std::setprecision(4) << std::scientific
304 << _result.
defect << std::endl;
306 std::cout <<
" Newton iteration "
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 "
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);
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;
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);
519 if (parameterTree.hasSub(
"Terminate")){
520 _terminate->setParameters(parameterTree.sub(
"Terminate"));
523 ParameterTree terminateTree;
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"));
532 ParameterTree lineSearchTree;
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);
622 const ParameterTree& parameterTree)
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;
643 std::shared_ptr<Domain> _previousSolution;
645 std::shared_ptr<TerminateInterface> _terminate;
646 std::shared_ptr<LineSearch> _lineSearch;
650 bool _resultValid =
false;
651 Real _previousDefect = 0.0;
654 bool _reassembled =
true;
655 Real _linearReduction = 0.0;
658 unsigned int _verbosity = 0;
659 Real _reduction = 1e-8;
660 Real _absoluteLimit = 1e-12;
661 bool _keepMatrix =
true;
662 bool _useMaxNorm =
false;
665 bool _hangingNodeModifications =
false;
668 Real _minLinearReduction = 1e-3;
669 bool _fixedLinearReduction =
false;
670 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 Solver::Domain > > createLineSearch(Solver &solver, LineSearchStrategy strategy)
fectory function to create an instace of a line-search
Definition: linesearch.hh:254
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
Abstract base class describing the line search interface.
Definition: linesearch.hh:15
Newton solver for solving non-linear problems.
Definition: solver/newton.hh:64
typename GridOperator::Traits::Jacobian Jacobian
Type of the Jacobian matrix.
Definition: solver/newton.hh:79
const Result & result() const
Return results.
Definition: solver/newton.hh:91
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver, const ParameterTree ¶meterTree)
Construct Newton passing a parameter tree.
Definition: solver/newton.hh:619
virtual ~NewtonMethod()
Definition: solver/newton.hh:633
void setMinLinearReduction(Real minLinearReduction)
Set the minimal reduction in the linear solver.
Definition: solver/newton.hh:447
void setTerminate(std::shared_ptr< TerminateInterface > terminate)
Set the termination criterion.
Definition: solver/newton.hh:541
virtual void linearSolve()
Definition: solver/newton.hh:145
Real getAbsoluteLimit() const
Definition: solver/newton.hh:403
bool keepMatrix() const
Return whether the jacobian matrix is kept across calls to apply().
Definition: solver/newton.hh:428
unsigned int getVerbosityLevel() const
Get verbosity level.
Definition: solver/newton.hh:380
std::shared_ptr< TerminateInterface > getTerminate() const
Return a pointer to the stored termination criterion.
Definition: solver/newton.hh:547
typename GridOperator::Traits::Domain Domain
Type of the domain (solution)
Definition: solver/newton.hh:73
void setParameters(const ParameterTree ¶meterTree)
Interpret a parameter tree as a set of options for the newton solver.
Definition: solver/newton.hh:502
virtual void prepareStep(Domain &solution)
Definition: solver/newton.hh:98
void setHangingNodeModifications(bool b)
Does the problem have hanging nodes.
Definition: solver/newton.hh:421
void setFixedLinearReduction(bool fixedLinearReduction)
Set wether to use a fixed reduction in the linear solver.
Definition: solver/newton.hh:457
std::shared_ptr< LineSearch > getLineSearch() const
Return a pointer to the stored line search.
Definition: solver/newton.hh:562
void discardMatrix()
Discard the stored Jacobian matrix.
Definition: solver/newton.hh:434
GridOperator_ GridOperator
Type of the grid operator.
Definition: solver/newton.hh:67
void setLineSearch(std::shared_ptr< LineSearch > lineSearch)
Set the line search.
Definition: solver/newton.hh:556
void setReduction(Real reduction)
Set reduction Newton needs to achieve.
Definition: solver/newton.hh:386
typename Dune::FieldTraits< typename Domain::ElementType >::real_type Real
Number type.
Definition: solver/newton.hh:82
void setKeepMatrix(bool b)
Set whether the jacobian matrix should be kept across calls to apply().
Definition: solver/newton.hh:409
LinearSolver_ LinearSolver
Type of the linear solver.
Definition: solver/newton.hh:70
void setUseMaxNorm(bool b)
Set whether to use the maximum norm for stopping criteria.
Definition: solver/newton.hh:415
typename GridOperator::Traits::Range Range
Type of the range (residual)
Definition: solver/newton.hh:76
void setAbsoluteLimit(Real absoluteLimit)
Set absolute convergence limit.
Definition: solver/newton.hh:398
virtual void updateDefect(Domain &solution)
Update _residual and defect in _result.
Definition: solver/newton.hh:342
void printParameters(const std::string &_name="NewtonMethod") const
Output NewtonMethod parameters.
Definition: solver/newton.hh:572
virtual void apply(Domain &solution)
Solve the nonlinear problem using solution as initial guess and for storing the result.
Definition: solver/newton.hh:177
PDESolverResult< Real > Result
Type of results.
Definition: solver/newton.hh:85
NewtonMethod(const GridOperator &gridOperator, LinearSolver &linearSolver)
Construct Newton using default parameters with default parameters.
Definition: solver/newton.hh:606
void setReassembleThreshold(Real reassembleThreshold)
Set a threshold, when the linear operator is reassembled.
Definition: solver/newton.hh:468
void setVerbosityLevel(unsigned int verbosity)
Set how much output you get.
Definition: solver/newton.hh:371
Real getReduction() const
Get reduction.
Definition: solver/newton.hh:392
RFType defect
Definition: solver/utility.hh:11
double assembler_time
Definition: solver/utility.hh:12
int linear_solver_iterations
Definition: solver/utility.hh:14
double linear_solver_time
Definition: solver/utility.hh:13
RFType first_defect
Definition: solver/utility.hh:10
void clear()
Definition: solver/utility.hh:21