1 #ifndef _RHEOLEF_CONTINUATION_STEP_H
2 #define _RHEOLEF_CONTINUATION_STEP_H
24 namespace rheolef {
namespace details {
26 template<
class Solver>
27 typename Solver::float_type
31 typename Solver::float_type delta_parameter_prev,
32 const typename Solver::value_type& uh_prev,
33 const typename Solver::value_type& duh_dparameter,
34 const typename Solver::float_type& duh_dparameter_sign,
37 typename Solver::value_type& uh_guess)
39 typedef typename Solver::value_type
value_type;
40 typedef typename Solver::float_type
float_type;
41 std::string
name = F.parameter_name();
42 float_type delta_parameter = delta_parameter_prev;
43 if (p_err) *p_err <<
"#[continuation] delta_"<<
name<<
"(0) = " << delta_parameter << std::endl;
44 for (
size_t k = 1, k_max = 10;
true; ++k) {
45 F.set_parameter (F.parameter() + delta_parameter);
48 uh0 = uh0 + (duh_dparameter_sign*delta_parameter)*duh_dparameter;
50 F.update_derivative (uh0);
52 value_type delta_uh0 = - F.derivative_solve (m_r0h);
54 F.update_derivative (uh1);
56 value_type delta_uh1 = - F.derivative_solve (m_r1h);
62 value_type uh_possible_guess = (r2 < r1) ? uh2 : uh1;
64 if (sqr(r0) < opts.
tol && sqr(r1) < opts.
tol) {
67 float_type new_delta_parameter = delta_parameter*theta;
70 delta_parameter = new_delta_parameter;
71 uh_guess = uh_possible_guess;
72 if (p_err) *p_err <<
"#[continuation] prediction: too small residues (r0="<<r0<<
", r1="<<r1<<
"): increase delta_"<<
name<<
"("<<k<<
") = " << delta_parameter<<std::endl;
73 return delta_parameter;
81 if (p_err) *p_err <<
"#[continuation] prediction: not-a-number in residues (r0="<<r0<<
", r1="<<r1<<
"): set theta="<<theta<<std::endl;
85 if (p_err) *p_err <<
"#[continuation] prediction: optimal rate (kappa="<<kappa0<<
", theta="<<theta<<
") with delta_"<<
name<<
"("<<k<<
") = " << delta_parameter<<std::endl;
86 uh_guess = uh_possible_guess;
87 return delta_parameter;
92 float_type new_delta_parameter = delta_parameter*theta;
97 new_delta_parameter = min (new_delta_parameter, opts.
theta_incr*delta_parameter_prev);
101 if (p_err) *p_err <<
"#[continuation] prediction: avoid infinite loop with delta_"<<
name<<
"("<<k<<
") = " << delta_parameter<<std::endl;
102 uh_guess = uh_possible_guess;
103 return delta_parameter;
105 F.set_parameter (F.parameter() - delta_parameter);
106 delta_parameter = new_delta_parameter;
107 if (p_err) *p_err <<
"#[continuation] prediction: loop with delta_"<<
name<<
"("<<k<<
") = " << new_delta_parameter<<std::endl;
112 #endif // _RHEOLEF_CONTINUATION_STEP_H