Rheolef  7.2
an efficient C++ finite element environment
harten0.h
Go to the documentation of this file.
1 // w=harten0(t,x) such that: w - g(w) = 0 where g(w) = sin(pi*(x-w*t))
26 struct 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; }
90 protected:
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)")
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
Float epsilon