Loading...
Searching...
No Matches
continuation.h
Go to the documentation of this file.
#ifndef _RHEOLEF_CONTINUATION_H
#define _RHEOLEF_CONTINUATION_H
//
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// =========================================================================
// AUTHOR: Pierre.Saramito@imag.fr
// DATE: 19 janvt 2018
namespace rheolef {
} // namespace rheolef
#include "rheolef/continuation_option.h"
#include "rheolef/continuation_step.h"
#include "rheolef/keller.h"
namespace rheolef { namespace details {
// ======================================================
// completion
// ======================================================
template<class Problem>
int
continuation_solve (
const Problem& F,
typename Problem::value_type& uh,
odiststream* p_err,
const continuation_option& opts)
{
typename Problem::float_type res = opts.tol;
size_t n_iter = opts.newton_max_iter;
int status = damped_newton(F, newton_identity_preconditioner(), uh, res, n_iter, p_err);
if (p_err) *p_err << std::endl << std::endl;
return ((status == 0 && res <= opts.tol) || sqr(res) <= opts.tol) ? 0 : 1;
}
// --------------------
// with mesh adaptation
// --------------------
template<class Problem>
void continuation_internal (
Problem& F,
typename Problem::value_type& uh,
odiststream* p_out,
odiststream* p_err,
const continuation_option& opts)
{
typedef typename Problem::value_type value_type;
typedef typename Problem::float_type float_type;
std::string name = F.parameter_name();
opts.check();
float_type delta_parameter = opts.ini_delta_parameter;
value_type duh_dparameter = F.direction(uh);
float_type duh_dparameter_sign = opts.ini_direction;
float_type duh_dparameter_norm = F.space_norm(duh_dparameter);
value_type uh_prev = uh;
size_t min_delta_parameter_successive_count = 0;
if (p_err) *p_err << "#[continuation] n parameter" << std::endl;
for (size_t n = 0; n < opts.max_iter; n++) {
float_type delta_parameter_prev = delta_parameter;
value_type uh_prev2 = uh_prev;
uh_prev = uh;
delta_parameter = step_adjust (F, n, delta_parameter_prev, uh_prev, duh_dparameter, duh_dparameter_sign, opts, p_err, uh);
using std::isnan;
if (isnan(delta_parameter)) {
if (p_err) *p_err << "#[continuation] stop, since the parameter step is not-a-number" << std::endl;
return;
}
size_t i_adapt = 0;
do { // delta_parameter loop
int status = 0;
value_type uh_old_mesh = uh;
value_type uh_prev_old_mesh = uh_prev;
value_type duh_dparameter_old_mesh = duh_dparameter;
for (i_adapt = 0; status == 0 && i_adapt <= opts.n_adapt; ++i_adapt) { // adapt loop
status = continuation_solve (F, uh, p_err, opts);
// check when going back (TODO: also when i_adapt > 0 : compare with guess as uh_prev)
if (opts.do_check_going_back && i_adapt == 0 && status == 0 && n >= 1) {
float_type cos_angle1 = F.space_dot (uh-uh_prev, uh_prev-uh_prev2);
if (cos_angle1 < 0) {
if (delta_parameter > opts.min_delta_parameter) {
if (p_err) *p_err << "#[continuation] check cos_angle1="<<cos_angle1<<" < 0 : treated as failed"<<std::endl;
status = 1;
} else {
if (p_err) *p_err << "#[continuation] check cos_angle1="<<cos_angle1<<" < 0 : reverse accepted since delta_"<<name<<"="<<delta_parameter << " is minimum" <<std::endl;
}
}
float_type cos_angle2 = duh_dparameter_sign*F.space_dot (uh-uh_prev, duh_dparameter);
if (cos_angle2 < 0) {
if (delta_parameter > opts.min_delta_parameter) {
if (p_err) *p_err << "#[continuation] check cos_angle2="<<cos_angle2<<" < 0 : treated as failed"<<std::endl;
status = 1;
} else {
if (p_err) *p_err << "#[continuation] check cos_angle2="<<cos_angle2<<" < 0 : reverse accepted since delta_"<<name<<"="<<delta_parameter << " is minimum" <<std::endl;
}
}
}
// check when going back, second part
if (i_adapt == 0 && status == 0) {
value_type new_duh_dparameter = F.direction(uh);
float_type new_duh_dparameter_norm = F.space_norm(new_duh_dparameter);
float_type cos_angle = F.space_dot (new_duh_dparameter, duh_dparameter)/(duh_dparameter_norm*new_duh_dparameter_norm);
if (opts.do_check_going_back) {
if (fabs(cos_angle) < 1 - opts.tol_cos_angle) {
if (delta_parameter > opts.min_delta_parameter) {
if (p_err) *p_err << "#[continuation] check cos(new_dir,dir)="<<cos_angle
<<" : too much curvature, treated as failed" << std::endl;
status = 1;
} else {
if (p_err) *p_err << "#[continuation] check cos(new_dir,dir)="<<cos_angle
<<" : too much curvature (HINT: jump to another branch ?), but accepted since delta_"<<name<<"="<<delta_parameter << " is minimum" <<std::endl;
}
} else {
if (cos_angle > 0) {
if (p_err) *p_err << "#[continuation] check cos(new_dir,dir)="<<cos_angle <<" : ok" <<std::endl;
} else {
if (p_err) *p_err << "#[continuation] check cos(new_dir,dir)="<<cos_angle <<" < 0 : HINT: cross a bifurcation point ?" <<std::endl;
}
}
}
if (status == 0) {
duh_dparameter = new_duh_dparameter;
duh_dparameter_norm = new_duh_dparameter_norm;
duh_dparameter_sign = (cos_angle >= 0 || !opts.do_check_going_back) ? duh_dparameter_sign : -duh_dparameter_sign;
}
}
if (status == 0 && i_adapt < opts.n_adapt) {
// prepare next adapt loop:
if (p_err) *p_err << "#[continuation] prepare i_adapt="<< i_adapt+1 << " <= n_adapt="<<opts.n_adapt<<" iteration..." << std::endl;
uh_old_mesh = uh;
uh_prev_old_mesh = uh_prev;
duh_dparameter_old_mesh = duh_dparameter;
F.adapt (uh, opts);
uh = F.reinterpolate (uh);
}
} // end adapt loop
if (i_adapt > 0) i_adapt--; // undo the last increment in the for(i_adapt) loop TODO: add while (status == 0)
if (status == 0) {
break; // success!
}
if (status != 0 && i_adapt > 0) {
i_adapt--;
if (p_err) *p_err << "#[continuation] failed, but the previous "<<i_adapt<<"th adaped mesh has a valid solution for this parameter" << std::endl;
// back on the previous (old) adapted mesh
uh = uh_old_mesh;
uh_prev = uh_prev_old_mesh;
duh_dparameter = duh_dparameter_old_mesh;
F.reset_geo (uh); // reset problem to the previous (old) adapted mesh
status = 0;
break; // relative success!
}
if (status != 0) {
if (delta_parameter == opts.min_delta_parameter) {
if (p_err) *p_err << "#[continuation] stop, since cannot decrease more the parameter step" << std::endl;
return;
}
// here, we have failed => uh is invalid
// go back to previous valid and redo prediction with smaller step
F.set_parameter (F.parameter() - delta_parameter);
delta_parameter = max (opts.theta_decr*delta_parameter, opts.min_delta_parameter);
F.set_parameter (F.parameter() + delta_parameter);
uh = uh_prev;
F.reset_geo (uh); // reset on the previous uh_prev mesh
if (opts.do_prediction) {
uh = uh + (duh_dparameter_sign*delta_parameter)*duh_dparameter;
}
if (p_err) *p_err << "#[continuation] solve failed: decreases delta_"<<name<<"="<<delta_parameter << std::endl;
}
} while (true); // end delta_parameter loop
if (p_err) *p_err << "[continuation] "<< n+1 << " " << F.parameter() << std::endl;
if (p_out) F.put (*p_out, uh);
if (F.stop(uh)) {
if (p_err) *p_err << "#[continuation] stop, from problem specific stopping criteria" << std::endl;
return;
}
if (delta_parameter == opts.min_delta_parameter) {
min_delta_parameter_successive_count++;
if (opts.min_delta_parameter_successive_count_max > 0 &&
min_delta_parameter_successive_count > opts.min_delta_parameter_successive_count_max) {
if (p_err) *p_err << "#[continuation] stop, since too much iteration with delta_"<<name<<"="<<delta_parameter << std::endl;
return;
}
} else {
min_delta_parameter_successive_count = 0;
}
if (i_adapt > 0) {
uh_prev = F.reinterpolate (uh_prev);
duh_dparameter = F.reinterpolate (duh_dparameter);
}
F.refresh (F.parameter(), uh, duh_dparameter);
}
}
} // namespace details
// ------------------------------------------------------
// algo selection: with or without mesh adapt
// ------------------------------------------------------
// [verbatim_continuation]
template<class Problem>
void continuation (
Problem& F,
typename Problem::value_type& uh,
odiststream* p_out,
odiststream* p_err,
const continuation_option& opts = continuation_option())
// [verbatim_continuation]
{
constexpr bool has_adapt_value = details::has_inherited_member_adapt<Problem>::value;
continuation_option opts1 = opts;
opts1.n_adapt = (has_adapt_value ? opts.n_adapt : 0); // force n_adapt=0 when !has_adapt
details::add_missing_continuation<Problem> G (F);
details::continuation_internal (G, uh, p_out, p_err, opts1);
F = G; // TODO: avoid this copy: how to use a reference to F in G ? and inherit the class interface ? move cstor ?
// requiered by gnutef/mossolov/mossolov_regularized_adapt.cc : get F.parameter() after continuation
}
template<class Problem>
void
continuation (
keller<Problem>& F,
typename keller<Problem>::value_type& uh,
odiststream* p_out,
odiststream* p_err,
const continuation_option& opts = continuation_option())
{
constexpr bool has_adapt_value = details::has_inherited_member_adapt<Problem>::value;
continuation_option opts1 = opts;
opts1.n_adapt = (has_adapt_value ? opts.n_adapt : 0); // force n_adapt=0 when !has_adapt
details::continuation_internal (F, uh, p_out, p_err, opts1);
}
} // namespace rheolef
#endif // _RHEOLEF_CONTINUATION_H
rheolef::details::add_missing_continuation
Definition newton_add_missing.h:275
rheolef::keller
see the continuation page for the full documentation
Definition keller.h:32
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::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::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
This file is part of Rheolef.
Definition compiler_eigen.h:39
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::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
see the continuation_option page for the full documentation
Definition continuation_option.h:77
rheolef::continuation_option::ini_delta_parameter
Float ini_delta_parameter
Definition continuation_option.h:85
rheolef::continuation_option::newton_max_iter
size_t newton_max_iter
Definition continuation_option.h:82
rheolef::continuation_option::n_adapt
size_t n_adapt
Definition continuation_option.h:93
rheolef::continuation_option::do_check_going_back
bool do_check_going_back
Definition continuation_option.h:92
rheolef::continuation_option::check
void check() const
Definition continuation_option.h:124
rheolef::continuation_option::do_prediction
bool do_prediction
Definition continuation_option.h:91
rheolef::continuation_option::tol
Float tol
Definition continuation_option.h:80
rheolef::continuation_option::tol_cos_angle
Float tol_cos_angle
Definition continuation_option.h:90
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::continuation_option::max_iter
size_t max_iter
Definition continuation_option.h:81
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
rheolef::newton_identity_preconditioner
Definition damped-newton-generic.h:77