newton-backtrack.h
Go to the documentation of this file.
# ifndef _RHEO_NEWTON_BACKTRACK_H
# define _RHEO_NEWTON_BACKTRACK_H
namespace rheolef {
template <class Problem, class Preconditioner, class Field, class Real>
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)
{
const Float alpha = 1e-4; // 1e-8 when strongly nonlinear
const Float eps_mach = std::numeric_limits<Float>::epsilon();
Float norm_delta_u = P.space_norm(delta_u);
if (norm_delta_u > norm_delta_u_max) {
Float c = norm_delta_u_max/norm_delta_u;
if (p_derr) *p_derr << "# damped-Newton/backtrack: warning: delta_u bounded by factor " << c << std::endl << std::flush;
delta_u = c*delta_u;
}
// compute lambda_min
Float norm_u = P.space_norm(u);
if (norm_u < norm_delta_u) norm_u = norm_delta_u;
Float lambda_min = eps_mach*norm_u/norm_delta_u;
if (lambda_min > 1) { // machine precision problem detected
u = u_old;
return 1;
}
lambda = std::max (Real(0.0), std::min (Real(1.0), lambda));
Float lambda_prev = 0;
Float Tu_prev = 0;
Float Tu_prev_old = 0;
for (size_t k = 0; true; k++) {
u = u_old + lambda*delta_u;
Fu = P.residue(u);
Tu = T(P,Fu);
Float lambda_next;
if (lambda < lambda_min) { // machine precision problem detected
u = u_old;
return 1;
} else if (Tu <= Tu_old + alpha*lambda*slope) {
return 0; // have a valid lambda
} else if (lambda == 1) {
// first iteration: first order recursion
lambda_next = - 0.5*slope/(Tu - Tu_old - slope);
} else {
// second and more iterations: second order recursion
Float z = Tu - Tu_old - lambda*slope;
Float z_prev = Tu_prev - Tu_prev_old - lambda_prev*slope;
Float a = ( z/sqr(lambda) - z_prev/sqr(lambda_prev))/(lambda - lambda_prev);
Float b = (- z*lambda_prev/sqr(lambda) + z_prev*lambda/sqr(lambda_prev))
/(lambda - lambda_prev);
if (a == 0) {
lambda_next = - slope/(2*b);
} else {
Float Delta = sqr(b) - 3*a*slope;
if (Delta < 0) {
if (p_derr) *p_derr << "# damped-Newton/backtrack: warning: machine precision reached" << std::endl << std::flush;
return 1;
}
lambda_next = (-b + sqrt(Delta))/(3*a);
}
lambda_next = std::min (lambda/2, lambda_next);
}
lambda_prev = lambda;
Tu_prev = Tu;
Tu_prev_old = Tu_old;
lambda = std::max (lambda/10, lambda_next);
}
return 0;
}
}// namespace rheolef
# endif // _RHEO_NEWTON_BACKTRACK_H
Float
see the Float page for the full documentation
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
T
Expr1::float_type T
Definition: field_expr.h:261
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
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
a
Definition: diffusion_isotropic.h:25
lambda
Definition: yield_slip_circle.h:34
u
Definition: leveque.h:25
epsilon
Float epsilon
Definition: transmission_error.cc:25