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