1#ifndef _RHEOLEF_KELLER_H
2#define _RHEOLEF_KELLER_H
24#include "rheolef/newton_add_missing.h"
25#include "rheolef/keller_details.h"
31template<class Problem, class Sfinae = typename details::has_inherited_member_adapt<Problem>::type>
37template<
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;
50 void refresh (float_type s,
const value_type& xh,
const value_type& d_xh_ds)
const;
51 value_type direction (value_type& xh)
const;
52 value_type derivative_versus_parameter (
const value_type& xh)
const;
53 value_type
residue (
const value_type& xh)
const;
55 solver::determinant_type update_derivative (
const value_type& xh)
const;
56 value_type derivative_solve (
const value_type& mrh)
const;
57 value_type derivative_trans_mult (
const value_type& mrh)
const;
58 float_type space_norm (
const value_type& xh)
const;
59 float_type dual_space_norm (
const value_type& mrh)
const;
60 float_type space_dot (
const value_type& xh,
const value_type& yh)
const;
76 mutable bool has_init, has_refresh, has_spherical;
79template<
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()
86template<
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";
96template<
class Problem>
100 check_macro (metric ==
"orthogonal" || metric ==
"spherical",
"invalid metric="<<metric);
101 has_spherical = (metric !=
"orthogonal");
103template<
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;
117template<
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);
131template<
class Problem>
132typename 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 },
156template<
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);
173template<
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);
190template<
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;
200template<
class Problem>
204 if (!has_refresh)
return d_xh0_ds;
205 update_derivative(xh);
206 return -derivative_solve(derivative_versus_parameter(xh));
208template<
class Problem>
216template<
class Problem>
221 if (!has_spherical) {
224 m_df_ds.first = -2*(s-s0);
229template<
class Problem>
233 return xh.first*yh.first + pb.space_dot(xh.second, yh.second);
235template<
class Problem>
239 return sqrt (space_dot(xh,xh));
241template<
class Problem>
245 return sqrt (sqr(mrh.first) + sqr(pb.dual_space_norm(mrh.second)));
247template<
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);
265template<
class Problem>
270 pb.get (is, xh.second);
271 xh.first = pb.parameter();
272 if (!has_init) { initial(); }
278template<
class Problem>
280 :
public keller<Problem,std::false_type>
286 keller (
const Problem& pb, std::string metric=
"orthogonal") :
base(pb,metric) {}
289 void reset_geo (
const value_type& xh);
290 value_type reinterpolate (
const value_type& xh);
292template<
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);
299template<
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);
306template<
class Problem>
311 yh.second = base::pb.reinterpolate (yh.second);
see the Float page for the full documentation
see the catchmark page for the full documentation
see the csr page for the full documentation
see the distributor page for the full documentation
idiststream: see the diststream page for the full documentation
bool stop(const value_type &xh) const
std::string parameter_name() const
details::add_missing_continuation< Problem > pb
Problem::value_type m_df_d_parameter
float_type parameter() const
void reset_geo(const value_type &)
details::pair_with_linear_algebra< float_type, typename Problem::value_type > value_type
void set_parameter(float_type s1)
value_type reinterpolate(const value_type &xh)
void adapt(const value_type &, const adapt_option &)
base::value_type value_type
keller(const Problem &pb, std::string metric="orthogonal")
base::float_type float_type
keller< Problem, std::false_type > base
see the continuation page for the full documentation
odiststream: see the diststream page for the full documentation
see the solver_option page for the full documentation
see the vec page for the full documentation
see the solver page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
size_t get_unknown(const T &x)
void set_unknown(T &x, const T &value)
This file is part of Rheolef.
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
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
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
field residue(Float p, const field &uh)
adapt_option: see the adapt page for the full documentation
see the continuation_option page for the full documentation