Loading...
Searching...
No Matches
keller.h
Go to the documentation of this file.
#ifndef _RHEOLEF_KELLER_H
#define _RHEOLEF_KELLER_H
#include "rheolef/newton_add_missing.h"
#include "rheolef/keller_details.h"
namespace rheolef {
// switch depending on existing Problem::adapt() method
template<class Problem, class Sfinae = typename details::has_inherited_member_adapt<Problem>::type>
class keller {};
// -----------------------------
// class without mesh adaptation
// -----------------------------
template<class Problem>
class keller<Problem,std::false_type> {
public:
typedef Float float_type;
typedef details::pair_with_linear_algebra<float_type,typename Problem::value_type> value_type;
keller(const Problem& pb, std::string metric="orthogonal");
keller(const keller<Problem,std::false_type>&);
void reset(const Problem& pb, std::string metric="orthogonal");
float_type parameter() const { return s; }
std::string parameter_name() const { return "s"; }
Problem& get_problem() { return pb; }
void set_parameter(float_type s1) { s = s1; }
value_type initial (std::string restart="") const;
void refresh (float_type s, const value_type& xh, const value_type& d_xh_ds) const;
value_type direction (value_type& xh) const;
value_type derivative_versus_parameter (const value_type& xh) const;
value_type residue (const value_type& xh) const;
int solve (value_type& uh, const continuation_option& opts) const;
solver::determinant_type update_derivative (const value_type& xh) const;
value_type derivative_solve (const value_type& mrh) const;
value_type derivative_trans_mult (const value_type& mrh) const;
float_type space_norm (const value_type& xh) const;
float_type dual_space_norm (const value_type& mrh) const;
float_type space_dot (const value_type& xh, const value_type& yh) const;
odiststream& put (odiststream& os, const value_type& xh) const;
idiststream& get (idiststream& is, value_type& xh);
bool stop (const value_type& xh) const { return pb.stop(xh.second); }
// adapt (dummy):
void adapt (const value_type&, const adapt_option&) {}
void reset_geo (const value_type&) {}
value_type reinterpolate (const value_type& xh) { return xh; }
protected:
float_type s;
mutable details::add_missing_continuation<Problem> pb;
mutable csr<float_type> A1;
mutable solver sA1;
mutable typename Problem::value_type m_df_d_parameter;
mutable float_type s0;
mutable value_type xh0, d_xh0_ds;
mutable bool has_init, has_refresh, has_spherical;
mutable distributor u_ownership;
};
template<class Problem>
keller<Problem,std::false_type>::keller (const Problem& pb1, std::string metric)
: s(0), pb(pb1), A1(), sA1(), m_df_d_parameter(), s0(0), xh0(), d_xh0_ds(),
has_init(false), has_refresh(false), has_spherical(false), u_ownership()
{
reset (pb1, metric);
}
template<class Problem>
keller<Problem,std::false_type>::keller (const keller<Problem,std::false_type>& x)
: s(x.s), pb(x.pb), A1(x.A1), sA1(x.sA1), m_df_d_parameter(x.m_df_d_parameter),
s0(x.s0), xh0(x.xh0), d_xh0_ds(x.d_xh0_ds),
has_init(x.has_init), has_refresh(x.has_refresh), has_spherical(x.has_spherical),
u_ownership(x.u_ownership)
{
std::string metric = (!has_spherical) ? "orthogonal" : "spherical";
reset (pb, metric);
}
template<class Problem>
void
keller<Problem,std::false_type>::reset (const Problem& pb1, std::string metric) {
pb = pb1;
check_macro (metric == "orthogonal" || metric == "spherical", "invalid metric="<<metric);
has_spherical = (metric != "orthogonal");
}
template<class Problem>
typename
keller<Problem,std::false_type>::value_type
keller<Problem,std::false_type>::initial (std::string restart) const {
s0 = 0;
xh0.second = pb.initial (restart);
xh0.first = pb.parameter();
typename Problem::value_type d_uh0_ds = pb.direction(xh0.second);
float_type c = 1/sqrt(1 + sqr(pb.space_norm(d_uh0_ds)));
d_xh0_ds.first = c;
d_xh0_ds.second = c*d_uh0_ds;
has_init = true;
return xh0;
}
template<class Problem>
typename
keller<Problem,std::false_type>::value_type
keller<Problem,std::false_type>::residue (const value_type& xh) const {
pb.set_parameter (xh.first);
value_type mrh;
if (!has_spherical) {
mrh.first = space_dot(d_xh0_ds, xh-xh0) - (s-s0);
} else {
mrh.first = sqr(space_norm(xh-xh0)) - sqr(s-s0);
}
mrh.second = pb.residue (xh.second);
return mrh;
}
template<class Problem>
typename solver::determinant_type
keller<Problem,std::false_type>::update_derivative (const value_type& xh) const {
pb.set_parameter (xh.first);
m_df_d_parameter = pb.derivative_versus_parameter(xh.second);
form df_du = pb.derivative(xh.second);
u_ownership = df_du.uu().row_ownership();
vec<float_type> b2 = details::get_unknown (m_df_d_parameter, u_ownership);
vec<float_type> b1;
Float c;
if (!has_spherical) {
b1 = details::get_unknown (pb.massify(d_xh0_ds.second), u_ownership);
c = d_xh0_ds.first;
} else {
typename Problem::value_type duh = xh.second - xh0.second;
b1 = details::get_unknown (pb.massify(Float(2)*duh), u_ownership);
c = 2*(xh.first - xh0.first);
}
A1 = {{ df_du.uu(), b2 },
{ trans(b1), c }};
solver_option sopt;
sopt.compute_determinant = true;
sA1 = solver_basic<float_type> (A1, sopt);
return sA1.det();
}
template<class Problem>
typename
keller<Problem,std::false_type>::value_type
keller<Problem,std::false_type>::derivative_solve (const value_type& mrh) const {
vec<float_type> MB = { details::get_unknown(mrh.second,u_ownership),
mrh.first };
vec<float_type> X = sA1.solve(MB);
value_type delta_xh;
delta_xh.second = mrh.second; details::reset (delta_xh.second);
details::set_unknown (delta_xh.second, X);
size_t iproc0 = 0;
delta_xh.first = (size_t(X.comm().rank()) == iproc0) ? X[X.size()-1] : 0;
#ifdef _RHEOLEF_HAVE_MPI
mpi::broadcast (X.comm(), delta_xh.first, iproc0);
#endif // _RHEOLEF_HAVE_MPI
return delta_xh;
}
template<class Problem>
typename
keller<Problem,std::false_type>::value_type
keller<Problem,std::false_type>::derivative_trans_mult (const value_type& mrh) const {
vec<float_type> MRH = { details::get_unknown(mrh.second,u_ownership),
mrh.first };
vec<float_type> MGH = A1.trans_mult(MRH);
value_type mgh = mrh; // for std::valarray alloc
mgh.second = mrh.second; details::reset (mgh.second);
details::set_unknown (mgh.second, MGH);
size_t iproc0 = 0;
mgh.first = (size_t(MGH.comm().rank()) == iproc0) ? MGH[MGH.size()-1] : 0;
#ifdef _RHEOLEF_HAVE_MPI
mpi::broadcast (MGH.comm(), mgh.first, iproc0);
#endif // _RHEOLEF_HAVE_MPI
return mgh;
}
template<class Problem>
int
keller<Problem,std::false_type>::solve (value_type& uh, const continuation_option& opts) const {
float_type res = opts.tol;
size_t n_iter = opts.max_iter;
odiststream* p_err = &derr;
int status = damped_newton (*this, 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;
}
template<class Problem>
typename
keller<Problem,std::false_type>::value_type
keller<Problem,std::false_type>::direction (value_type& xh) const {
if (!has_refresh) return d_xh0_ds;
update_derivative(xh);
return -derivative_solve(derivative_versus_parameter(xh));
}
template<class Problem>
void
keller<Problem,std::false_type>::refresh (float_type s, const value_type& xh, const value_type& d_xh_ds) const {
has_refresh = true;
d_xh0_ds = d_xh_ds;
xh0 = xh;
s0 = s;
}
template<class Problem>
typename
keller<Problem,std::false_type>::value_type
keller<Problem,std::false_type>::derivative_versus_parameter (const value_type& xh) const {
value_type m_df_ds;
if (!has_spherical) {
m_df_ds.first = -1;
} else {
m_df_ds.first = -2*(s-s0);
}
m_df_ds.second = xh.second; details::reset(m_df_ds.second); // allocate+reset
return m_df_ds;
}
template<class Problem>
typename
keller<Problem,std::false_type>::float_type
keller<Problem,std::false_type>::space_dot (const value_type& xh, const value_type& yh) const {
return xh.first*yh.first + pb.space_dot(xh.second, yh.second);
}
template<class Problem>
typename
keller<Problem,std::false_type>::float_type
keller<Problem,std::false_type>::space_norm (const value_type& xh) const {
return sqrt (space_dot(xh,xh));
}
template<class Problem>
typename
keller<Problem,std::false_type>::float_type
keller<Problem,std::false_type>::dual_space_norm (const value_type& mrh) const {
return sqrt (sqr(mrh.first) + sqr(pb.dual_space_norm(mrh.second)));
}
template<class Problem>
odiststream&
keller<Problem,std::false_type>::put (odiststream& os, const value_type& xh) const {
static const bool trace = true;
float_type prev_param = pb.parameter();
pb.set_parameter(xh.first);
if (trace) {
static size_t n = 0;
if (n==0) derr << "#[keller] n s "<<pb.parameter_name()<<" |x|"<<std::endl;
derr << "[keller] " << n << " " << s << " " << xh.first << " "
<< pb.space_norm(xh.second) << std::endl;
n++;
}
os << std::endl << "#s " << s << std::endl;
pb.put(os,xh.second);
pb.set_parameter(prev_param);
return os;
}
template<class Problem>
idiststream&
keller<Problem,std::false_type>::get (idiststream& is, value_type& xh) {
is >> catchmark("s") >> s;
if (!is) return is;
pb.get (is, xh.second);
xh.first = pb.parameter();
if (!has_init) { initial(); }
return is;
}
// -----------------------------
// class with mesh adaptation
// -----------------------------
template<class Problem>
class keller<Problem,std::true_type>
: public keller<Problem,std::false_type>
{
public:
typedef keller<Problem,std::false_type> base;
typedef typename base::float_type float_type;
typedef typename base::value_type value_type;
keller (const Problem& pb, std::string metric="orthogonal") : base(pb,metric) {}
// adapt:
void adapt (const value_type& xh, const adapt_option& aopt);
void reset_geo (const value_type& xh);
value_type reinterpolate (const value_type& xh);
};
template<class Problem>
void
keller<Problem,std::true_type>::adapt (const value_type& xh, const adapt_option& aopt) {
base::pb.adapt (xh.second, aopt);
base::xh0.second = base::pb.reinterpolate (base::xh0.second);
base::d_xh0_ds.second = base::pb.reinterpolate (base::d_xh0_ds.second);
}
template<class Problem>
void
keller<Problem,std::true_type>::reset_geo (const value_type& xh) {
base::pb.reset_geo (xh.second);
base::xh0.second = base::pb.reinterpolate (base::xh0.second);
base::d_xh0_ds.second = base::pb.reinterpolate (base::d_xh0_ds.second);
}
template<class Problem>
typename
keller<Problem,std::true_type>::value_type
keller<Problem,std::true_type>::reinterpolate (const value_type& xh) {
typename keller<Problem>::value_type yh = xh;
yh.second = base::pb.reinterpolate (yh.second);
return yh;
}
} // namespace rheolef
#endif // _RHEOLEF_KELLER_H
Float
see the Float page for the full documentation
form
see the form page for the full documentation
rheolef::catchmark
see the catchmark page for the full documentation
Definition catchmark.h:67
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::details::add_missing_continuation
Definition newton_add_missing.h:275
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::idiststream
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
rheolef::keller< Problem, std::false_type >
Definition keller.h:38
rheolef::keller< Problem, std::false_type >::stop
bool stop(const value_type &xh) const
Definition keller.h:63
rheolef::keller< Problem, std::false_type >::float_type
Float float_type
Definition keller.h:40
rheolef::keller< Problem, std::false_type >::parameter_name
std::string parameter_name() const
Definition keller.h:46
rheolef::keller< Problem, std::false_type >::A1
csr< float_type > A1
Definition keller.h:71
rheolef::keller< Problem, std::false_type >::has_init
bool has_init
Definition keller.h:76
rheolef::keller< Problem, std::false_type >::get_problem
Problem & get_problem()
Definition keller.h:47
rheolef::keller< Problem, std::false_type >::pb
details::add_missing_continuation< Problem > pb
Definition keller.h:70
rheolef::keller< Problem, std::false_type >::m_df_d_parameter
Problem::value_type m_df_d_parameter
Definition keller.h:73
rheolef::keller< Problem, std::false_type >::parameter
float_type parameter() const
Definition keller.h:45
rheolef::keller< Problem, std::false_type >::reset_geo
void reset_geo(const value_type &)
Definition keller.h:66
rheolef::keller< Problem, std::false_type >::sA1
solver sA1
Definition keller.h:72
rheolef::keller< Problem, std::false_type >::s
float_type s
Definition keller.h:69
rheolef::keller< Problem, std::false_type >::s0
float_type s0
Definition keller.h:74
rheolef::keller< Problem, std::false_type >::value_type
details::pair_with_linear_algebra< float_type, typename Problem::value_type > value_type
Definition keller.h:41
rheolef::keller< Problem, std::false_type >::set_parameter
void set_parameter(float_type s1)
Definition keller.h:48
rheolef::keller< Problem, std::false_type >::u_ownership
distributor u_ownership
Definition keller.h:77
rheolef::keller< Problem, std::false_type >::d_xh0_ds
value_type d_xh0_ds
Definition keller.h:75
rheolef::keller< Problem, std::false_type >::reinterpolate
value_type reinterpolate(const value_type &xh)
Definition keller.h:67
rheolef::keller< Problem, std::false_type >::adapt
void adapt(const value_type &, const adapt_option &)
Definition keller.h:65
rheolef::keller< Problem, std::true_type >::value_type
base::value_type value_type
Definition keller.h:285
rheolef::keller< Problem, std::true_type >::keller
keller(const Problem &pb, std::string metric="orthogonal")
Definition keller.h:286
rheolef::keller< Problem, std::true_type >::float_type
base::float_type float_type
Definition keller.h:284
rheolef::keller< Problem, std::true_type >::base
keller< Problem, std::false_type > base
Definition keller.h:283
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::solver_basic
Definition solver.h:264
rheolef::solver_option
see the solver_option page for the full documentation
Definition solver_option.h:155
rheolef::solver_option::compute_determinant
bool compute_determinant
Definition solver_option.h:169
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
solver
see the solver page for the full documentation
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::details::reset
void reset(T &x)
Definition keller_details.h:31
rheolef::details::get_unknown
size_t get_unknown(const T &x)
Definition keller_details.h:65
rheolef::details::set_unknown
void set_unknown(T &x, const T &value)
Definition keller_details.h:119
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::adapt
geo_basic< T, M > adapt(const field_basic< T, M > &uh, const adapt_option &opts)
adapt(uh,opts): see the adapt page for the full documentation
Definition adapt.cc:179
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::solve
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
Definition tiny_lu.h:92
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
std
STL namespace.
residue
field residue(Float p, const field &uh)
Definition p_laplacian_post.cc:35
rheolef::adapt_option
adapt_option: see the adapt page for the full documentation
Definition adapt.h:147
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::max_iter
size_t max_iter
Definition continuation_option.h:81
rheolef::details::pair_with_linear_algebra
Definition pair_with_linear_algebra.h:27