1 #ifndef _RHEOLEF_KELLER_H
2 #define _RHEOLEF_KELLER_H
24 #include "rheolef/newton_add_missing.h"
25 #include "rheolef/keller_details.h"
31 template<class Problem, class Sfinae = typename details::has_inherited_member_adapt<Problem>::type>
37 template<
class Problem>
42 keller(
const Problem& pb, std::string metric=
"orthogonal");
44 void reset(
const Problem& pb, std::string metric=
"orthogonal");
49 value_type initial (std::string restart=
"")
const;
55 solver::determinant_type update_derivative (
const value_type& xh)
const;
62 idiststream& get (idiststream& is,
value_type& xh);
79 template<
class Problem>
81 : s(0), pb(pb1), A1(), sA1(), m_df_d_parameter(), s0(0), xh0(), d_xh0_ds(),
82 has_init(false), has_refresh(false), has_spherical(false), u_ownership()
86 template<
class Problem>
88 : s(x.s), pb(x.pb), A1(x.A1), sA1(x.sA1), m_df_d_parameter(x.m_df_d_parameter),
89 s0(x.s0), xh0(x.xh0), d_xh0_ds(x.d_xh0_ds),
90 has_init(x.has_init), has_refresh(x.has_refresh), has_spherical(x.has_spherical),
91 u_ownership(x.u_ownership)
93 std::string metric = (!has_spherical) ?
"orthogonal" :
"spherical";
96 template<
class Problem>
100 check_macro (metric ==
"orthogonal" || metric ==
"spherical",
"invalid metric="<<metric);
101 has_spherical = (metric !=
"orthogonal");
103 template<
class Problem>
108 xh0.second = pb.initial (restart);
109 xh0.first = pb.parameter();
110 typename Problem::value_type d_uh0_ds = pb.direction(xh0.second);
111 float_type c = 1/sqrt(1 + sqr(pb.space_norm(d_uh0_ds)));
113 d_xh0_ds.second =
c*d_uh0_ds;
117 template<
class Problem>
121 pb.set_parameter (xh.first);
123 if (!has_spherical) {
124 mrh.first = space_dot(d_xh0_ds, xh-xh0) - (s-s0);
126 mrh.first = sqr(space_norm(xh-xh0)) - sqr(s-s0);
128 mrh.second = pb.residue (xh.second);
131 template<
class Problem>
132 typename solver::determinant_type
134 pb.set_parameter (xh.first);
135 m_df_d_parameter = pb.derivative_versus_parameter(xh.second);
136 form df_du = pb.derivative(xh.second);
137 u_ownership = df_du.uu().row_ownership();
141 if (!has_spherical) {
145 typename Problem::value_type duh = xh.second - xh0.second;
147 c = 2*(xh.first - xh0.first);
149 A1 = {{ df_du.uu(), b2 },
156 template<
class Problem>
167 delta_xh.first = (size_t(X.comm().rank()) == iproc0) ? X[X.size()-1] : 0;
168 #ifdef _RHEOLEF_HAVE_MPI
169 mpi::broadcast (X.comm(), delta_xh.first, iproc0);
170 #endif // _RHEOLEF_HAVE_MPI
173 template<
class Problem>
184 mgh.first = (size_t(MGH.comm().rank()) == iproc0) ? MGH[MGH.size()-1] : 0;
185 #ifdef _RHEOLEF_HAVE_MPI
186 mpi::broadcast (MGH.comm(), mgh.first, iproc0);
187 #endif // _RHEOLEF_HAVE_MPI
190 template<
class Problem>
197 if (p_err) *p_err << std::endl << std::endl;
198 return ((
status == 0 && res <= opts.
tol) || sqr(res) <= opts.
tol) ? 0 : 1;
200 template<
class Problem>
204 if (!has_refresh)
return d_xh0_ds;
205 update_derivative(xh);
206 return -derivative_solve(derivative_versus_parameter(xh));
208 template<
class Problem>
216 template<
class Problem>
221 if (!has_spherical) {
224 m_df_ds.first = -2*(s-s0);
229 template<
class Problem>
233 return xh.first*yh.first + pb.space_dot(xh.second, yh.second);
235 template<
class Problem>
239 return sqrt (space_dot(xh,xh));
241 template<
class Problem>
245 return sqrt (sqr(mrh.first) + sqr(pb.dual_space_norm(mrh.second)));
247 template<
class Problem>
250 static const bool trace =
true;
252 pb.set_parameter(xh.first);
255 if (
n==0)
derr <<
"#[keller] n s "<<pb.parameter_name()<<
" |x|"<<std::endl;
256 derr <<
"[keller] " <<
n <<
" " << s <<
" " << xh.first <<
" "
257 << pb.space_norm(xh.second) << std::endl;
260 os << std::endl <<
"#s " << s << std::endl;
261 pb.put(os,xh.second);
262 pb.set_parameter(prev_param);
265 template<
class Problem>
270 pb.get (is, xh.second);
271 xh.first = pb.parameter();
272 if (!has_init) { initial(); }
278 template<
class Problem>
280 :
public keller<Problem,std::false_type>
286 keller (
const Problem& pb, std::string metric=
"orthogonal") :
base(pb,metric) {}
292 template<
class Problem>
295 base::pb.adapt (xh.second, aopt);
296 base::xh0.second = base::pb.reinterpolate (base::xh0.second);
297 base::d_xh0_ds.second = base::pb.reinterpolate (base::d_xh0_ds.second);
299 template<
class Problem>
302 base::pb.reset_geo (xh.second);
303 base::xh0.second = base::pb.reinterpolate (base::xh0.second);
304 base::d_xh0_ds.second = base::pb.reinterpolate (base::d_xh0_ds.second);
306 template<
class Problem>
311 yh.second = base::pb.reinterpolate (yh.second);
316 #endif // _RHEOLEF_KELLER_H