Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
harten0.h
Go to the documentation of this file.
1
25// w=harten0(t,x) such that: w - g(w) = 0 where g(w) = sin(pi*(x-w*t))
26struct harten0 {
27 harten0 (Float t1) : t0(t1), pi(acos(Float(-1))), verbose(false) {}
28 Float operator() (const point& x) const {
29 Float x0 = fabs(x[0]);
30 // find a good Newton's method starting point w0 vs w_star / g'(w_star) = 0
31 Float wmin = 0, wmax = 1, w0 = 0.5;
32 if (t0 > 0) {
33 Float x_star = (x0 < 0.5) ? x0+0.5 : x0-0.5;
34 Float w_star = x_star/t0;
35 Float w_star2 = (x_star+1)/t0;
36 Float g_star = sin(pi*(x0-w_star*t0));
37 if (w_star > wmax) {
38 w0 = 0.5;
39 } else {
40 if (g_star < w_star) {
41 w0 = (w_star + wmin)/2;
42 } else {
43 if (w_star2 > wmax) {
44 w0 = (w_star + wmax)/2;
45 } else {
46 w0 = (w_star + w_star2)/2;
47 }
48 }
49 }
50 }
51 // use a list of starting points for Newton
52 Float wini_list[] = {w0, wmax, wmin};
53 // loop
54 Float r = 0;
55 size_t i = 0;
56 int status = 0;
57 const Float tol = 1e2*numeric_limits<Float>::epsilon();
58 Float w = 0;
59 for (Float wini : wini_list) {
60 if (verbose) derr << "[harten] wmin="<<wmin<<", wmax="<<wmax<<", wini="<<wini<<endl;
61 if (verbose) derr << "[harten] i r dw w, for t0=" << t0 << ", x0="<<x0<< endl;
62 w = wini;
63 i = 0;
64 status = 0;
65 do {
66 r = f(w,t0,x0);
67 Float dw = -r/df_dw(w,t0,x0);
68 if (verbose) derr << "[harten] " << i << " " << fabs(r)
69 << " " << dw << " " << w << endl;
70 if (fabs(r) <= tol && fabs(dw) <= tol) { break; }
71 if (++i >= max_iter) { status = 1; break; }
72 if (w+dw > 0 && w+dw < wmax) {
73 w += dw;
74 } else if (w+dw <= 0) {
75 w = (w + wmin)/2;
76 check_macro (1+w != w, "Newton iteration: machine precision problem.");
77 } else {
78 w = (w + wmax)/2;
79 }
80 } while (true);
81 if (status == 0) break;
82 if (verbose) derr << "[harten] failed: restart with another starting point..." << endl;
83 }
84 check_macro (status == 0, "t = " << t0 << ", x = " << x0 << " : Newton iteration " << i
85 << ": precision " << tol << " not reached: " << fabs(r));
86 return (x[0] >= 0) ? w : -w;
87 }
88 void set_verbose() { verbose = true; }
89 void set_noverbose() { verbose = false; }
90protected:
91 Float f (Float w, Float t, Float x) const { return w - sin(pi*(x-w*t)); }
92 Float df_dw (Float w, Float t, Float x) const { return 1 + pi*t*cos(pi*(x-w*t)); }
93 static const size_t max_iter = 100;
95 bool verbose;
96};
see the Float page for the full documentation
see the point 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)")
Definition cavity_dg.h:29
Float t0
Definition harten0.h:94
static const size_t max_iter
Definition harten0.h:93
void set_noverbose()
Definition harten0.h:89
harten0(Float t1)
Definition harten0.h:27
bool verbose
Definition harten0.h:95
Float operator()(const point &x) const
Definition harten0.h:28
Float f(Float w, Float t, Float x) const
Definition harten0.h:91
void set_verbose()
Definition harten0.h:88
Float df_dw(Float w, Float t, Float x) const
Definition harten0.h:92
Float pi
Definition harten0.h:94