Rheolef
7.1
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;
94
Float
t0
,
pi
;
95
bool
verbose
;
96
};
harten0::df_dw
Float df_dw(Float w, Float t, Float x) const
Definition:
harten0.h:92
harten0::verbose
bool verbose
Definition:
harten0.h:95
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
harten0::max_iter
static const size_t max_iter
Definition:
harten0.h:93
harten0::t0
Float t0
Definition:
harten0.h:94
harten0::f
Float f(Float w, Float t, Float x) const
Definition:
harten0.h:91
harten0::set_verbose
void set_verbose()
Definition:
harten0.h:88
Float
see the Float page for the full documentation
point
see the point page for the full documentation
harten0::set_noverbose
void set_noverbose()
Definition:
harten0.h:89
harten0::operator()
Float operator()(const point &x) const
Definition:
harten0.h:28
harten0::harten0
harten0(Float t1)
Definition:
harten0.h:27
epsilon
Float epsilon
Definition:
transmission_error.cc:25
harten0::pi
Float pi
Definition:
harten0.h:94
mkgeo_contraction.status
status
Definition:
mkgeo_contraction.sh:290
harten0
Definition:
harten0.h:26