Rheolef  7.1
an efficient C++ finite element environment
continuation.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_CONTINUATION_H
2 #define _RHEOLEF_CONTINUATION_H
3 //
4 // This file is part of Rheolef.
5 //
6 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7 //
8 // Rheolef is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // Rheolef is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with Rheolef; if not, write to the Free Software
20 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // =========================================================================
23 // AUTHOR: Pierre.Saramito@imag.fr
24 // DATE: 19 janvt 2018
25 
26 namespace rheolef {
113 } // namespace rheolef
114 
115 #include "rheolef/continuation_option.h"
116 #include "rheolef/continuation_step.h"
117 #include "rheolef/keller.h"
118 
119 namespace rheolef { namespace details {
120 
121 // ======================================================
122 // completion
123 // ======================================================
124 template<class Problem>
125 int
127  const Problem& F,
128  typename Problem::value_type& uh,
129  odiststream* p_err,
130  const continuation_option& opts)
131 {
132  typename Problem::float_type res = opts.tol;
133  size_t n_iter = opts.max_iter;
134  int status = damped_newton(F, newton_identity_preconditioner(), uh, res, n_iter, p_err);
135  if (p_err) *p_err << std::endl << std::endl;
136  return ((status == 0 && res <= opts.tol) || sqr(res) <= opts.tol) ? 0 : 1;
137 }
138 
139 // --------------------
140 // with mesh adaptation
141 // --------------------
142 template<class Problem>
144  Problem& F,
145  typename Problem::value_type& uh,
146  odiststream* p_out,
147  odiststream* p_err,
148  const continuation_option& opts)
149 {
150  typedef typename Problem::value_type value_type;
151  typedef typename Problem::float_type float_type;
152  std::string name = F.parameter_name();
153  opts.check();
154  float_type delta_parameter = opts.ini_delta_parameter;
155  value_type duh_dparameter = F.direction(uh);
156  float_type duh_dparameter_sign = opts.ini_direction;
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;
164  uh_prev = uh;
165  delta_parameter = step_adjust (F, n, delta_parameter_prev, uh_prev, duh_dparameter, duh_dparameter_sign, opts, p_err, uh);
166  using std::isnan;
167  if (isnan(delta_parameter)) {
168  if (p_err) *p_err << "#[continuation] stop, since the parameter step is not-a-number" << std::endl;
169  return;
170  }
171  size_t i_adapt = 0;
172  do { // delta_parameter loop
173  int status = 0;
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) { // adapt loop
178  status = continuation_solve (F, uh, p_err, opts);
179  // check when going back (TODO: also when i_adapt > 0 : compare with guess as uh_prev)
180  if (opts.do_check_going_back && i_adapt == 0 && status == 0 && n >= 1) {
181  float_type cos_angle1 = F.space_dot (uh-uh_prev, uh_prev-uh_prev2);
182  if (cos_angle1 < 0) {
183  if (delta_parameter > opts.min_delta_parameter) {
184  if (p_err) *p_err << "#[continuation] check cos_angle1="<<cos_angle1<<" < 0 : treated as failed"<<std::endl;
185  status = 1;
186  } else {
187  if (p_err) *p_err << "#[continuation] check cos_angle1="<<cos_angle1<<" < 0 : reverse accepted since delta_"<<name<<"="<<delta_parameter << " is minimum" <<std::endl;
188  }
189  }
190  float_type cos_angle2 = duh_dparameter_sign*F.space_dot (uh-uh_prev, duh_dparameter);
191  if (cos_angle2 < 0) {
192  if (delta_parameter > opts.min_delta_parameter) {
193  if (p_err) *p_err << "#[continuation] check cos_angle2="<<cos_angle2<<" < 0 : treated as failed"<<std::endl;
194  status = 1;
195  } else {
196  if (p_err) *p_err << "#[continuation] check cos_angle2="<<cos_angle2<<" < 0 : reverse accepted since delta_"<<name<<"="<<delta_parameter << " is minimum" <<std::endl;
197  }
198  }
199  }
200  // check when going back, second part
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);
205  if (opts.do_check_going_back) {
206  if (fabs(cos_angle) < 1 - opts.tol_cos_angle) {
207  if (delta_parameter > opts.min_delta_parameter) {
208  if (p_err) *p_err << "#[continuation] check cos(new_dir,dir)="<<cos_angle
209  <<" : too much curvature, treated as failed" << std::endl;
210  status = 1;
211  } else {
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;
214  }
215  } else {
216  if (cos_angle > 0) {
217  if (p_err) *p_err << "#[continuation] check cos(new_dir,dir)="<<cos_angle <<" : ok" <<std::endl;
218  } else {
219  if (p_err) *p_err << "#[continuation] check cos(new_dir,dir)="<<cos_angle <<" < 0 : HINT: cross a bifurcation point ?" <<std::endl;
220  }
221  }
222  }
223  if (status == 0) {
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;
227  }
228  }
229  if (status == 0 && i_adapt < opts.n_adapt) {
230  // prepare next adapt loop:
231  if (p_err) *p_err << "#[continuation] prepare i_adapt="<< i_adapt+1 << " <= n_adapt="<<opts.n_adapt<<" iteration..." << std::endl;
232  uh_old_mesh = uh;
233  uh_prev_old_mesh = uh_prev;
234  duh_dparameter_old_mesh = duh_dparameter;
235  F.adapt (uh, opts);
236  uh = F.reinterpolate (uh);
237  }
238  } // end adapt loop
239  if (i_adapt > 0) i_adapt--; // undo the last increment in the for(i_adapt) loop TODO: add while (status == 0)
240  if (status == 0) {
241  break; // success!
242  }
243  if (status != 0 && i_adapt > 0) {
244  i_adapt--;
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;
246  // back on the previous (old) adapted mesh
247  uh = uh_old_mesh;
248  uh_prev = uh_prev_old_mesh;
249  duh_dparameter = duh_dparameter_old_mesh;
250  F.reset_geo (uh); // reset problem to the previous (old) adapted mesh
251  status = 0;
252  break; // relative success!
253  }
254  if (status != 0) {
255  if (delta_parameter == opts.min_delta_parameter) {
256  if (p_err) *p_err << "#[continuation] stop, since cannot decrease more the parameter step" << std::endl;
257  return;
258  }
259  // here, we have failed => uh is invalid
260  // go back to previous valid and redo prediction with smaller step
261  F.set_parameter (F.parameter() - delta_parameter);
262  delta_parameter = max (opts.theta_decr*delta_parameter, opts.min_delta_parameter);
263  F.set_parameter (F.parameter() + delta_parameter);
264  uh = uh_prev;
265  F.reset_geo (uh); // reset on the previous uh_prev mesh
266  if (opts.do_prediction) {
267  uh = uh + (duh_dparameter_sign*delta_parameter)*duh_dparameter;
268  }
269  if (p_err) *p_err << "#[continuation] solve failed: decreases delta_"<<name<<"="<<delta_parameter << std::endl;
270  }
271  } while (true); // end delta_parameter loop
272  if (p_err) *p_err << "[continuation] "<< n+1 << " " << F.parameter() << std::endl;
273  if (p_out) F.put (*p_out, uh);
274  if (F.stop(uh)) {
275  if (p_err) *p_err << "#[continuation] stop, from problem specific stopping criteria" << std::endl;
276  return;
277  }
278  if (delta_parameter == opts.min_delta_parameter) {
279  min_delta_parameter_successive_count++;
281  min_delta_parameter_successive_count > opts.min_delta_parameter_successive_count_max) {
282  if (p_err) *p_err << "#[continuation] stop, since too much iteration with delta_"<<name<<"="<<delta_parameter << std::endl;
283  return;
284  }
285  } else {
286  min_delta_parameter_successive_count = 0;
287  }
288  if (i_adapt > 0) {
289  uh_prev = F.reinterpolate (uh_prev);
290  duh_dparameter = F.reinterpolate (duh_dparameter);
291  }
292  F.refresh (F.parameter(), uh, duh_dparameter);
293  }
294 }
295 
296 } // namespace details
297 // ------------------------------------------------------
298 // algo selection: with or without mesh adapt
299 // ------------------------------------------------------
300 // [verbatim_continuation]
302 template<class Problem>
304  Problem& F,
305  typename Problem::value_type& uh,
306  odiststream* p_out,
307  odiststream* p_err,
309 // [verbatim_continuation]
310 {
311  constexpr bool has_adapt_value = details::has_inherited_member_adapt<Problem>::value;
312  continuation_option opts1 = opts;
313  opts1.n_adapt = (has_adapt_value ? opts.n_adapt : 0); // force n_adapt=0 when !has_adapt
315  details::continuation_internal (G, uh, p_out, p_err, opts1);
316  F = G; // TODO: avoid this copy: how to use a reference to F in G ? and inherit the class interface ? move cstor ?
317  // requiered by gnutef/mossolov/mossolov_regularized_adapt.cc : get F.parameter() after continuation
318 }
319 template<class Problem>
320 void
322  keller<Problem>& F,
323  typename keller<Problem>::value_type& uh,
324  odiststream* p_out,
325  odiststream* p_err,
327 {
328  constexpr bool has_adapt_value = details::has_inherited_member_adapt<Problem>::value;
329  continuation_option opts1 = opts;
330  opts1.n_adapt = (has_adapt_value ? opts.n_adapt : 0); // force n_adapt=0 when !has_adapt
331  details::continuation_internal (F, uh, p_out, p_err, opts1);
332 }
333 
334 } // namespace rheolef
335 #endif // _RHEOLEF_CONTINUATION_H
rheolef::continuation_option::theta_decr
Float theta_decr
Definition: continuation_option.h:86
rheolef::continuation_option::ini_delta_parameter
Float ini_delta_parameter
Definition: continuation_option.h:85
rheolef::details::add_missing_continuation
Definition: newton_add_missing.h:269
rheolef::continuation_option
see the continuation_option page for the full documentation
Definition: continuation_option.h:77
rheolef::continuation_option::tol
Float tol
Definition: continuation_option.h:80
rheolef::continuation_option::ini_direction
Float ini_direction
Definition: continuation_option.h:78
rheolef::continuation_option::min_delta_parameter_successive_count_max
size_t min_delta_parameter_successive_count_max
Definition: continuation_option.h:89
rheolef::float_type
typename float_traits< value_type >::type float_type
Definition: field_expr_recursive.h:501
rheolef::continuation_option::check
void check() const
Definition: continuation_option.h:124
rheolef::details::continuation_solve
int continuation_solve(const Problem &F, typename Problem::value_type &uh, odiststream *p_err, const continuation_option &opts)
Definition: continuation.h:126
rheolef::continuation_option::min_delta_parameter
Float min_delta_parameter
Definition: continuation_option.h:83
rheolef::continuation_option::max_iter
size_t max_iter
Definition: continuation_option.h:81
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::value_type
result_type value_type
Definition: field_expr_recursive.h:499
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::details::continuation_internal
void continuation_internal(Problem &F, typename Problem::value_type &uh, odiststream *p_out, odiststream *p_err, const continuation_option &opts)
Definition: continuation.h:143
rheolef::continuation
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
Definition: continuation.h:303
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
rheolef::damped_newton
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
Definition: damped-newton-generic.h:29
rheolef::continuation_option::n_adapt
size_t n_adapt
Definition: continuation_option.h:93
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::continuation_option::do_prediction
bool do_prediction
Definition: continuation_option.h:91
rheolef::newton_identity_preconditioner
Definition: damped-newton-generic.h:77
rheolef::continuation_option::tol_cos_angle
Float tol_cos_angle
Definition: continuation_option.h:90
rheolef::keller
see the continuation page for the full documentation
Definition: keller.h:32
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
rheolef::continuation_option::do_check_going_back
bool do_check_going_back
Definition: continuation_option.h:92