dune-pdelab 2.7-git
Loading...
Searching...
No Matches
solver/newton.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_PDELAB_SOLVER_NEWTON_HH
4#define DUNE_PDELAB_SOLVER_NEWTON_HH
5
6#include <dune/common/exceptions.hh>
7#include <dune/common/ios_state.hh>
8
13
14namespace Dune::PDELab
15{
16 namespace Impl
17 {
18 // Some SFinae magic to execute setReuse(bool) on a backend
19 template<typename T1, typename = void>
20 struct HasSetReuse
21 : std::false_type
22 {};
23
24 template<typename T>
25 struct HasSetReuse<T, decltype(std::declval<T>().setReuse(true), void())>
26 : std::true_type
27 {};
28
29 template<typename T>
30 inline void setLinearSystemReuse(T& solver_backend, bool reuse, std::true_type)
31 {
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);
35 }
36
37 template<typename T>
38 inline void setLinearSystemReuse(T&, bool, std::false_type)
39 {}
40
41 template<typename T>
42 inline void setLinearSystemReuse(T& solver_backend, bool reuse)
43 {
44 setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
45 }
46 }
47
48
61 template <typename GridOperator_, typename LinearSolver_>
63 {
64 public:
66 using GridOperator = GridOperator_;
67
69 using LinearSolver = LinearSolver_;
70
73
76
79
81 using Real = typename Dune::FieldTraits<typename Domain::ElementType>::real_type;
82
85
87 const Result& result() const
88 {
89 if (not _resultValid)
90 DUNE_THROW(NewtonError, "NewtonMethod::result() called before NewtonMethod::solve()");
91 return _result;
92 }
93
94 virtual void prepareStep(Domain& solution)
95 {
96 _reassembled = false;
97 if (_result.defect/_previousDefect > _reassembleThreshold){
98 if (_hangingNodeModifications){
99 auto dirichletValues = solution;
100 // Set all non dirichlet values to zero
101 Dune::PDELab::set_shifted_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, dirichletValues);
102 // Set all constrained DOFs to zero in solution
103 Dune::PDELab::set_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, solution);
104 // Copy correct Dirichlet values back into solution vector
105 Dune::PDELab::copy_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), dirichletValues, solution);
106 // Interpolate periodic constraints / hanging nodes
107 _gridOperator.localAssembler().backtransform(solution);
108 }
109 if (_verbosity>=3)
110 std::cout << " Reassembling matrix..." << std::endl;
111 *_jacobian = Real(0.0);
112 _gridOperator.jacobian(solution, *_jacobian);
113 _reassembled = true;
114 }
115
116 _linearReduction = _minLinearReduction;
117 if (not _fixedLinearReduction){
118 // Determine maximum defect, where Newton is converged.
119 using std::min;
120 using std::max;
121 auto stop_defect = max(_result.first_defect*_reduction, _absoluteLimit);
122
123 // To achieve second order convergence of newton we need a linear
124 // reduction of at least current_defect^2/prev_defect^2. For the
125 // last newton step a linear reduction of
126 // 1/10*end_defect/current_defect is sufficient for convergence.
127 if (stop_defect/(10*_result.defect) > _result.defect*_result.defect/(_previousDefect*_previousDefect))
128 _linearReduction = stop_defect/(10*_result.defect);
129 else
130 _linearReduction = min(_minLinearReduction, _result.defect*_result.defect/(_previousDefect*_previousDefect));
131 }
132
133 if (_verbosity >= 3)
134 std::cout << " requested linear reduction: "
135 << std::setw(12) << std::setprecision(4) << std::scientific
136 << _linearReduction << std::endl;
137 }
138
139 virtual void linearSolve()
140 {
141 if (_verbosity >= 4)
142 std::cout << " Solving linear system..." << std::endl;
143
144 // If the jacobian was not reassembled we might save some work in the solver backend
145 Impl::setLinearSystemReuse(_linearSolver, not _reassembled);
146
147 // Solve the linear system
148 _correction = 0.0;
149 _linearSolver.apply(*_jacobian, _correction, _residual, _linearReduction);
150
151 if (not _linearSolver.result().converged)
152 DUNE_THROW(NewtonLinearSolverError,
153 "NewtonMethod::linearSolve(): Linear solver did not converge "
154 "in " << _linearSolver.result().iterations << " iterations");
155 if (_verbosity >= 4)
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;
161 }
162
164 virtual void apply(Domain& solution)
165 {
166 // Reset solver statistics
167 _result.clear();
168 _resultValid = true;
169
170 // Store old ios flags (will be reset when this goes out of scope)
171 ios_base_all_saver restorer(std::cout);
172
173 // Prepare time measuring
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();
179 auto to_seconds =
180 [](Duration duration){
181 return std::chrono::duration<double>(duration).count();
182 };
183 auto start_solve = Clock::now();
184
185 //=========================
186 // Calculate initial defect
187 //=========================
188 updateDefect(solution);
189 _result.first_defect = _result.defect;
190 _previousDefect = _result.defect;
191
192 if (_verbosity >= 2)
193 std::cout << " Initial defect: "
194 << std::setw(12) << std::setprecision(4) << std::scientific
195 << _result.defect << std::endl;
196
197 //==========================
198 // Calculate Jacobian matrix
199 //==========================
200 if (not _jacobian)
201 _jacobian = std::make_shared<Jacobian>(_gridOperator);
202
203 //=========================
204 // Nonlinear iteration loop
205 //=========================
206 while (not _terminate->terminate()){
207 if(_verbosity >= 3)
208 std::cout << " Newton iteration " << _result.iterations
209 << " --------------------------------" << std::endl;
210
211 //=============
212 // Prepare step
213 //=============
214 auto start = Clock::now();
215 prepareStep(solution);
216 auto end = Clock::now();
217 assembler_time += end -start;
218
219 // Store defect
220 _previousDefect = _result.defect;
221
222 //====================
223 // Solve linear system
224 //====================
225 start = Clock::now();
226 linearSolve();
227 end = Clock::now();
228 linear_solver_time += end -start;
229
230 //===================================
231 // Do line search and update solution
232 //===================================
233 start = Clock::now();
234 _lineSearch->lineSearch(solution, _correction);
235 // lineSearch(solution);
236 end = Clock::now();
237 line_search_time += end -start;
238
239 //========================================
240 // Store statistics and create some output
241 //========================================
242 _result.reduction = _result.defect/_result.first_defect;
243 _result.iterations++;
244 _result.conv_rate = std::pow(_result.reduction, 1.0/_result.iterations);
245
246 if (_verbosity >= 3)
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
255 << _result.reduction << std::endl
256 << " new defect: "
257 << std::setw(12) << std::setprecision(4) << std::scientific
258 << _result.defect << std::endl;
259 if (_verbosity == 2)
260 std::cout << " Newton iteration "
261 << std::setw(2)
262 << _result.iterations
263 << ". New defect: "
264 << std::setw(12) << std::setprecision(4) << std::scientific
265 << _result.defect
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
271 << _result.reduction
272 << std::endl;
273 } // end while loop of nonlinear iterations
274
275 auto end_solve = Clock::now();
276 auto solve_time = end_solve - start_solve;
277 _result.elapsed = to_seconds(solve_time);
278
279 if (_verbosity == 1)
280 std::cout << " Newton converged after "
281 << std::setw(2)
282 << _result.iterations
283 << " iterations. Reduction: "
284 << std::setw(12) << std::setprecision(4) << std::scientific
285 << _result.reduction
286 << " (" << std::setprecision(4) << _result.elapsed << "s)"
287 << std::endl;
288
289 if (not _keepMatrix)
290 _jacobian.reset();
291 }
292
294 virtual void updateDefect(Domain& solution)
295 {
296 if (_hangingNodeModifications){
297 auto dirichletValues = solution;
298 // Set all non dirichlet values to zero
299 Dune::PDELab::set_shifted_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, dirichletValues);
300 // Set all constrained DOFs to zero in solution
301 Dune::PDELab::set_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, solution);
302 // Copy correct Dirichlet values back into solution vector
303 Dune::PDELab::copy_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), dirichletValues, solution);
304 // Interpolate periodic constraints / hanging nodes
305 _gridOperator.localAssembler().backtransform(solution);
306 }
307
308 _residual = 0.0;
309 _gridOperator.residual(solution, _residual);
310
311 // Use the maximum norm as a stopping criterion. This helps loosen the tolerance
312 // when solving for stationary solutions of nonlinear time-dependent problems.
313 // The default is to use the Euclidean norm (in the else-block) as before
314 if (_useMaxNorm){
315 auto rankMax = Backend::native(_residual).infinity_norm();
316 _result.defect = _gridOperator.testGridFunctionSpace().gridView().comm().max(rankMax);
317 }
318 else
319 _result.defect = _linearSolver.norm(_residual);
320 }
321
323 void setVerbosityLevel(unsigned int verbosity)
324 {
325 if (_gridOperator.trialGridFunctionSpace().gridView().comm().rank()>0)
326 _verbosity = 0;
327 else
328 _verbosity = verbosity;
329 }
330
332 unsigned int getVerbosityLevel() const
333 {
334 return _verbosity;
335 }
336
338 void setReduction(Real reduction)
339 {
340 _reduction = reduction;
341 }
342
345 {
346 return _reduction;
347 }
348
350 void setAbsoluteLimit(Real absoluteLimit)
351 {
352 _absoluteLimit = absoluteLimit;
353 }
354
356 {
357 return _absoluteLimit;
358 }
359
361 void setKeepMatrix(bool b)
362 {
363 _keepMatrix = b;
364 }
365
367 void setUseMaxNorm(bool b)
368 {
369 _useMaxNorm = b;
370 }
371
374 {
375 _hangingNodeModifications = b;
376 }
377
378
380 bool keepMatrix() const
381 {
382 return _keepMatrix;
383 }
384
387 {
388 if(_jacobian)
389 _jacobian.reset();
390 }
391
399 void setMinLinearReduction(Real minLinearReduction)
400 {
401 _minLinearReduction = minLinearReduction;
402 }
403
409 void setFixedLinearReduction(bool fixedLinearReduction)
410 {
411 _fixedLinearReduction = fixedLinearReduction;
412 }
413
420 void setReassembleThreshold(Real reassembleThreshold)
421 {
422 _reassembleThreshold = reassembleThreshold;
423 }
424
458 void setParameters(const ParameterTree& parameterTree){
459 if (parameterTree.hasKey("verbosity"))
460 setVerbosityLevel(parameterTree.get<unsigned int>("verbosity"));
461 if (parameterTree.hasKey("reduction"))
462 setReduction(parameterTree.get<Real>("reduction"));
463 if (parameterTree.hasKey("absolute_limit"))
464 setAbsoluteLimit(parameterTree.get<Real>("absolute_limit"));
465 if (parameterTree.hasKey("keeep_matrix"))
466 setKeepMatrix(parameterTree.get<bool>("keep_matrix"));
467 if (parameterTree.hasKey("use_max_norm"))
468 setUseMaxNorm(parameterTree.get<bool>("use_max_norm"));
469 if (parameterTree.hasKey("hanging_node_modifications"))
470 _hangingNodeModifications = parameterTree.get<bool>("hanging_node_modifications");
471
472 if (parameterTree.hasKey("min_linear_reduction"))
473 setMinLinearReduction(parameterTree.get<Real>("min_linear_reduction"));
474 if (parameterTree.hasKey("fixed_linear_reduction"))
475 setFixedLinearReduction(parameterTree.get<bool>("fixed_linear_reduction"));
476 if (parameterTree.hasKey("reassemble_threshold"))
477 setReassembleThreshold(parameterTree.get<Real>("reassemble_threshold"));
478
479 _terminate->setParameters(parameterTree.sub("terminate"));
480 _lineSearch->setParameters(parameterTree.sub("line_search"));
481 }
482
484 void setTerminate(std::shared_ptr<TerminateInterface> terminate)
485 {
486 _terminate = terminate;
487 }
488
493 void setLineSearch(std::shared_ptr<LineSearchInterface<Domain>> lineSearch)
494 {
495 _lineSearch = lineSearch;
496 }
497
500 const GridOperator& gridOperator,
501 LinearSolver& linearSolver,
502 const std::string& lineSearchStrategy="hackbusch_reusken")
503 : _gridOperator(gridOperator)
504 , _linearSolver(linearSolver)
505 , _residual(gridOperator.testGridFunctionSpace())
506 , _correction(gridOperator.trialGridFunctionSpace())
507 {
508 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
509 _lineSearch = getLineSearch(*this, lineSearchStrategy);
510 }
511
514 const GridOperator& gridOperator,
515 LinearSolver& linearSolver,
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())
522
523 {
524 setParameters(parameterTree);
525 _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
526 _lineSearch = getLineSearch(*this, lineSearchStrategy);
527 }
528
529 private:
530 const GridOperator& _gridOperator;
531 LinearSolver& _linearSolver;
532
533 // Vectors and Jacobi matrix we set up only once
534 Range _residual;
535 Domain _correction;
536 std::shared_ptr<Jacobian> _jacobian;
537 std::shared_ptr<Domain> _previousSolution;
538
539 std::shared_ptr<TerminateInterface> _terminate;
540 std::shared_ptr<LineSearchInterface<Domain>> _lineSearch;
541
542 // Class for storing results
543 Result _result;
544 bool _resultValid = false; // result class only valid after calling apply
545 Real _previousDefect = 0.0;
546
547 // Remember if jacobian was reassembled in prepareStep
548 bool _reassembled = true; // will be set in prepare step
549 Real _linearReduction = 0.0; // will be set in prepare step
550
551 // User parameters
552 unsigned int _verbosity = 0;
553 Real _reduction = 1e-8;
554 Real _absoluteLimit = 1e-12;
555 bool _keepMatrix = true;
556 bool _useMaxNorm = false;
557
558 // Special treatment if we have hanging nodes
559 bool _hangingNodeModifications = false;
560
561 // User parameters for prepareStep()
562 Real _minLinearReduction = 1e-3;
563 bool _fixedLinearReduction = false;
564 Real _reassembleThreshold = 0.0;
565 };
566
567} // namespace Dune::PDELab
568
569#endif
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
STL namespace.
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 &parameterTree, 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 &parameterTree)
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