Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
26namespace rheolef {
113} // namespace rheolef
114
115#include "rheolef/continuation_option.h"
116#include "rheolef/continuation_step.h"
117#include "rheolef/keller.h"
118
119namespace rheolef { namespace details {
120
121// ======================================================
122// completion
123// ======================================================
124template<class Problem>
125int
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.newton_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// --------------------
142template<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]
302template<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}
319template<class Problem>
320void
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
see the continuation page for the full documentation
Definition keller.h:32
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
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