DUNE PDELab (git)

newton.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4//
5// Note: This Newton is deprecated and will be removed for PDELab 2.8. Please
6// use the Newton from <dune/pdelab/solver/newton.hh>
7//
8#ifndef DUNE_PDELAB_NEWTON_NEWTON_HH
9#define DUNE_PDELAB_NEWTON_NEWTON_HH
10
11#include <iostream>
12#include <iomanip>
13#include <cmath>
14#include <memory>
15
16#include <type_traits>
17
18#include <math.h>
19
23#include <dune/common/timer.hh>
25
26#include <dune/pdelab/backend/solver.hh>
27#include <dune/pdelab/constraints/common/constraints.hh>
28#include <dune/pdelab/solver/newton.hh>
29
30
31// Note: Before introducing the new Newton there was no doxygen documentation
32// for Newton at all. This means that it doesn't make any sense to have doxygen
33// documentation for the old Newton now.
34#ifndef DOXYGEN
35
36namespace Dune
37{
38 namespace PDELab
39 {
40
41 class NewtonLineSearchError : public NewtonError {};
42 class NewtonNotConverged : public NewtonError {};
43
44 // Status information of Newton's method
45 template<class RFType>
46 struct NewtonResult : LinearSolverResult<RFType>
47 {
48 RFType first_defect; // the first defect
49 RFType defect; // the final defect
50 double assembler_time; // Cumulative time for matrix assembly
51 double linear_solver_time; // Cumulative time for linear solver
52 int linear_solver_iterations; // Total number of linear iterations
53
54 NewtonResult() :
55 first_defect(0.0), defect(0.0), assembler_time(0.0), linear_solver_time(0.0),
56 linear_solver_iterations(0) {}
57 };
58
59 template<class GOS, class TrlV, class TstV>
60 class NewtonBase
61 {
62 typedef GOS GridOperator;
63 typedef TrlV TrialVector;
64 typedef TstV TestVector;
65
66 typedef typename TestVector::ElementType RFType;
67 typedef typename GOS::Traits::Jacobian Matrix;
68
69
70 public:
71 // export result type
72 typedef NewtonResult<RFType> Result;
73
74 void setVerbosityLevel(unsigned int verbosity_level)
75 {
76 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
77 verbosity_level_ = 0;
78 else
79 verbosity_level_ = verbosity_level;
80 }
81
83 void setKeepMatrix(bool b)
84 {
85 keep_matrix_ = b;
86 }
87
89 bool keepMatrix() const
90 {
91 return keep_matrix_;
92 }
93
95 void discardMatrix()
96 {
97 if(A_)
98 A_.reset();
99 }
100
101 protected:
102 const GridOperator& gridoperator_;
103 TrialVector *u_;
104 std::shared_ptr<TrialVector> z_;
105 std::shared_ptr<TestVector> r_;
106 std::shared_ptr<Matrix> A_;
107 Result res_;
108 unsigned int verbosity_level_;
109 RFType prev_defect_;
110 RFType linear_reduction_;
111 bool reassembled_;
112 RFType reduction_;
113 RFType abs_limit_;
114 bool keep_matrix_;
115
116 NewtonBase(const GridOperator& go, TrialVector& u)
117 : gridoperator_(go)
118 , u_(&u)
119 , verbosity_level_(1)
120 , keep_matrix_(true)
121 {
122 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
123 verbosity_level_ = 0;
124 }
125
126 NewtonBase(const GridOperator& go)
127 : gridoperator_(go)
128 , u_(0)
129 , verbosity_level_(1)
130 , keep_matrix_(true)
131 {
132 if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
133 verbosity_level_ = 0;
134 }
135
136 virtual ~NewtonBase() { }
137
138 virtual bool terminate() = 0;
139 virtual void prepare_step(Matrix& A, TestVector& r) = 0;
140 virtual void line_search(TrialVector& z, TestVector& r) = 0;
141 virtual void defect(TestVector& r) = 0;
142 }; // end class NewtonBase
143
144 template<class GOS, class S, class TrlV, class TstV>
145 class NewtonSolver : public virtual NewtonBase<GOS,TrlV,TstV>
146 {
147 typedef S Solver;
148 typedef GOS GridOperator;
149 typedef TrlV TrialVector;
150 typedef TstV TestVector;
151
152 typedef typename TestVector::ElementType RFType;
153 typedef typename GOS::Traits::Jacobian Matrix;
154
155 public:
156 typedef NewtonResult<RFType> Result;
157
158 NewtonSolver(const GridOperator& go, TrialVector& u_, Solver& solver)
159 : NewtonBase<GOS,TrlV,TstV>(go,u_)
160 , solver_(solver)
161 , result_valid_(false)
162 , use_maxnorm_(false)
163 {}
164
165 NewtonSolver(const GridOperator& go, Solver& solver)
166 : NewtonBase<GOS,TrlV,TstV>(go)
167 , solver_(solver)
168 , result_valid_(false)
169 , use_maxnorm_(false)
170 , hanging_node_modifications_(false)
171 {}
172
173 void apply();
174
175 void apply(TrialVector& u_);
176
177 const Result& result() const
178 {
179 if (!result_valid_)
180 DUNE_THROW(NewtonError,
181 "NewtonSolver::result() called before NewtonSolver::solve()");
182 return this->res_;
183 }
184
186 void setUseMaxNorm(bool b)
187 {
188 use_maxnorm_ = b;
189 }
190
191 void setHangingNodeModifications(bool b)
192 {
193 hanging_node_modifications_ = b;
194 }
195
196 protected:
197 virtual void defect(TestVector& r) override
198 {
199 if (hanging_node_modifications_){
200 auto dirichletValues = *this->u_;
201 // Set all non dirichlet values to zero
202 Dune::PDELab::set_shifted_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, dirichletValues);
203 // Set all constrained DOFs to zero in solution
204 Dune::PDELab::set_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, *this->u_);
205 // Copy correct Dirichlet values back into solution vector
206 Dune::PDELab::copy_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), dirichletValues, *this->u_);
207 // Interpolate periodic constraints / hanging nodes
208 this->gridoperator_.localAssembler().backtransform(*this->u_);
209 }
210
211 r = 0.0;
212 this->gridoperator_.residual(*this->u_, r);
213 // use maximum norm as a stopping criterion, this helps loosen the tolerance
214 // when solving for stationary solutions of nonlinear time-dependent problems
215 // the default is to use the Euclidean norm (in the else-block) as before
216 if(use_maxnorm_)
217 {
218 auto rank_max = Backend::native(r).infinity_norm();
219 this->res_.defect = this->gridoperator_.testGridFunctionSpace().gridView().comm().max(rank_max);
220 }
221 else
222 {
223 this->res_.defect = this->solver_.norm(r);
224 }
225 if (!std::isfinite(this->res_.defect))
226 DUNE_THROW(NewtonDefectError,
227 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
228 }
229
230
231 private:
232 void linearSolve(Matrix& A, TrialVector& z, TestVector& r) const
233 {
234 if (this->verbosity_level_ >= 4)
235 std::cout << " Solving linear system..." << std::endl;
236 z = 0.0;
237 // If possible tell solver backend to reuse linear system when we did not reassemble.
238 Impl::setLinearSystemReuse(this->solver_, !this->reassembled_);
239 this->solver_.apply(A, z, r, this->linear_reduction_);
240
241 ios_base_all_saver restorer(std::cout); // store old ios flags
242
243 if (!this->solver_.result().converged)
244 DUNE_THROW(NewtonLinearSolverError,
245 "NewtonSolver::linearSolve(): Linear solver did not converge "
246 "in " << this->solver_.result().iterations << " iterations");
247 if (this->verbosity_level_ >= 4)
248 std::cout << " linear solver iterations: "
249 << std::setw(12) << solver_.result().iterations << std::endl
250 << " linear defect reduction: "
251 << std::setw(12) << std::setprecision(4) << std::scientific
252 << solver_.result().reduction << std::endl;
253 }
254
255 Solver& solver_;
256 bool result_valid_;
257 bool use_maxnorm_;
258 bool hanging_node_modifications_;
259 }; // end class NewtonSolver
260
261 template<class GOS, class S, class TrlV, class TstV>
262 void NewtonSolver<GOS,S,TrlV,TstV>::apply(TrialVector& u)
263 {
264 this->u_ = &u;
265 apply();
266 }
267
268 template<class GOS, class S, class TrlV, class TstV>
269 void NewtonSolver<GOS,S,TrlV,TstV>::apply()
270 {
271 this->res_.iterations = 0;
272 this->res_.converged = false;
273 this->res_.reduction = 1.0;
274 this->res_.conv_rate = 1.0;
275 this->res_.elapsed = 0.0;
276 this->res_.assembler_time = 0.0;
277 this->res_.linear_solver_time = 0.0;
278 this->res_.linear_solver_iterations = 0;
279 result_valid_ = true;
280 Timer timer;
281
282 try
283 {
284 if(!this->r_) {
285 // std::cout << "=== Setting up residual vector ..." << std::endl;
286 this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
287 }
288 // residual calculation in member function "defect":
289 //--------------------------------------------------
290 // - set residual vector to zero
291 // - calculate new residual
292 // - store norm of residual in "this->res_.defect"
293 this->defect(*this->r_);
294 this->res_.first_defect = this->res_.defect;
295 this->prev_defect_ = this->res_.defect;
296
297 if (this->verbosity_level_ >= 2)
298 {
299 // store old ios flags
300 ios_base_all_saver restorer(std::cout);
301 std::cout << " Initial defect: "
302 << std::setw(12) << std::setprecision(4) << std::scientific
303 << this->res_.defect << std::endl;
304 }
305
306 if(!this->A_) {
307 // std::cout << "==== Setting up jacobian matrix ... " << std::endl;
308 this->A_ = std::make_shared<Matrix>(this->gridoperator_);
309 }
310 if(!this->z_) {
311 // std::cout << "==== Setting up correction vector ... " << std::endl;
312 this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
313 }
314
315 while (!this->terminate())
316 {
317 if (this->verbosity_level_ >= 3)
318 std::cout << " Newton iteration " << this->res_.iterations
319 << " --------------------------------" << std::endl;
320
321 Timer assembler_timer;
322 try
323 {
324 // jacobian calculation in member function "prepare_step"
325 //-------------------------------------------------------
326 // - if above reassemble threshold
327 // - set jacobian to zero
328 // - calculate new jacobian
329 // - set linear reduction
330 this->prepare_step(*this->A_,*this->r_);
331 }
332 catch (...)
333 {
334 this->res_.assembler_time += assembler_timer.elapsed();
335 throw;
336 }
337 double assembler_time = assembler_timer.elapsed();
338 this->res_.assembler_time += assembler_time;
339 if (this->verbosity_level_ >= 3)
340 std::cout << " matrix assembly time: "
341 << std::setw(12) << std::setprecision(4) << std::scientific
342 << assembler_time << std::endl;
343
344 Timer linear_solver_timer;
345 try
346 {
347 // solution of linear system in member function "linearSolve"
348 //-----------------------------------------------------------
349 // - set initial guess for correction z to zero
350 // - call linear solver
351 this->linearSolve(*this->A_, *this->z_, *this->r_);
352 }
353 catch (...)
354 {
355 this->res_.linear_solver_time += linear_solver_timer.elapsed();
356 this->res_.linear_solver_iterations += this->solver_.result().iterations;
357 throw;
358 }
359 double linear_solver_time = linear_solver_timer.elapsed();
360 this->res_.linear_solver_time += linear_solver_time;
361 this->res_.linear_solver_iterations += this->solver_.result().iterations;
362
363 try
364 {
365 // line search with correction z
366 // the undamped version is also integrated in here
367 this->line_search(*this->z_, *this->r_);
368 }
369 catch (NewtonLineSearchError&)
370 {
371 if (this->reassembled_)
372 throw;
373 if (this->verbosity_level_ >= 3)
374 std::cout << " line search failed - trying again with reassembled matrix" << std::endl;
375 continue;
376 }
377
378 this->res_.reduction = this->res_.defect/this->res_.first_defect;
379 this->res_.iterations++;
380 this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
381
382 // store old ios flags
383 ios_base_all_saver restorer(std::cout);
384
385 if (this->verbosity_level_ >= 3)
386 std::cout << " linear solver time: "
387 << std::setw(12) << std::setprecision(4) << std::scientific
388 << linear_solver_time << std::endl
389 << " defect reduction (this iteration):"
390 << std::setw(12) << std::setprecision(4) << std::scientific
391 << this->res_.defect/this->prev_defect_ << std::endl
392 << " defect reduction (total): "
393 << std::setw(12) << std::setprecision(4) << std::scientific
394 << this->res_.reduction << std::endl
395 << " new defect: "
396 << std::setw(12) << std::setprecision(4) << std::scientific
397 << this->res_.defect << std::endl;
398 if (this->verbosity_level_ == 2)
399 std::cout << " Newton iteration " << std::setw(2) << this->res_.iterations
400 << ". New defect: "
401 << std::setw(12) << std::setprecision(4) << std::scientific
402 << this->res_.defect
403 << ". Reduction (this): "
404 << std::setw(12) << std::setprecision(4) << std::scientific
405 << this->res_.defect/this->prev_defect_
406 << ". Reduction (total): "
407 << std::setw(12) << std::setprecision(4) << std::scientific
408 << this->res_.reduction << std::endl;
409 } // end while
410 } // end try
411 catch(...)
412 {
413 this->res_.elapsed = timer.elapsed();
414 throw;
415 }
416 this->res_.elapsed = timer.elapsed();
417
418 ios_base_all_saver restorer(std::cout); // store old ios flags
419
420 if (this->verbosity_level_ == 1)
421 std::cout << " Newton converged after " << std::setw(2) << this->res_.iterations
422 << " iterations. Reduction: "
423 << std::setw(12) << std::setprecision(4) << std::scientific
424 << this->res_.reduction
425 << " (" << std::setprecision(4) << this->res_.elapsed << "s)"
426 << std::endl;
427
428 if(!this->keep_matrix_)
429 this->A_.reset();
430 } // end apply
431
432 template<class GOS, class TrlV, class TstV>
433 class NewtonTerminate : public virtual NewtonBase<GOS,TrlV,TstV>
434 {
435 typedef GOS GridOperator;
436 typedef TrlV TrialVector;
437
438 typedef typename TstV::ElementType RFType;
439
440 public:
441 NewtonTerminate(const GridOperator& go, TrialVector& u_)
442 : NewtonBase<GOS,TrlV,TstV>(go,u_)
443 , maxit_(40)
444 , force_iteration_(false)
445 {
446 this->reduction_ = 1e-8;
447 this->abs_limit_ = 1e-12;
448 }
449
450 NewtonTerminate(const GridOperator& go)
451 : NewtonBase<GOS,TrlV,TstV>(go)
452 , maxit_(40)
453 , force_iteration_(false)
454 {
455 this->reduction_ = 1e-8;
456 this->abs_limit_ = 1e-12;
457 }
458
459 void setReduction(RFType reduction)
460 {
461 this->reduction_ = reduction;
462 }
463
464 void setMaxIterations(unsigned int maxit)
465 {
466 maxit_ = maxit;
467 }
468
469 void setForceIteration(bool force_iteration)
470 {
471 force_iteration_ = force_iteration;
472 }
473
474 void setAbsoluteLimit(RFType abs_limit_)
475 {
476 this->abs_limit_ = abs_limit_;
477 }
478
479 virtual bool terminate() override
480 {
481 if (force_iteration_ && this->res_.iterations == 0)
482 return false;
483 this->res_.converged = this->res_.defect < this->abs_limit_
484 || this->res_.defect < this->res_.first_defect * this->reduction_;
485 if (this->res_.iterations >= maxit_ && !this->res_.converged)
486 DUNE_THROW(NewtonNotConverged,
487 "NewtonTerminate::terminate(): Maximum iteration count reached");
488 return this->res_.converged;
489 }
490
491 private:
492 unsigned int maxit_;
493 bool force_iteration_;
494 }; // end class NewtonTerminate
495
496 template<class GOS, class TrlV, class TstV>
497 class NewtonPrepareStep : public virtual NewtonBase<GOS,TrlV,TstV>
498 {
499 typedef GOS GridOperator;
500 typedef TrlV TrialVector;
501
502 typedef typename TstV::ElementType RFType;
503 typedef typename GOS::Traits::Jacobian Matrix;
504
505 public:
506 NewtonPrepareStep(const GridOperator& go, TrialVector& u_)
507 : NewtonBase<GOS,TrlV,TstV>(go,u_)
508 , min_linear_reduction_(1e-3)
509 , fixed_linear_reduction_(0.0)
510 , reassemble_threshold_(0.0)
511 {}
512
513 NewtonPrepareStep(const GridOperator& go)
514 : NewtonBase<GOS,TrlV,TstV>(go)
515 , min_linear_reduction_(1e-3)
516 , fixed_linear_reduction_(0.0)
517 , reassemble_threshold_(0.0)
518 , hanging_node_modifications_(false)
519 {}
520
527 void setMinLinearReduction(RFType min_linear_reduction)
528 {
529 min_linear_reduction_ = min_linear_reduction;
530 }
531
536 void setFixedLinearReduction(bool fixed_linear_reduction)
537 {
538 fixed_linear_reduction_ = fixed_linear_reduction;
539 }
540
548 void setReassembleThreshold(RFType reassemble_threshold)
549 {
550 reassemble_threshold_ = reassemble_threshold;
551 }
552
553 void setHangingNodeModifications(bool b)
554 {
555 hanging_node_modifications_ = b;
556 }
557
558 virtual void prepare_step(Matrix& A, TstV& ) override
559 {
560 this->reassembled_ = false;
561 if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
562 {
563 if (hanging_node_modifications_){
564 auto dirichletValues = *this->u_;
565 // Set all non dirichlet values to zero
566 Dune::PDELab::set_shifted_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, dirichletValues);
567 // Set all constrained DOFs to zero in solution
568 Dune::PDELab::set_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, *this->u_);
569 // Copy correct Dirichlet values back into solution vector
570 Dune::PDELab::copy_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), dirichletValues, *this->u_);
571 // Interpolate periodic constraints / hanging nodes
572 this->gridoperator_.localAssembler().backtransform(*this->u_);
573 }
574
575 if (this->verbosity_level_ >= 3)
576 std::cout << " Reassembling matrix..." << std::endl;
577 A = 0.0;
578 this->gridoperator_.jacobian(*this->u_, A);
579 this->reassembled_ = true;
580 }
581
582 // corner case: force_iteration_==true. Use min_linear_reduction when initial defect is less than AbsoluteLimit.
583 if (fixed_linear_reduction_ == true || (this->res_.iterations==0 && this->res_.first_defect<this->abs_limit_))
584 this->linear_reduction_ = min_linear_reduction_;
585 else {
586 // determine maximum defect, where Newton is converged.
587 RFType stop_defect =
588 std::max(this->res_.first_defect * this->reduction_,
589 this->abs_limit_);
590
591 /*
592 To achieve second order convergence of newton
593 we need a linear reduction of at least
594 current_defect^2/prev_defect^2.
595 For the last newton step a linear reduction of
596 1/10*end_defect/current_defect
597 is sufficient for convergence.
598 */
599 this->linear_reduction_ =
600 std::max( stop_defect/(10*this->res_.defect),
601 std::min(min_linear_reduction_,this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_)) );
602 }
603
604 this->prev_defect_ = this->res_.defect;
605
606 ios_base_all_saver restorer(std::cout); // store old ios flags
607
608 if (this->verbosity_level_ >= 3)
609 std::cout << " requested linear reduction: "
610 << std::setw(12) << std::setprecision(4) << std::scientific
611 << this->linear_reduction_ << std::endl;
612 }
613
614 private:
615 RFType min_linear_reduction_;
616 bool fixed_linear_reduction_;
617 RFType reassemble_threshold_;
618 bool hanging_node_modifications_;
619 }; // end class NewtonPrepareStep
620
621 template<class GOS, class TrlV, class TstV>
622 class NewtonLineSearch : public virtual NewtonBase<GOS,TrlV,TstV>
623 {
624 typedef GOS GridOperator;
625 typedef TrlV TrialVector;
626 typedef TstV TestVector;
627
628 typedef typename TestVector::ElementType RFType;
629
630 public:
631 enum Strategy {
633 noLineSearch,
637 hackbuschReusken,
639 hackbuschReuskenAcceptBest };
640
641 NewtonLineSearch(const GridOperator& go, TrialVector& u_)
642 : NewtonBase<GOS,TrlV,TstV>(go,u_)
643 , strategy_(hackbuschReusken)
644 , maxit_(10)
645 , damping_factor_(0.5)
646 {}
647
648 NewtonLineSearch(const GridOperator& go)
649 : NewtonBase<GOS,TrlV,TstV>(go)
650 , strategy_(hackbuschReusken)
651 , maxit_(10)
652 , damping_factor_(0.5)
653 {}
654
655 void setLineSearchStrategy(Strategy strategy)
656 {
657 strategy_ = strategy;
658 }
659
660 void setLineSearchStrategy(std::string strategy)
661 {
662 strategy_ = strategyFromName(strategy);
663 }
664
665 void setLineSearchMaxIterations(unsigned int maxit)
666 {
667 maxit_ = maxit;
668 }
669
670 void setLineSearchDampingFactor(RFType damping_factor)
671 {
672 damping_factor_ = damping_factor;
673 }
674
675 virtual void line_search(TrialVector& z, TestVector& r) override
676 {
677 if ((strategy_ == noLineSearch) || (this->res_.defect < this->abs_limit_))
678 {
679 this->u_->axpy(-1.0, z);
680 this->defect(r);
681 return;
682 }
683
684 if (this->verbosity_level_ >= 4)
685 std::cout << " Performing line search..." << std::endl;
686 RFType lambda = 1.0;
687 RFType best_lambda = 0.0;
688 RFType best_defect = this->res_.defect;
689 TrialVector prev_u(*this->u_);
690 unsigned int i = 0;
691 ios_base_all_saver restorer(std::cout); // store old ios flags
692
693 while (1)
694 {
695 if (this->verbosity_level_ >= 4)
696 std::cout << " trying line search damping factor: "
697 << std::setw(12) << std::setprecision(4) << std::scientific
698 << lambda
699 << std::endl;
700
701 this->u_->axpy(-lambda, z);
702 try {
703 this->defect(r);
704 }
705 catch (NewtonDefectError&)
706 {
707 if (this->verbosity_level_ >= 4)
708 std::cout << " NaNs detected" << std::endl;
709 } // ignore NaNs and try again with lower lambda
710
711 if (this->res_.defect <= (1.0 - lambda/4) * this->prev_defect_)
712 {
713 if (this->verbosity_level_ >= 4)
714 std::cout << " line search converged" << std::endl;
715 break;
716 }
717
718 if (this->res_.defect < best_defect)
719 {
720 best_defect = this->res_.defect;
721 best_lambda = lambda;
722 }
723
724 if (++i >= maxit_)
725 {
726 if (this->verbosity_level_ >= 4)
727 std::cout << " max line search iterations exceeded" << std::endl;
728 switch (strategy_)
729 {
730 case hackbuschReusken:
731 *this->u_ = prev_u;
732 this->defect(r);
733 DUNE_THROW(NewtonLineSearchError,
734 "NewtonLineSearch::line_search(): line search failed, "
735 "max iteration count reached, "
736 "defect did not improve enough");
737 case hackbuschReuskenAcceptBest:
738 if (best_lambda == 0.0)
739 {
740 *this->u_ = prev_u;
741 this->defect(r);
742 DUNE_THROW(NewtonLineSearchError,
743 "NewtonLineSearch::line_search(): line search failed, "
744 "max iteration count reached, "
745 "defect did not improve in any of the iterations");
746 }
747 if (best_lambda != lambda)
748 {
749 *this->u_ = prev_u;
750 this->u_->axpy(-best_lambda, z);
751 this->defect(r);
752 }
753 break;
754 case noLineSearch:
755 break;
756 }
757 break;
758 }
759
760 lambda *= damping_factor_;
761 *this->u_ = prev_u;
762 }
763 if (this->verbosity_level_ >= 4)
764 std::cout << " line search damping factor: "
765 << std::setw(12) << std::setprecision(4) << std::scientific
766 << lambda << std::endl;
767 } // end line_search
768
769 protected:
771 Strategy strategyFromName(const std::string & s) {
772 if (s == "noLineSearch")
773 return noLineSearch;
774 if (s == "hackbuschReusken")
775 return hackbuschReusken;
776 if (s == "hackbuschReuskenAcceptBest")
777 return hackbuschReuskenAcceptBest;
778 DUNE_THROW(Exception, "unknown line search strategy" << s);
779 }
780
781 private:
782 Strategy strategy_;
783 unsigned int maxit_;
784 RFType damping_factor_;
785 }; // end class NewtonLineSearch
786
787 template<class GOS, class S, class TrlV, class TstV = TrlV>
788 class Newton : public NewtonSolver<GOS,S,TrlV,TstV>
789 , public NewtonTerminate<GOS,TrlV,TstV>
790 , public NewtonLineSearch<GOS,TrlV,TstV>
791 , public NewtonPrepareStep<GOS,TrlV,TstV>
792 {
793 typedef GOS GridOperator;
794 typedef S Solver;
795 typedef TrlV TrialVector;
796
797 public:
798 [[deprecated("This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")]]
799 Newton(const GridOperator& go, TrialVector& u_, Solver& solver_)
800 : NewtonBase<GOS,TrlV,TstV>(go,u_)
801 , NewtonSolver<GOS,S,TrlV,TstV>(go,u_,solver_)
802 , NewtonTerminate<GOS,TrlV,TstV>(go,u_)
803 , NewtonLineSearch<GOS,TrlV,TstV>(go,u_)
804 , NewtonPrepareStep<GOS,TrlV,TstV>(go,u_)
805 {}
806 [[deprecated("This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")]]
807 Newton(const GridOperator& go, Solver& solver_)
808 : NewtonBase<GOS,TrlV,TstV>(go)
809 , NewtonSolver<GOS,S,TrlV,TstV>(go,solver_)
810 , NewtonTerminate<GOS,TrlV,TstV>(go)
811 , NewtonLineSearch<GOS,TrlV,TstV>(go)
812 , NewtonPrepareStep<GOS,TrlV,TstV>(go)
813 {}
815
836 void setParameters(const Dune::ParameterTree & param)
837 {
838 typedef typename TstV::ElementType RFType;
839 if (param.hasKey("VerbosityLevel"))
840 this->setVerbosityLevel(
841 param.get<unsigned int>("VerbosityLevel"));
842 if (param.hasKey("Reduction"))
843 this->setReduction(
844 param.get<RFType>("Reduction"));
845 if (param.hasKey("MaxIterations"))
846 this->setMaxIterations(
847 param.get<unsigned int>("MaxIterations"));
848 if (param.hasKey("ForceIteration"))
849 this->setForceIteration(
850 param.get<bool>("ForceIteration"));
851 if (param.hasKey("AbsoluteLimit"))
852 this->setAbsoluteLimit(
853 param.get<RFType>("AbsoluteLimit"));
854 if (param.hasKey("MinLinearReduction"))
855 this->setMinLinearReduction(
856 param.get<RFType>("MinLinearReduction"));
857 if (param.hasKey("FixedLinearReduction"))
858 this->setFixedLinearReduction(
859 param.get<bool>("FixedLinearReduction"));
860 if (param.hasKey("ReassembleThreshold"))
861 this->setReassembleThreshold(
862 param.get<RFType>("ReassembleThreshold"));
863 if (param.hasKey("LineSearchStrategy"))
864 this->setLineSearchStrategy(
865 param.get<std::string>("LineSearchStrategy"));
866 if (param.hasKey("LineSearchMaxIterations"))
867 this->setLineSearchMaxIterations(
868 param.get<unsigned int>("LineSearchMaxIterations"));
869 if (param.hasKey("LineSearchDampingFactor"))
870 this->setLineSearchDampingFactor(
871 param.get<RFType>("LineSearchDampingFactor"));
872 if (param.hasKey("KeepMatrix"))
873 this->setKeepMatrix(
874 param.get<bool>("KeepMatrix"));
875 if (param.hasKey("UseMaxNorm"))
876 this->setUseMaxNorm(
877 param.get<bool>("UseMaxNorm"));
878 }
879 }; // end class Newton
880 } // end namespace PDELab
881} // end namespace Dune
882
883#endif // DOXYGEN
884
885#endif // DUNE_PDELAB_NEWTON_NEWTON_HH
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
bool hasKey(const std::string &key) const
test for key
Definition: parametertree.cc:51
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
Utility class for storing and resetting stream attributes.
Dune namespace.
Definition: alignedallocator.hh:13
A hierarchical structure of string parameters.
Standard Dune debug streams.
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)