Rheolef  7.1
an efficient C++ finite element environment
newton-backtrack.h
Go to the documentation of this file.
1 # ifndef _RHEO_NEWTON_BACKTRACK_H
2 # define _RHEO_NEWTON_BACKTRACK_H
3 namespace rheolef {
24 
25 template <class Problem, class Preconditioner, class Field, class Real>
27  const Problem& P, const Preconditioner& T,
28  const Field& u_old, Float Tu_old, Field& delta_u, Real slope, Real norm_delta_u_max,
29  Field& u, Field& Fu, Real& Tu, Real& lambda, odiststream *p_derr = 0)
30 {
31  const Float alpha = 1e-4; // 1e-8 when strongly nonlinear
32  const Float eps_mach = std::numeric_limits<Float>::epsilon();
33  Float norm_delta_u = P.space_norm(delta_u);
34  if (norm_delta_u > norm_delta_u_max) {
35  Float c = norm_delta_u_max/norm_delta_u;
36  if (p_derr) *p_derr << "# damped-Newton/backtrack: warning: delta_u bounded by factor " << c << std::endl << std::flush;
37  delta_u = c*delta_u;
38  }
39  // compute lambda_min
40  Float norm_u = P.space_norm(u);
41  if (norm_u < norm_delta_u) norm_u = norm_delta_u;
42  Float lambda_min = eps_mach*norm_u/norm_delta_u;
43  if (lambda_min > 1) { // machine precision problem detected
44  u = u_old;
45  return 1;
46  }
47  lambda = std::max (Real(0.0), std::min (Real(1.0), lambda));
48  Float lambda_prev = 0;
49  Float Tu_prev = 0;
50  Float Tu_prev_old = 0;
51  for (size_t k = 0; true; k++) {
52  u = u_old + lambda*delta_u;
53  Fu = P.residue(u);
54  Tu = T(P,Fu);
55  Float lambda_next;
56  if (lambda < lambda_min) { // machine precision problem detected
57  u = u_old;
58  return 1;
59  } else if (Tu <= Tu_old + alpha*lambda*slope) {
60  return 0; // have a valid lambda
61  } else if (lambda == 1) {
62  // first iteration: first order recursion
63  lambda_next = - 0.5*slope/(Tu - Tu_old - slope);
64  } else {
65  // second and more iterations: second order recursion
66  Float z = Tu - Tu_old - lambda*slope;
67  Float z_prev = Tu_prev - Tu_prev_old - lambda_prev*slope;
68  Float a = ( z/sqr(lambda) - z_prev/sqr(lambda_prev))/(lambda - lambda_prev);
69  Float b = (- z*lambda_prev/sqr(lambda) + z_prev*lambda/sqr(lambda_prev))
70  /(lambda - lambda_prev);
71  if (a == 0) {
72  lambda_next = - slope/(2*b);
73  } else {
74  Float Delta = sqr(b) - 3*a*slope;
75  if (Delta < 0) {
76  if (p_derr) *p_derr << "# damped-Newton/backtrack: warning: machine precision reached" << std::endl << std::flush;
77  return 1;
78  }
79  lambda_next = (-b + sqrt(Delta))/(3*a);
80  }
81  lambda_next = std::min (lambda/2, lambda_next);
82  }
83  lambda_prev = lambda;
84  Tu_prev = Tu;
85  Tu_prev_old = Tu_old;
86  lambda = std::max (lambda/10, lambda_next);
87  }
88  return 0;
89 }
90 }// namespace rheolef
91 # endif // _RHEO_NEWTON_BACKTRACK_H
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
a
Definition: diffusion_isotropic.h:25
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::newton_backtrack
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)
Definition: newton-backtrack.h:26
u
Definition: leveque.h:25
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
Float
see the Float page for the full documentation
epsilon
Float epsilon
Definition: transmission_error.cc:25
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
T
Expr1::float_type T
Definition: field_expr.h:261
lambda
Definition: yield_slip_circle.h:34