Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
keller.h
Go to the documentation of this file.
1#ifndef _RHEOLEF_KELLER_H
2#define _RHEOLEF_KELLER_H
23
24#include "rheolef/newton_add_missing.h"
25#include "rheolef/keller_details.h"
26
27namespace rheolef {
28
29// switch depending on existing Problem::adapt() method
31template<class Problem, class Sfinae = typename details::has_inherited_member_adapt<Problem>::type>
32class keller {};
33
34// -----------------------------
35// class without mesh adaptation
36// -----------------------------
37template<class Problem>
38class keller<Problem,std::false_type> {
39public:
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; }
68protected:
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};
79template<class Problem>
80keller<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}
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)
92{
93 std::string metric = (!has_spherical) ? "orthogonal" : "spherical";
94 reset (pb, metric);
95}
96template<class Problem>
97void
98keller<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}
103template<class Problem>
104typename
106keller<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}
117template<class Problem>
118typename
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}
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();
138 vec<float_type> b2 = details::get_unknown (m_df_d_parameter, u_ownership);
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}
156template<class Problem>
157typename
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}
173template<class Problem>
174typename
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}
190template<class Problem>
191int
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}
200template<class Problem>
201typename
204 if (!has_refresh) return d_xh0_ds;
205 update_derivative(xh);
206 return -derivative_solve(derivative_versus_parameter(xh));
207}
208template<class Problem>
209void
211 has_refresh = true;
212 d_xh0_ds = d_xh_ds;
213 xh0 = xh;
214 s0 = s;
215}
216template<class Problem>
217typename
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}
229template<class Problem>
230typename
233 return xh.first*yh.first + pb.space_dot(xh.second, yh.second);
234}
235template<class Problem>
236typename
239 return sqrt (space_dot(xh,xh));
240}
241template<class Problem>
242typename
245 return sqrt (sqr(mrh.first) + sqr(pb.dual_space_norm(mrh.second)));
246}
247template<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}
265template<class Problem>
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// -----------------------------
278template<class Problem>
279class keller<Problem,std::true_type>
280 : public keller<Problem,std::false_type>
281{
282public:
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};
292template<class Problem>
293void
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}
299template<class Problem>
300void
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}
306template<class Problem>
307typename
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
see the Float page for the full documentation
see the form page for the full documentation
see the catchmark page for the full documentation
Definition catchmark.h:67
see the csr page for the full documentation
Definition csr.h:317
see the distributor page for the full documentation
Definition distributor.h:69
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
bool stop(const value_type &xh) const
Definition keller.h:63
details::add_missing_continuation< Problem > pb
Definition keller.h:70
void reset_geo(const value_type &)
Definition keller.h:66
details::pair_with_linear_algebra< float_type, typename Problem::value_type > value_type
Definition keller.h:41
value_type reinterpolate(const value_type &xh)
Definition keller.h:67
void adapt(const value_type &, const adapt_option &)
Definition keller.h:65
keller(const Problem &pb, std::string metric="orthogonal")
Definition keller.h:286
keller< Problem, std::false_type > base
Definition keller.h:283
see the continuation page for the full documentation
Definition keller.h:32
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
see the solver_option page for the full documentation
see the vec page for the full documentation
Definition vec.h:79
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
Definition adapt.cc:179
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)
Definition tiny_lu.h:92
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
STL namespace.
field residue(Float p, const field &uh)
adapt_option: see the adapt page for the full documentation
Definition adapt.h:147
see the continuation_option page for the full documentation