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 = 1
e-8;
554 Real _absoluteLimit = 1
e-12;
555 bool _keepMatrix =
true;
556 bool _useMaxNorm =
false;
559 bool _hangingNodeModifications =
false;
562 Real _minLinearReduction = 1
e-3;
563 bool _fixedLinearReduction =
false;
564 Real _reassembleThreshold = 0.0;