Rheolef  7.2
an efficient C++ finite element environment
continuation_step.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_CONTINUATION_STEP_H
2 #define _RHEOLEF_CONTINUATION_STEP_H
23 
24 namespace rheolef { namespace details {
25 
26 template<class Solver>
27 typename 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;
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
typename float_traits< value_type >::type float_type
result_type value_type
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
Float epsilon