Rheolef  7.1
an efficient C++ finite element environment
keller.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_KELLER_H
2 #define _RHEOLEF_KELLER_H
3 
24 #include "rheolef/newton_add_missing.h"
25 #include "rheolef/keller_details.h"
26 
27 namespace rheolef {
28 
29 // switch depending on existing Problem::adapt() method
31 template<class Problem, class Sfinae = typename details::has_inherited_member_adapt<Problem>::type>
32 class keller {};
33 
34 // -----------------------------
35 // class without mesh adaptation
36 // -----------------------------
37 template<class Problem>
38 class keller<Problem,std::false_type> {
39 public:
40  typedef Float float_type;
42  keller(const Problem& pb, std::string metric="orthogonal");
44  void reset(const Problem& pb, std::string metric="orthogonal");
45  float_type parameter() const { return s; }
46  std::string parameter_name() const { return "s"; }
47  Problem& get_problem() { return pb; }
48  void set_parameter(float_type s1) { s = s1; }
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;
54  int solve (value_type& uh, const continuation_option& opts) 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;
61  odiststream& put (odiststream& os, const value_type& xh) const;
62  idiststream& get (idiststream& is, value_type& xh);
63  bool stop (const value_type& xh) const { return pb.stop(xh.second); }
64 // adapt (dummy):
65  void adapt (const value_type&, const adapt_option&) {}
66  void reset_geo (const value_type&) {}
67  value_type reinterpolate (const value_type& xh) { return xh; }
68 protected:
72  mutable solver sA1;
73  mutable typename Problem::value_type m_df_d_parameter;
74  mutable float_type s0;
75  mutable value_type xh0, d_xh0_ds;
76  mutable bool has_init, has_refresh, has_spherical;
78 };
79 template<class Problem>
80 keller<Problem,std::false_type>::keller (const Problem& pb1, std::string metric)
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()
83 {
84  reset (pb1, metric);
85 }
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)
92 {
93  std::string metric = (!has_spherical) ? "orthogonal" : "spherical";
94  reset (pb, metric);
95 }
96 template<class Problem>
97 void
98 keller<Problem,std::false_type>::reset (const Problem& pb1, std::string metric) {
99  pb = pb1;
100  check_macro (metric == "orthogonal" || metric == "spherical", "invalid metric="<<metric);
101  has_spherical = (metric != "orthogonal");
102 }
103 template<class Problem>
104 typename
106 keller<Problem,std::false_type>::initial (std::string restart) const {
107  s0 = 0;
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)));
112  d_xh0_ds.first = c;
113  d_xh0_ds.second = c*d_uh0_ds;
114  has_init = true;
115  return xh0;
116 }
117 template<class Problem>
118 typename
121  pb.set_parameter (xh.first);
122  value_type mrh;
123  if (!has_spherical) {
124  mrh.first = space_dot(d_xh0_ds, xh-xh0) - (s-s0);
125  } else {
126  mrh.first = sqr(space_norm(xh-xh0)) - sqr(s-s0);
127  }
128  mrh.second = pb.residue (xh.second);
129  return mrh;
130 }
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();
138  vec<float_type> b2 = details::get_unknown (m_df_d_parameter, u_ownership);
139  vec<float_type> b1;
140  Float c;
141  if (!has_spherical) {
142  b1 = details::get_unknown (pb.massify(d_xh0_ds.second), u_ownership);
143  c = d_xh0_ds.first;
144  } else {
145  typename Problem::value_type duh = xh.second - xh0.second;
146  b1 = details::get_unknown (pb.massify(Float(2)*duh), u_ownership);
147  c = 2*(xh.first - xh0.first);
148  }
149  A1 = {{ df_du.uu(), b2 },
150  { trans(b1), c }};
151  solver_option sopt;
152  sopt.compute_determinant = true;
153  sA1 = solver_basic<float_type> (A1, sopt);
154  return sA1.det();
155 }
156 template<class Problem>
157 typename
160  vec<float_type> MB = { details::get_unknown(mrh.second,u_ownership),
161  mrh.first };
162  vec<float_type> X = sA1.solve(MB);
163  value_type delta_xh;
164  delta_xh.second = mrh.second; details::reset (delta_xh.second);
165  details::set_unknown (delta_xh.second, X);
166  size_t iproc0 = 0;
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
171  return delta_xh;
172 }
173 template<class Problem>
174 typename
177  vec<float_type> MRH = { details::get_unknown(mrh.second,u_ownership),
178  mrh.first };
179  vec<float_type> MGH = A1.trans_mult(MRH);
180  value_type mgh = mrh; // for std::valarray alloc
181  mgh.second = mrh.second; details::reset (mgh.second);
182  details::set_unknown (mgh.second, MGH);
183  size_t iproc0 = 0;
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
188  return mgh;
189 }
190 template<class Problem>
191 int
193  float_type res = opts.tol;
194  size_t n_iter = opts.max_iter;
195  odiststream* p_err = &derr;
196  int status = damped_newton (*this, uh, res, n_iter, p_err);
197  if (p_err) *p_err << std::endl << std::endl;
198  return ((status == 0 && res <= opts.tol) || sqr(res) <= opts.tol) ? 0 : 1;
199 }
200 template<class Problem>
201 typename
204  if (!has_refresh) return d_xh0_ds;
205  update_derivative(xh);
206  return -derivative_solve(derivative_versus_parameter(xh));
207 }
208 template<class Problem>
209 void
211  has_refresh = true;
212  d_xh0_ds = d_xh_ds;
213  xh0 = xh;
214  s0 = s;
215 }
216 template<class Problem>
217 typename
220  value_type m_df_ds;
221  if (!has_spherical) {
222  m_df_ds.first = -1;
223  } else {
224  m_df_ds.first = -2*(s-s0);
225  }
226  m_df_ds.second = xh.second; details::reset(m_df_ds.second); // allocate+reset
227  return m_df_ds;
228 }
229 template<class Problem>
230 typename
233  return xh.first*yh.first + pb.space_dot(xh.second, yh.second);
234 }
235 template<class Problem>
236 typename
239  return sqrt (space_dot(xh,xh));
240 }
241 template<class Problem>
242 typename
245  return sqrt (sqr(mrh.first) + sqr(pb.dual_space_norm(mrh.second)));
246 }
247 template<class Problem>
250  static const bool trace = true;
251  float_type prev_param = pb.parameter();
252  pb.set_parameter(xh.first);
253  if (trace) {
254  static size_t n = 0;
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;
258  n++;
259  }
260  os << std::endl << "#s " << s << std::endl;
261  pb.put(os,xh.second);
262  pb.set_parameter(prev_param);
263  return os;
264 }
265 template<class Problem>
266 idiststream&
268  is >> catchmark("s") >> s;
269  if (!is) return is;
270  pb.get (is, xh.second);
271  xh.first = pb.parameter();
272  if (!has_init) { initial(); }
273  return is;
274 }
275 // -----------------------------
276 // class with mesh adaptation
277 // -----------------------------
278 template<class Problem>
279 class keller<Problem,std::true_type>
280  : public keller<Problem,std::false_type>
281 {
282 public:
284  typedef typename base::float_type float_type;
285  typedef typename base::value_type value_type;
286  keller (const Problem& pb, std::string metric="orthogonal") : base(pb,metric) {}
287 // adapt:
288  void adapt (const value_type& xh, const adapt_option& aopt);
289  void reset_geo (const value_type& xh);
290  value_type reinterpolate (const value_type& xh);
291 };
292 template<class Problem>
293 void
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);
298 }
299 template<class Problem>
300 void
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);
305 }
306 template<class Problem>
307 typename
310  typename keller<Problem>::value_type yh = xh;
311  yh.second = base::pb.reinterpolate (yh.second);
312  return yh;
313 }
314 
315 } // namespace rheolef
316 #endif // _RHEOLEF_KELLER_H
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
put
void put(idiststream &in, odiststream &out, bool do_proj, bool do_lumped_mass, bool def_fill_opt, size_type extract_id, const Float &scale_value, const std::pair< Float, Float > &u_range, render_type render)
Definition: branch.cc:500
rheolef::keller< Problem, std::false_type >
Definition: keller.h:38
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:172
form
see the form page for the full documentation
rheolef::keller< Problem, std::true_type >::value_type
base::value_type value_type
Definition: keller.h:285
rheolef::keller< Problem, std::false_type >::pb
details::add_missing_continuation< Problem > pb
Definition: keller.h:70
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
solver
see the solver page for the full documentation
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition: tiny_lu.h:155
rheolef::details::add_missing_continuation
Definition: newton_add_missing.h:275
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::catchmark
see the catchmark page for the full documentation
Definition: catchmark.h:67
rheolef::details::reset
void reset(T &x)
Definition: keller_details.h:31
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::keller< Problem, std::false_type >::u_ownership
distributor u_ownership
Definition: keller.h:77
rheolef::keller< Problem, std::false_type >::parameter
float_type parameter() const
Definition: keller.h:45
rheolef::solver_basic
Definition: solver.h:264
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::float_type
typename float_traits< value_type >::type float_type
Definition: field_expr_recursive.h:501
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 >::get_problem
Problem & get_problem()
Definition: keller.h:47
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
rheolef::keller< Problem, std::false_type >::set_parameter
void set_parameter(float_type s1)
Definition: keller.h:48
rheolef::keller< Problem, std::false_type >::stop
bool stop(const value_type &xh) const
Definition: keller.h:63
rheolef::continuation_option::max_iter
size_t max_iter
Definition: continuation_option.h:81
rheolef::value_type
result_type value_type
Definition: field_expr_recursive.h:499
rheolef::csr
see the csr page for the full documentation
Definition: csr.h:317
rheolef::keller< Problem, std::false_type >::sA1
solver sA1
Definition: keller.h:72
rheolef::keller< Problem, std::true_type >::float_type
base::float_type float_type
Definition: keller.h:284
rheolef::keller< Problem, std::false_type >::reinterpolate
value_type reinterpolate(const value_type &xh)
Definition: keller.h:67
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::details::get_unknown
size_t get_unknown(const T &x)
Definition: keller_details.h:65
rheolef::keller< Problem, std::false_type >::adapt
void adapt(const value_type &, const adapt_option &)
Definition: keller.h:65
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::derr
odiststream derr(cerr)
see the diststream page for the full documentation
Definition: diststream.h:436
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
Float
see the Float page for the full documentation
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::details::pair_with_linear_algebra< float_type, typename Problem::value_type >
rheolef::keller< Problem, std::false_type >::reset_geo
void reset_geo(const value_type &)
Definition: keller.h:66
rheolef::keller< Problem, std::true_type >::base
keller< Problem, std::false_type > base
Definition: keller.h:283
rheolef::keller< Problem, std::false_type >::s
float_type s
Definition: keller.h:69
rheolef::solver_option::compute_determinant
bool compute_determinant
Definition: solver_option.h:169
rheolef::keller< Problem, std::false_type >::s0
float_type s0
Definition: keller.h:74
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 >::float_type
Float float_type
Definition: keller.h:40
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
rheolef::keller
see the continuation page for the full documentation
Definition: keller.h:32
rheolef::keller< Problem, std::false_type >::A1
csr< float_type > A1
Definition: keller.h:71
rheolef::keller< Problem, std::false_type >::has_spherical
bool has_spherical
Definition: keller.h:76
rheolef::std
Definition: vec_expr_v2.h:402
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
rheolef::details::set_unknown
void set_unknown(T &x, const T &value)
Definition: keller_details.h:119
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
rheolef::keller< Problem, std::false_type >::xh0
value_type xh0
Definition: keller.h:75
rheolef::keller< Problem, std::true_type >::keller
keller(const Problem &pb, std::string metric="orthogonal")
Definition: keller.h:286
rheolef::keller< Problem, std::false_type >::parameter_name
std::string parameter_name() const
Definition: keller.h:46
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155