-
Notifications
You must be signed in to change notification settings - Fork 8
/
interpolate.H
128 lines (90 loc) · 4.04 KB
/
interpolate.H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#ifndef INTERPOLATE_H
#define INTERPOLATE_H
#include <read_model.H>
AMREX_INLINE AMREX_GPU_HOST_DEVICE
int
locate(const Real r, const initial_model_t& initial_model) {
int loc;
if (r <= initial_model.r(0)) {
loc = 0;
} else if (r > initial_model.r(initial_model.npts-2)) {
loc = initial_model.npts-1;
} else {
int ilo = 0;
int ihi = initial_model.npts-2;
while (ilo+1 != ihi) {
int imid = (ilo + ihi) / 2;
if (r <= initial_model.r(imid)) {
ihi = imid;
} else {
ilo = imid;
}
}
loc = ihi;
}
return loc;
}
AMREX_INLINE AMREX_GPU_HOST_DEVICE
Real
interpolate(const Real r, const int var_index, const initial_model_t& initial_model, const bool extrapolate_top = false) {
// find the value of model_state component var_index at point r
// using linear interpolation. Eventually, we can do something
// fancier here.
int id = locate(r, initial_model);
Real slope;
Real interp;
if (id == 0) {
// base of the model
slope = (initial_model.state(id+1, var_index) - initial_model.state(id, var_index)) /
(initial_model.r(id+1) - initial_model.r(id));
interp = slope * (r - initial_model.r(id)) + initial_model.state(id, var_index);
// safety check to make sure interp lies within the bounding points
Real minvar = amrex::min(initial_model.state(id+1, var_index),
initial_model.state(id, var_index));
Real maxvar = amrex::max(initial_model.state(id+1, var_index),
initial_model.state(id, var_index));
interp = amrex::max(interp, minvar);
interp = amrex::min(interp, maxvar);
} else if (id == initial_model.npts-1) {
// top of the model
slope = (initial_model.state(id, var_index) - initial_model.state(id-1, var_index)) /
(initial_model.r(id) - initial_model.r(id-1));
interp = slope * (r - initial_model.r(id)) + initial_model.state(id, var_index);
// safety check to make sure interp lies within the bounding points -- only do this
// if !extrapolate_top
if (! extrapolate_top) {
Real minvar = amrex::min(initial_model.state(id-1, var_index),
initial_model.state(id, var_index));
Real maxvar = amrex::max(initial_model.state(id-1, var_index),
initial_model.state(id, var_index));
interp = amrex::max(interp, minvar);
interp = amrex::min(interp, maxvar);
}
} else {
if (r >= initial_model.r(id)) {
slope = (initial_model.state(id+1, var_index) -
initial_model.state(id, var_index)) /
(initial_model.r(id+1) - initial_model.r(id));
interp = slope * (r - initial_model.r(id)) + initial_model.state(id, var_index);
Real minvar = amrex::min(initial_model.state(id+1, var_index),
initial_model.state(id, var_index));
Real maxvar = amrex::max(initial_model.state(id+1, var_index),
initial_model.state(id, var_index));
interp = amrex::max(interp, minvar);
interp = amrex::min(interp, maxvar);
} else {
slope = (initial_model.state(id, var_index) -
initial_model.state(id-1, var_index)) /
(initial_model.r(id) - initial_model.r(id-1));
interp = slope * (r - initial_model.r(id)) + initial_model.state(id, var_index);
Real minvar = amrex::min(initial_model.state(id-1, var_index),
initial_model.state(id, var_index));
Real maxvar = amrex::max(initial_model.state(id-1, var_index),
initial_model.state(id, var_index));
interp = amrex::max(interp, minvar);
interp = amrex::min(interp, maxvar);
}
}
return interp;
}
#endif