Rheolef  7.1
an efficient C++ finite element environment
damped-newton-generic.h
Go to the documentation of this file.
1 # ifndef _RHEO_DAMPED_NEWTON_GENERIC_H
2 # define _RHEO_DAMPED_NEWTON_GENERIC_H
23 #include "rheolef/dis_cpu_time.h"
24 #include "rheolef/newton-backtrack.h"
25 namespace rheolef {
26 
28 template <class Problem, class Preconditioner, class Field, class Real, class Size>
29 int damped_newton (const Problem& P, const Preconditioner& T, Field& u, Real& tol, Size& max_iter, odiststream *p_derr = 0) {
30  const Real delta_u_max_factor = 100;
31  Real norm_delta_u_max = delta_u_max_factor*std::max(P.space_norm(u), Real(1.));
32  Real lambda = 1;
33  Real Tu = 0;
34  double t0 = dis_wall_time();
35  double c0 = dis_cpu_time();
36  if (p_derr) *p_derr << "# damped-Newton: n r T lambda wall-time cpu" << std::endl << std::flush;
37  for (Size n = 0; true; n++) {
38  P.update_derivative (u);
39  Field Fu = P.residue(u);
40  Field delta_u = -P.derivative_solve (Fu);
41  Tu = T(P,Fu,delta_u);
42  if (p_derr) *p_derr << n << " " << P.dual_space_norm(Fu) << " " << sqrt(2.*Tu)
43  << " " << lambda << " " << dis_wall_time()-t0 << " " << dis_cpu_time()-c0 << std::endl << std::flush;
44  if (2.*Tu <= sqr(tol) || n >= max_iter) {
45  max_iter = n;
46  tol = sqrt(2.*Tu);
47  return 0;
48  }
49  Real slope = T.slope(P, Fu, delta_u);
50  Field u_old = u;
51  Real Tu_old = Tu;
52  lambda = 1;
53  int status = newton_backtrack (
54  P, T, u_old, Tu_old, delta_u, slope, norm_delta_u_max, u, Fu, Tu, lambda, p_derr);
55  if (status != 0) {
56  // check if grad(T)(u) approx 0 ?
57  // let T(u) = 0.5*|A*F(u)|_M
58  // compute grad(T) = F(u)^T A^T M A F(u)
59  // and compare to F(u) : |grad(T)(u)| / |F(u)| > tol ?
60  Field Gu = T.grad(P,Fu);
61  const Float eps_mach = std::numeric_limits<Float>::epsilon();
62  max_iter = n;
63  tol = sqrt(2.*Tu);
64  if (P.space_norm(Gu) > eps_mach*P.dual_space_norm(Fu)) {
65  if (p_derr) *p_derr << "# damped-Newton: warning: machine precision reached" << std::endl << std::flush;
66  return 0;
67  } else {
68  if (p_derr) *p_derr << "# damped-Newton: warning: gradient is zero up to machine precision" << std::endl << std::flush;
69  return 1;
70  }
71  }
72  }
73  tol = sqrt(2*Tu);
74  return 1;
75 }
76 
78  template <class Problem>
79  Float operator() (const Problem& P, const typename Problem::value_type& MFv) const {
80  return 0.5*sqr(P.dual_space_norm(MFv));
81  }
82  template <class Problem>
83  Float operator() (const Problem& P, const typename Problem::value_type& MFu, const typename Problem::value_type& delta_u) const {
84  return operator() (P, MFu);
85  }
86  template <class Problem>
87  typename Problem::value_type grad (const Problem& P, const typename Problem::value_type& MFu) const {
88  return P.derivative_trans_mult (MFu);
89  }
90  template <class Problem>
91  Float slope (const Problem& P, const typename Problem::value_type& MFu, const typename Problem::value_type& delta_u) const {
92  return -sqr(P.dual_space_norm(MFu));
93  }
94 };
95 
96 }// namespace rheolef
97 # endif // _RHEO_DAMPED_NEWTON_GENERIC_H
see the Float page for the full documentation
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
Expr1::float_type T
Definition: field_expr.h:261
This file is part of Rheolef.
double dis_wall_time()
double dis_cpu_time()
Definition: dis_cpu_time.cc:45
int damped_newton(const Problem &P, const Preconditioner &T, Field &u, Real &tol, Size &max_iter, odiststream *p_derr=0)
see the damped_newton page for the full documentation
int newton_backtrack(const Problem &P, const Preconditioner &T, const Field &u_old, Float Tu_old, Field &delta_u, Real slope, Real norm_delta_u_max, Field &u, Field &Fu, Real &Tu, Real &lambda, odiststream *p_derr=0)
Float operator()(const Problem &P, const typename Problem::value_type &MFv) const
Float slope(const Problem &P, const typename Problem::value_type &MFu, const typename Problem::value_type &delta_u) const
Problem::value_type grad(const Problem &P, const typename Problem::value_type &MFu) const
Definition: leveque.h:25
Float u(const point &x)
Float epsilon