1 # ifndef _RHEO_DAMPED_NEWTON_GENERIC_H
2 # define _RHEO_DAMPED_NEWTON_GENERIC_H
3 #include "rheolef/dis_cpu_time.h"
24 #include "rheolef/newton-backtrack.h"
28 template <
class Problem,
class Preconditioner,
class Field,
class Real,
class Size>
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.));
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);
42 if (p_derr) *p_derr <<
n <<
" " << P.dual_space_norm(Fu) <<
" " << sqrt(2.*Tu)
44 if (2.*Tu <= sqr(tol) ||
n >= max_iter) {
49 Real slope =
T.slope(P, Fu, delta_u);
54 P,
T, u_old, Tu_old, delta_u, slope, norm_delta_u_max,
u, Fu, Tu,
lambda, p_derr);
60 Field Gu =
T.grad(P,Fu);
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;
68 if (p_derr) *p_derr <<
"# damped-Newton: warning: gradient is zero up to machine precision" << std::endl << std::flush;
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));
82 template <
class Problem>
83 Float operator() (
const Problem& P,
const typename Problem::value_type& MFu,
const typename Problem::value_type& delta_u)
const {
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);
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));
97 # endif // _RHEO_DAMPED_NEWTON_GENERIC_H