3 #ifndef DUNE_PDELAB_NEWTON_NEWTON_HH
4 #define DUNE_PDELAB_NEWTON_NEWTON_HH
11 #include <type_traits>
15 #include <dune/common/stdstreams.hh>
16 #include <dune/common/exceptions.hh>
17 #include <dune/common/ios_state.hh>
18 #include <dune/common/timer.hh>
19 #include <dune/common/parametertree.hh>
30 template<
typename T1,
typename =
void>
36 struct HasSetReuse<T, decltype(std::declval<T>().setReuse(true), void())>
41 inline void setLinearSystemReuse(T& solver_backend,
bool reuse, std::true_type)
43 if (!solver_backend.getReuse() && reuse)
44 dwarn <<
"WARNING: Newton needed to override your choice to reuse the linear system in order to work!" << std::endl;
45 solver_backend.setReuse(reuse);
49 inline void setLinearSystemReuse(T&,
bool, std::false_type)
53 inline void setLinearSystemReuse(T& solver_backend,
bool reuse)
55 setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
67 template<
class RFType>
81 template<
class GOS,
class TrlV,
class TstV>
84 typedef GOS GridOperator;
85 typedef TrlV TrialVector;
86 typedef TstV TestVector;
88 typedef typename TestVector::ElementType RFType;
89 typedef typename GOS::Traits::Jacobian Matrix;
98 if (
gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
126 std::shared_ptr<TrialVector>
z_;
127 std::shared_ptr<TestVector>
r_;
128 std::shared_ptr<Matrix>
A_;
144 if (
gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
154 if (
gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
161 virtual void prepare_step(Matrix& A, TestVector& r) = 0;
162 virtual void line_search(TrialVector& z, TestVector& r) = 0;
163 virtual void defect(TestVector& r) = 0;
166 template<
class GOS,
class S,
class TrlV,
class TstV>
170 typedef GOS GridOperator;
171 typedef TrlV TrialVector;
172 typedef TstV TestVector;
174 typedef typename TestVector::ElementType RFType;
175 typedef typename GOS::Traits::Jacobian Matrix;
183 , result_valid_(false)
189 , result_valid_(false)
200 "NewtonSolver::result() called before NewtonSolver::solve()");
205 virtual void defect(TestVector& r)
override
212 "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
217 void linearSolve(Matrix& A, TrialVector& z, TestVector& r)
const
220 std::cout <<
" Solving linear system..." << std::endl;
223 Impl::setLinearSystemReuse(this->solver_, !this->
reassembled_);
226 ios_base_all_saver restorer(std::cout);
228 if (!this->solver_.result().converged)
230 "NewtonSolver::linearSolve(): Linear solver did not converge "
231 "in " << this->solver_.result().iterations <<
" iterations");
233 std::cout <<
" linear solver iterations: "
234 << std::setw(12) << solver_.result().iterations << std::endl
235 <<
" linear defect reduction: "
236 << std::setw(12) << std::setprecision(4) << std::scientific
237 << solver_.result().reduction << std::endl;
244 template<
class GOS,
class S,
class TrlV,
class TstV>
251 template<
class GOS,
class S,
class TrlV,
class TstV>
254 this->res_.iterations = 0;
255 this->res_.converged =
false;
256 this->res_.reduction = 1.0;
257 this->res_.conv_rate = 1.0;
258 this->res_.elapsed = 0.0;
259 this->res_.assembler_time = 0.0;
260 this->res_.linear_solver_time = 0.0;
261 this->res_.linear_solver_iterations = 0;
262 result_valid_ =
true;
269 this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
276 this->defect(*this->r_);
277 this->res_.first_defect = this->res_.defect;
278 this->prev_defect_ = this->res_.defect;
280 if (this->verbosity_level_ >= 2)
283 ios_base_all_saver restorer(std::cout);
284 std::cout <<
" Initial defect: "
285 << std::setw(12) << std::setprecision(4) << std::scientific
286 << this->res_.defect << std::endl;
291 this->A_ = std::make_shared<Matrix>(this->gridoperator_);
295 this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
298 while (!this->terminate())
300 if (this->verbosity_level_ >= 3)
301 std::cout <<
" Newton iteration " << this->res_.iterations
302 <<
" --------------------------------" << std::endl;
304 Timer assembler_timer;
313 this->prepare_step(*this->A_,*this->r_);
317 this->res_.assembler_time += assembler_timer.elapsed();
320 double assembler_time = assembler_timer.elapsed();
321 this->res_.assembler_time += assembler_time;
322 if (this->verbosity_level_ >= 3)
323 std::cout <<
" matrix assembly time: "
324 << std::setw(12) << std::setprecision(4) << std::scientific
325 << assembler_time << std::endl;
327 Timer linear_solver_timer;
334 this->linearSolve(*this->A_, *this->z_, *this->r_);
338 this->res_.linear_solver_time += linear_solver_timer.elapsed();
339 this->res_.linear_solver_iterations += this->solver_.result().iterations;
342 double linear_solver_time = linear_solver_timer.elapsed();
343 this->res_.linear_solver_time += linear_solver_time;
344 this->res_.linear_solver_iterations += this->solver_.result().iterations;
350 this->line_search(*this->z_, *this->r_);
354 if (this->reassembled_)
356 if (this->verbosity_level_ >= 3)
357 std::cout <<
" line search failed - trying again with reassembled matrix" << std::endl;
361 this->res_.reduction = this->res_.defect/this->res_.first_defect;
362 this->res_.iterations++;
363 this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
366 ios_base_all_saver restorer(std::cout);
368 if (this->verbosity_level_ >= 3)
369 std::cout <<
" linear solver time: "
370 << std::setw(12) << std::setprecision(4) << std::scientific
371 << linear_solver_time << std::endl
372 <<
" defect reduction (this iteration):"
373 << std::setw(12) << std::setprecision(4) << std::scientific
374 << this->res_.defect/this->prev_defect_ << std::endl
375 <<
" defect reduction (total): "
376 << std::setw(12) << std::setprecision(4) << std::scientific
377 << this->res_.reduction << std::endl
379 << std::setw(12) << std::setprecision(4) << std::scientific
380 << this->res_.defect << std::endl;
381 if (this->verbosity_level_ == 2)
382 std::cout <<
" Newton iteration " << std::setw(2) << this->res_.iterations
384 << std::setw(12) << std::setprecision(4) << std::scientific
386 <<
". Reduction (this): "
387 << std::setw(12) << std::setprecision(4) << std::scientific
388 << this->res_.defect/this->prev_defect_
389 <<
". Reduction (total): "
390 << std::setw(12) << std::setprecision(4) << std::scientific
391 << this->res_.reduction << std::endl;
396 this->res_.elapsed = timer.elapsed();
399 this->res_.elapsed = timer.elapsed();
401 ios_base_all_saver restorer(std::cout);
403 if (this->verbosity_level_ == 1)
404 std::cout <<
" Newton converged after " << std::setw(2) << this->res_.iterations
405 <<
" iterations. Reduction: "
406 << std::setw(12) << std::setprecision(4) << std::scientific
407 << this->res_.reduction
408 <<
" (" << std::setprecision(4) << this->res_.elapsed <<
"s)"
411 if(!this->keep_matrix_)
415 template<
class GOS,
class TrlV,
class TstV>
418 typedef GOS GridOperator;
419 typedef TrlV TrialVector;
421 typedef typename TstV::ElementType RFType;
427 , force_iteration_(false)
436 , force_iteration_(false)
454 force_iteration_ = force_iteration;
470 "NewtonTerminate::terminate(): Maximum iteration count reached");
476 bool force_iteration_;
479 template<
class GOS,
class TrlV,
class TstV>
482 typedef GOS GridOperator;
483 typedef TrlV TrialVector;
485 typedef typename TstV::ElementType RFType;
486 typedef typename GOS::Traits::Jacobian Matrix;
491 , min_linear_reduction_(1
e-3)
492 , fixed_linear_reduction_(0.0)
493 , reassemble_threshold_(0.0)
498 , min_linear_reduction_(1
e-3)
499 , fixed_linear_reduction_(0.0)
500 , reassemble_threshold_(0.0)
511 min_linear_reduction_ = min_linear_reduction;
520 fixed_linear_reduction_ = fixed_linear_reduction;
532 reassemble_threshold_ = reassemble_threshold;
538 if (this->
res_.
defect/this->prev_defect_ > reassemble_threshold_)
541 std::cout <<
" Reassembling matrix..." << std::endl;
547 if (fixed_linear_reduction_ ==
true)
564 this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_) )
569 std::min(min_linear_reduction_,this->
res_.
defect*this->res_.defect/(this->prev_defect_*this->prev_defect_));
574 ios_base_all_saver restorer(std::cout);
577 std::cout <<
" requested linear reduction: "
578 << std::setw(12) << std::setprecision(4) << std::scientific
583 RFType min_linear_reduction_;
584 bool fixed_linear_reduction_;
585 RFType reassemble_threshold_;
588 template<
class GOS,
class TrlV,
class TstV>
591 typedef GOS GridOperator;
592 typedef TrlV TrialVector;
593 typedef TstV TestVector;
595 typedef typename TestVector::ElementType RFType;
612 , damping_factor_(0.5)
619 , damping_factor_(0.5)
624 strategy_ = strategy;
639 damping_factor_ = damping_factor;
646 this->
u_->axpy(-1.0, z);
652 std::cout <<
" Performing line search..." << std::endl;
654 RFType best_lambda = 0.0;
656 TrialVector prev_u(*this->
u_);
658 ios_base_all_saver restorer(std::cout);
663 std::cout <<
" trying line search damping factor: "
664 << std::setw(12) << std::setprecision(4) << std::scientific
668 this->
u_->axpy(-lambda, z);
675 std::cout <<
" NaNs detected" << std::endl;
678 if (this->
res_.
defect <= (1.0 - lambda/4) * this->prev_defect_)
681 std::cout <<
" line search converged" << std::endl;
688 best_lambda = lambda;
694 std::cout <<
" max line search iterations exceeded" << std::endl;
701 "NewtonLineSearch::line_search(): line search failed, "
702 "max iteration count reached, "
703 "defect did not improve enough");
705 if (best_lambda == 0.0)
710 "NewtonLineSearch::line_search(): line search failed, "
711 "max iteration count reached, "
712 "defect did not improve in any of the iterations");
714 if (best_lambda != lambda)
717 this->
u_->axpy(-best_lambda, z);
727 lambda *= damping_factor_;
731 std::cout <<
" line search damping factor: "
732 << std::setw(12) << std::setprecision(4) << std::scientific
733 << lambda << std::endl;
739 if (
s ==
"noLineSearch")
741 if (
s ==
"hackbuschReusken")
743 if (
s ==
"hackbuschReuskenAcceptBest")
745 DUNE_THROW(
Exception,
"unknown line search strategy" <<
s);
751 RFType damping_factor_;
754 template<
class GOS,
class S,
class TrlV,
class TstV = TrlV>
760 typedef GOS GridOperator;
762 typedef TrlV TrialVector;
765 Newton(
const GridOperator& go, TrialVector&
u_, Solver& solver_)
772 Newton(
const GridOperator& go, Solver& solver_)
803 typedef typename TstV::ElementType RFType;
804 if (param.hasKey(
"VerbosityLevel"))
806 param.get<
unsigned int>(
"VerbosityLevel"));
807 if (param.hasKey(
"Reduction"))
809 param.get<RFType>(
"Reduction"));
810 if (param.hasKey(
"MaxIterations"))
812 param.get<
unsigned int>(
"MaxIterations"));
813 if (param.hasKey(
"ForceIteration"))
815 param.get<
bool>(
"ForceIteration"));
816 if (param.hasKey(
"AbsoluteLimit"))
818 param.get<RFType>(
"AbsoluteLimit"));
819 if (param.hasKey(
"MinLinearReduction"))
821 param.get<RFType>(
"MinLinearReduction"));
822 if (param.hasKey(
"FixedLinearReduction"))
824 param.get<
bool>(
"FixedLinearReduction"));
825 if (param.hasKey(
"ReassembleThreshold"))
827 param.get<RFType>(
"ReassembleThreshold"));
828 if (param.hasKey(
"LineSearchStrategy"))
830 param.get<std::string>(
"LineSearchStrategy"));
831 if (param.hasKey(
"LineSearchMaxIterations"))
833 param.get<
unsigned int>(
"LineSearchMaxIterations"));
834 if (param.hasKey(
"LineSearchDampingFactor"))
836 param.get<RFType>(
"LineSearchDampingFactor"));
837 if (param.hasKey(
"KeepMatrix"))
839 param.get<
bool>(
"KeepMatrix"));
845 #endif // DUNE_PDELAB_NEWTON_NEWTON_HH