1#ifndef _RHEOLEF_CONTINUATION_H
2#define _RHEOLEF_CONTINUATION_H
115#include "rheolef/continuation_option.h"
116#include "rheolef/continuation_step.h"
117#include "rheolef/keller.h"
119namespace rheolef {
namespace details {
124template<
class Problem>
128 typename Problem::value_type& uh,
132 typename Problem::float_type res = opts.
tol;
135 if (p_err) *p_err << std::endl << std::endl;
136 return ((
status == 0 && res <= opts.
tol) || sqr(res) <= opts.
tol) ? 0 : 1;
142template<
class Problem>
145 typename Problem::value_type& uh,
150 typedef typename Problem::value_type
value_type;
151 typedef typename Problem::float_type
float_type;
152 std::string
name = F.parameter_name();
157 float_type duh_dparameter_norm = F.space_norm(duh_dparameter);
159 size_t min_delta_parameter_successive_count = 0;
160 if (p_err) *p_err <<
"#[continuation] n parameter" << std::endl;
162 float_type delta_parameter_prev = delta_parameter;
165 delta_parameter =
step_adjust (F,
n, delta_parameter_prev, uh_prev, duh_dparameter, duh_dparameter_sign, opts, p_err, uh);
167 if (isnan(delta_parameter)) {
168 if (p_err) *p_err <<
"#[continuation] stop, since the parameter step is not-a-number" << std::endl;
176 value_type duh_dparameter_old_mesh = duh_dparameter;
177 for (i_adapt = 0;
status == 0 && i_adapt <= opts.
n_adapt; ++i_adapt) {
181 float_type cos_angle1 = F.space_dot (uh-uh_prev, uh_prev-uh_prev2);
182 if (cos_angle1 < 0) {
184 if (p_err) *p_err <<
"#[continuation] check cos_angle1="<<cos_angle1<<
" < 0 : treated as failed"<<std::endl;
187 if (p_err) *p_err <<
"#[continuation] check cos_angle1="<<cos_angle1<<
" < 0 : reverse accepted since delta_"<<
name<<
"="<<delta_parameter <<
" is minimum" <<std::endl;
190 float_type cos_angle2 = duh_dparameter_sign*F.space_dot (uh-uh_prev, duh_dparameter);
191 if (cos_angle2 < 0) {
193 if (p_err) *p_err <<
"#[continuation] check cos_angle2="<<cos_angle2<<
" < 0 : treated as failed"<<std::endl;
196 if (p_err) *p_err <<
"#[continuation] check cos_angle2="<<cos_angle2<<
" < 0 : reverse accepted since delta_"<<
name<<
"="<<delta_parameter <<
" is minimum" <<std::endl;
201 if (i_adapt == 0 &&
status == 0) {
202 value_type new_duh_dparameter = F.direction(uh);
203 float_type new_duh_dparameter_norm = F.space_norm(new_duh_dparameter);
204 float_type cos_angle = F.space_dot (new_duh_dparameter, duh_dparameter)/(duh_dparameter_norm*new_duh_dparameter_norm);
208 if (p_err) *p_err <<
"#[continuation] check cos(new_dir,dir)="<<cos_angle
209 <<
" : too much curvature, treated as failed" << std::endl;
212 if (p_err) *p_err <<
"#[continuation] check cos(new_dir,dir)="<<cos_angle
213 <<
" : too much curvature (HINT: jump to another branch ?), but accepted since delta_"<<
name<<
"="<<delta_parameter <<
" is minimum" <<std::endl;
217 if (p_err) *p_err <<
"#[continuation] check cos(new_dir,dir)="<<cos_angle <<
" : ok" <<std::endl;
219 if (p_err) *p_err <<
"#[continuation] check cos(new_dir,dir)="<<cos_angle <<
" < 0 : HINT: cross a bifurcation point ?" <<std::endl;
224 duh_dparameter = new_duh_dparameter;
225 duh_dparameter_norm = new_duh_dparameter_norm;
226 duh_dparameter_sign = (cos_angle >= 0 || !opts.
do_check_going_back) ? duh_dparameter_sign : -duh_dparameter_sign;
231 if (p_err) *p_err <<
"#[continuation] prepare i_adapt="<< i_adapt+1 <<
" <= n_adapt="<<opts.
n_adapt<<
" iteration..." << std::endl;
233 uh_prev_old_mesh = uh_prev;
234 duh_dparameter_old_mesh = duh_dparameter;
236 uh = F.reinterpolate (uh);
239 if (i_adapt > 0) i_adapt--;
243 if (
status != 0 && i_adapt > 0) {
245 if (p_err) *p_err <<
"#[continuation] failed, but the previous "<<i_adapt<<
"th adaped mesh has a valid solution for this parameter" << std::endl;
248 uh_prev = uh_prev_old_mesh;
249 duh_dparameter = duh_dparameter_old_mesh;
256 if (p_err) *p_err <<
"#[continuation] stop, since cannot decrease more the parameter step" << std::endl;
261 F.set_parameter (F.parameter() - delta_parameter);
263 F.set_parameter (F.parameter() + delta_parameter);
267 uh = uh + (duh_dparameter_sign*delta_parameter)*duh_dparameter;
269 if (p_err) *p_err <<
"#[continuation] solve failed: decreases delta_"<<
name<<
"="<<delta_parameter << std::endl;
272 if (p_err) *p_err <<
"[continuation] "<<
n+1 <<
" " << F.parameter() << std::endl;
273 if (p_out) F.put (*p_out, uh);
275 if (p_err) *p_err <<
"#[continuation] stop, from problem specific stopping criteria" << std::endl;
279 min_delta_parameter_successive_count++;
282 if (p_err) *p_err <<
"#[continuation] stop, since too much iteration with delta_"<<
name<<
"="<<delta_parameter << std::endl;
286 min_delta_parameter_successive_count = 0;
289 uh_prev = F.reinterpolate (uh_prev);
290 duh_dparameter = F.reinterpolate (duh_dparameter);
292 F.refresh (F.parameter(), uh, duh_dparameter);
302template<
class Problem>
305 typename Problem::value_type& uh,
311 constexpr bool has_adapt_value = details::has_inherited_member_adapt<Problem>::value;
313 opts1.
n_adapt = (has_adapt_value ? opts.n_adapt : 0);
319template<
class Problem>
328 constexpr bool has_adapt_value = details::has_inherited_member_adapt<Problem>::value;
330 opts1.
n_adapt = (has_adapt_value ? opts.n_adapt : 0);
see the continuation page for the full documentation
odiststream: see the diststream page for the full documentation
typename float_traits< value_type >::type float_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)
void continuation_internal(Problem &F, typename Problem::value_type &uh, odiststream *p_out, odiststream *p_err, const continuation_option &opts)
int continuation_solve(const Problem &F, typename Problem::value_type &uh, odiststream *p_err, const continuation_option &opts)
This file is part of Rheolef.
void continuation(Problem &F, typename Problem::value_type &uh, odiststream *p_out, odiststream *p_err, const continuation_option &opts=continuation_option())
see the continuation page for the full documentation
int damped_newton(const Problem &P, const Preconditioner &T, Field &u, Real &tol, Size &max_iter, odiststream *p_derr=0)
see the damped_newton page for the full documentation
see the continuation_option page for the full documentation
Float ini_delta_parameter
size_t min_delta_parameter_successive_count_max
Float min_delta_parameter