Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
continuation_step.h
Go to the documentation of this file.
1#ifndef _RHEOLEF_CONTINUATION_STEP_H
2#define _RHEOLEF_CONTINUATION_STEP_H
23
24namespace rheolef { namespace details {
25
26template<class Solver>
27typename Solver::float_type
29 Solver& F,
30 size_t n,
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,
35 const continuation_option& opts,
36 odiststream* p_err,
37 typename Solver::value_type& uh_guess)
38{
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);
46 value_type uh0 = uh_prev;
47 if (opts.do_prediction) {
48 uh0 = uh0 + (duh_dparameter_sign*delta_parameter)*duh_dparameter;
49 }
50 F.update_derivative (uh0);
51 value_type m_r0h = F.residue (uh0);
52 value_type delta_uh0 = - F.derivative_solve (m_r0h);
53 value_type uh1 = uh0 + delta_uh0;
54 F.update_derivative (uh1);
55 value_type m_r1h = F.residue (uh1);
56 value_type delta_uh1 = - F.derivative_solve (m_r1h);
57 value_type uh2 = uh1 + delta_uh1;
58 value_type m_r2h = F.residue (uh2);
59 float_type r0 = F.dual_space_norm (m_r0h);
60 float_type r1 = F.dual_space_norm (m_r1h);
61 float_type r2 = F.dual_space_norm (m_r2h);
62 value_type uh_possible_guess = (r2 < r1) ? uh2 : uh1;
63 r0 = max(r0, std::numeric_limits<float_type>::epsilon());
64 if (sqr(r0) < opts.tol && sqr(r1) < opts.tol) {
65 // r0 and r1 are close to machine precision: step is probably very small: increase it strongly
66 float_type theta = std::max (opts.theta_incr, 1/opts.theta_decr);
67 float_type new_delta_parameter = delta_parameter*theta;
68 new_delta_parameter = min (new_delta_parameter, opts.max_delta_parameter);
69 new_delta_parameter = max (new_delta_parameter, opts.min_delta_parameter);
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;
74 }
75 float_type kappa0 = r1/r0;
76 float_type theta = sqrt(opts.kappa/kappa0);
77 using std::isnan;
78 if (isnan(theta)) {
79 // solver generates not-a-number => decrease delta_s
80 theta = opts.theta_decr;
81 if (p_err) *p_err << "#[continuation] prediction: not-a-number in residues (r0="<<r0<<", r1="<<r1<<"): set theta="<<theta<<std::endl;
82 }
83 if (kappa0 < 1 && fabs(1 - theta) < opts.theta_variation) {
84 // optimal convergence with kappa rate
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;
88 }
89 // otherwise: converge too fast: increase delta_parameter
90 // or diverge: decrease delta_parameter
91 theta = max(theta, opts.theta_decr);
92 float_type new_delta_parameter = delta_parameter*theta;
93 new_delta_parameter = min (new_delta_parameter, opts.max_delta_parameter);
94 new_delta_parameter = max (new_delta_parameter, opts.min_delta_parameter);
95 if (kappa0 < 1) {
96 // when converge, avoid to increase too fast
97 new_delta_parameter = min (new_delta_parameter, opts.theta_incr*delta_parameter_prev);
98 }
99 if (k >= k_max || fabs(new_delta_parameter - delta_parameter) < std::numeric_limits<float_type>::epsilon()) {
100 // avoid infinite loop
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;
104 }
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;
108 }
109}
110
111}} // namespace rheolef::details
112#endif // _RHEOLEF_CONTINUATION_STEP_H
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
Solver::float_type step_adjust(Solver &F, size_t n, typename Solver::float_type delta_parameter_prev, const typename Solver::value_type &uh_prev, const typename Solver::value_type &duh_dparameter, const typename Solver::float_type &duh_dparameter_sign, const continuation_option &opts, odiststream *p_err, typename Solver::value_type &uh_guess)
This file is part of Rheolef.
see the continuation_option page for the full documentation