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();
155 value_type duh_dparameter = F.direction(uh);
157 float_type duh_dparameter_norm = F.space_norm(duh_dparameter);
158 value_type uh_prev = uh;
159 size_t min_delta_parameter_successive_count = 0;
160 if (p_err) *p_err <<
"#[continuation] n parameter" << std::endl;
161 for (
size_t n = 0; n < opts.
max_iter; n++) {
162 float_type delta_parameter_prev = delta_parameter;
163 value_type uh_prev2 = uh_prev;
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;
174 value_type uh_old_mesh = uh;
175 value_type uh_prev_old_mesh = uh_prev;
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;
229 if (status == 0 && i_adapt < opts.
n_adapt) {
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);
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)