Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
newton-backtrack.h
Go to the documentation of this file.
1# ifndef _RHEO_NEWTON_BACKTRACK_H
2# define _RHEO_NEWTON_BACKTRACK_H
23namespace rheolef {
24
25template <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
see the Float page for the full documentation
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
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 leveque.h:25