Loading...
Searching...
No Matches
continuation_step.h
Go to the documentation of this file.
#ifndef _RHEOLEF_CONTINUATION_STEP_H
#define _RHEOLEF_CONTINUATION_STEP_H
namespace rheolef { namespace details {
template<class Solver>
typename 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)
{
typedef typename Solver::value_type value_type;
typedef typename Solver::float_type float_type;
std::string name = F.parameter_name();
float_type delta_parameter = delta_parameter_prev;
if (p_err) *p_err << "#[continuation] delta_"<<name<<"(0) = " << delta_parameter << std::endl;
for (size_t k = 1, k_max = 10; true; ++k) {
F.set_parameter (F.parameter() + delta_parameter);
value_type uh0 = uh_prev;
if (opts.do_prediction) {
uh0 = uh0 + (duh_dparameter_sign*delta_parameter)*duh_dparameter;
}
F.update_derivative (uh0);
value_type m_r0h = F.residue (uh0);
value_type delta_uh0 = - F.derivative_solve (m_r0h);
value_type uh1 = uh0 + delta_uh0;
F.update_derivative (uh1);
value_type m_r1h = F.residue (uh1);
value_type delta_uh1 = - F.derivative_solve (m_r1h);
value_type uh2 = uh1 + delta_uh1;
value_type m_r2h = F.residue (uh2);
float_type r0 = F.dual_space_norm (m_r0h);
float_type r1 = F.dual_space_norm (m_r1h);
float_type r2 = F.dual_space_norm (m_r2h);
value_type uh_possible_guess = (r2 < r1) ? uh2 : uh1;
r0 = max(r0, std::numeric_limits<float_type>::epsilon());
if (sqr(r0) < opts.tol && sqr(r1) < opts.tol) {
// r0 and r1 are close to machine precision: step is probably very small: increase it strongly
float_type theta = std::max (opts.theta_incr, 1/opts.theta_decr);
float_type new_delta_parameter = delta_parameter*theta;
new_delta_parameter = min (new_delta_parameter, opts.max_delta_parameter);
new_delta_parameter = max (new_delta_parameter, opts.min_delta_parameter);
delta_parameter = new_delta_parameter;
uh_guess = uh_possible_guess;
if (p_err) *p_err << "#[continuation] prediction: too small residues (r0="<<r0<<", r1="<<r1<<"): increase delta_"<<name<<"("<<k<<") = " << delta_parameter<<std::endl;
return delta_parameter;
}
float_type kappa0 = r1/r0;
float_type theta = sqrt(opts.kappa/kappa0);
using std::isnan;
if (isnan(theta)) {
// solver generates not-a-number => decrease delta_s
theta = opts.theta_decr;
if (p_err) *p_err << "#[continuation] prediction: not-a-number in residues (r0="<<r0<<", r1="<<r1<<"): set theta="<<theta<<std::endl;
}
if (kappa0 < 1 && fabs(1 - theta) < opts.theta_variation) {
// optimal convergence with kappa rate
if (p_err) *p_err << "#[continuation] prediction: optimal rate (kappa="<<kappa0<<", theta="<<theta<<") with delta_"<<name<<"("<<k<<") = " << delta_parameter<<std::endl;
uh_guess = uh_possible_guess;
return delta_parameter;
}
// otherwise: converge too fast: increase delta_parameter
// or diverge: decrease delta_parameter
theta = max(theta, opts.theta_decr);
float_type new_delta_parameter = delta_parameter*theta;
new_delta_parameter = min (new_delta_parameter, opts.max_delta_parameter);
new_delta_parameter = max (new_delta_parameter, opts.min_delta_parameter);
if (kappa0 < 1) {
// when converge, avoid to increase too fast
new_delta_parameter = min (new_delta_parameter, opts.theta_incr*delta_parameter_prev);
}
if (k >= k_max || fabs(new_delta_parameter - delta_parameter) < std::numeric_limits<float_type>::epsilon()) {
// avoid infinite loop
if (p_err) *p_err << "#[continuation] prediction: avoid infinite loop with delta_"<<name<<"("<<k<<") = " << delta_parameter<<std::endl;
uh_guess = uh_possible_guess;
return delta_parameter;
}
F.set_parameter (F.parameter() - delta_parameter);
delta_parameter = new_delta_parameter;
if (p_err) *p_err << "#[continuation] prediction: loop with delta_"<<name<<"("<<k<<") = " << new_delta_parameter<<std::endl;
}
}
}} // namespace rheolef::details
#endif // _RHEOLEF_CONTINUATION_STEP_H
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
rheolef::details::step_adjust
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)
Definition continuation_step.h:28
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::continuation_option
see the continuation_option page for the full documentation
Definition continuation_option.h:77
rheolef::continuation_option::kappa
Float kappa
Definition continuation_option.h:79
rheolef::continuation_option::theta_incr
Float theta_incr
Definition continuation_option.h:87
rheolef::continuation_option::do_prediction
bool do_prediction
Definition continuation_option.h:91
rheolef::continuation_option::theta_variation
Float theta_variation
Definition continuation_option.h:88
rheolef::continuation_option::tol
Float tol
Definition continuation_option.h:80
rheolef::continuation_option::max_delta_parameter
Float max_delta_parameter
Definition continuation_option.h:84
rheolef::continuation_option::theta_decr
Float theta_decr
Definition continuation_option.h:86
rheolef::continuation_option::min_delta_parameter
Float min_delta_parameter
Definition continuation_option.h:83