tColumnInterpolation.cc - pism - [fork] customized build of PISM, the parallel … | |
git clone git://src.adamsgaard.dk/pism | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
tColumnInterpolation.cc (7168B) | |
--- | |
1 /* Copyright (C) 2014, 2015 PISM Authors | |
2 * | |
3 * This file is part of PISM. | |
4 * | |
5 * PISM is free software; you can redistribute it and/or modify it under… | |
6 * terms of the GNU General Public License as published by the Free Soft… | |
7 * Foundation; either version 3 of the License, or (at your option) any … | |
8 * version. | |
9 * | |
10 * PISM is distributed in the hope that it will be useful, but WITHOUT A… | |
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FIT… | |
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more | |
13 * details. | |
14 * | |
15 * You should have received a copy of the GNU General Public License | |
16 * along with PISM; if not, write to the Free Software | |
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301… | |
18 */ | |
19 | |
20 #include "ColumnInterpolation.hh" | |
21 | |
22 #include <cmath> | |
23 | |
24 namespace pism { | |
25 | |
26 ColumnInterpolation::ColumnInterpolation(const std::vector<double> &new_… | |
27 const std::vector<double> &new_… | |
28 : m_z_fine(new_z_fine), m_z_coarse(new_z_coarse) { | |
29 init_interpolation(); | |
30 } | |
31 | |
32 std::vector<double> ColumnInterpolation::coarse_to_fine(const std::vecto… | |
33 unsigned int ks)… | |
34 std::vector<double> result(Mz_fine()); | |
35 coarse_to_fine(&input[0], ks, &result[0]); | |
36 return result; | |
37 } | |
38 | |
39 void ColumnInterpolation::coarse_to_fine(const double *input, unsigned i… | |
40 if (m_use_linear_interpolation) { | |
41 coarse_to_fine_linear(input, ks, result); | |
42 } else { | |
43 coarse_to_fine_quadratic(input, ks, result); | |
44 } | |
45 } | |
46 | |
47 void ColumnInterpolation::coarse_to_fine_linear(const double *input, uns… | |
48 double *result) const { | |
49 const unsigned int Mzfine = Mz_fine(); | |
50 const unsigned int Mzcoarse = Mz_coarse(); | |
51 | |
52 for (unsigned int k = 0; k < Mzfine; ++k) { | |
53 if (k > ks) { | |
54 result[k] = input[m_coarse2fine[k]]; | |
55 continue; | |
56 } | |
57 | |
58 unsigned int m = m_coarse2fine[k]; | |
59 | |
60 // extrapolate (if necessary): | |
61 if (m == Mzcoarse - 1) { | |
62 result[k] = input[Mzcoarse - 1]; | |
63 continue; | |
64 } | |
65 | |
66 const double incr = (m_z_fine[k] - m_z_coarse[m]) / (m_z_coarse[m + … | |
67 result[k] = input[m] + incr * (input[m + 1] - input[m]); | |
68 } | |
69 } | |
70 | |
71 void ColumnInterpolation::coarse_to_fine_quadratic(const double *input, … | |
72 double *result) const… | |
73 unsigned int k = 0, m = 0; | |
74 const unsigned int Mz = Mz_coarse(); | |
75 for (m = 0; m < Mz - 2 and k <= ks; ++m) { | |
76 | |
77 const double | |
78 z0 = m_z_coarse[m], | |
79 z1 = m_z_coarse[m + 1], | |
80 dz_inv = m_constants[3 * m + 0], // = 1.0 / (z1 - z0) | |
81 dz1_inv = m_constants[3 * m + 1], // = 1.0 / (z2 - z0) | |
82 dz2_inv = m_constants[3 * m + 2], // = 1.0 / (z2 - z1) | |
83 f0 = input[m], | |
84 f1 = input[m + 1], | |
85 f2 = input[m + 2]; | |
86 | |
87 const double | |
88 d1 = (f1 - f0) * dz_inv, | |
89 d2 = (f2 - f0) * dz1_inv, | |
90 b = (d2 - d1) * dz2_inv, | |
91 a = d1 - b * (z1 - z0), | |
92 c = f0; | |
93 | |
94 for (; m_z_fine[k] < z1 and k <= ks; ++k) { | |
95 const double s = m_z_fine[k] - z0; | |
96 | |
97 result[k] = s * (a + b * s) + c; | |
98 } | |
99 } // m-loop | |
100 | |
101 // check if we got to the end of the m-loop and use linear | |
102 // interpolation between the remaining 2 coarse levels | |
103 if (m == Mz - 2) { | |
104 const double | |
105 z0 = m_z_coarse[m], | |
106 z1 = m_z_coarse[m + 1], | |
107 f0 = input[m], | |
108 f1 = input[m + 1], | |
109 lambda = (f1 - f0) / (z1 - z0); | |
110 | |
111 for (; m_z_fine[k] < z1; ++k) { | |
112 result[k] = f0 + lambda * (m_z_fine[k] - z0); | |
113 } | |
114 } | |
115 | |
116 // fill the rest using constant extrapolation | |
117 const double f0 = input[Mz - 1]; | |
118 for (; k <= ks; ++k) { | |
119 result[k] = f0; | |
120 } | |
121 } | |
122 | |
123 std::vector<double> ColumnInterpolation::fine_to_coarse(const std::vecto… | |
124 std::vector<double> result(Mz_coarse()); | |
125 fine_to_coarse(&input[0], &result[0]); | |
126 return result; | |
127 } | |
128 | |
129 void ColumnInterpolation::fine_to_coarse(const double *input, double *re… | |
130 const unsigned int N = Mz_coarse(); | |
131 | |
132 for (unsigned int k = 0; k < N - 1; ++k) { | |
133 const int m = m_fine2coarse[k]; | |
134 | |
135 const double increment = (m_z_coarse[k] - m_z_fine[m]) / (m_z_fine[m… | |
136 result[k] = input[m] + increment * (input[m + 1] - input[m]); | |
137 } | |
138 | |
139 result[N - 1] = input[m_fine2coarse[N - 1]]; | |
140 } | |
141 | |
142 unsigned int ColumnInterpolation::Mz_coarse() const { | |
143 return m_z_coarse.size(); | |
144 } | |
145 | |
146 unsigned int ColumnInterpolation::Mz_fine() const { | |
147 return m_z_fine.size(); | |
148 } | |
149 | |
150 double ColumnInterpolation::dz_fine() const { | |
151 return m_z_fine[1] - m_z_fine[0]; | |
152 } | |
153 | |
154 const std::vector<double>& ColumnInterpolation::z_fine() const { | |
155 return m_z_fine; | |
156 } | |
157 | |
158 const std::vector<double>& ColumnInterpolation::z_coarse() const { | |
159 return m_z_coarse; | |
160 } | |
161 | |
162 /*! | |
163 * Given two 1D grids, `z_input` and `z_output`, for each `z_output` | |
164 * index `k` we find an index `m` so that | |
165 * | |
166 * `z_input[m] < z_output[k] <= z_input[m+1]` | |
167 * | |
168 * In other words, we look for two consecutive points in the input | |
169 * grid that bracket a point in the output grid. | |
170 * | |
171 * This function sets `result[k] = m`. This information is then used | |
172 * to interpolate from the grid defined by `z_input` to the one | |
173 * defined by `z_output`. | |
174 * | |
175 * We use constant extrapolation outside the range defined by `z_input`. | |
176 */ | |
177 static std::vector<unsigned int> init_interpolation_indexes(const std::v… | |
178 const std::v… | |
179 std::vector<unsigned int> result(z_output.size()); | |
180 | |
181 unsigned int m = 0; | |
182 for (unsigned int k = 0; k < z_output.size(); ++k) { | |
183 | |
184 if (z_output[k] <= z_input.front()) { | |
185 result[k] = 0; | |
186 continue; | |
187 } else if (z_output[k] >= z_input.back()) { | |
188 result[k] = z_input.size() - 1; | |
189 continue; | |
190 } | |
191 | |
192 while (z_input[m + 1] < z_output[k]) { | |
193 ++m; | |
194 } | |
195 | |
196 result[k] = m; | |
197 } | |
198 | |
199 return result; | |
200 } | |
201 | |
202 void ColumnInterpolation::init_interpolation() { | |
203 | |
204 // coarse -> fine | |
205 m_coarse2fine = init_interpolation_indexes(m_z_coarse, m_z_fine); | |
206 | |
207 // fine -> coarse | |
208 m_fine2coarse = init_interpolation_indexes(m_z_fine, m_z_coarse); | |
209 | |
210 // decide if we're going to use linear or quadratic interpolation | |
211 double dz_min = m_z_coarse.back(); | |
212 double dz_max = 0.0; | |
213 for (unsigned int k = 0; k < Mz_coarse() - 1; ++k) { | |
214 const double dz = m_z_coarse[k + 1] - m_z_coarse[k]; | |
215 dz_min = std::min(dz, dz_min); | |
216 dz_max = std::max(dz, dz_max); | |
217 } | |
218 | |
219 if (fabs(dz_max - dz_min) <= 1.0e-8) { | |
220 m_use_linear_interpolation = true; | |
221 } else { | |
222 m_use_linear_interpolation = false; | |
223 } | |
224 | |
225 // initialize quadratic interpolation constants | |
226 if (not m_use_linear_interpolation) { | |
227 const unsigned int N = Mz_coarse() - 2; | |
228 m_constants.resize(3 * N); | |
229 for (unsigned int m = 0; m < N; ++m) { | |
230 const double | |
231 z0 = m_z_coarse[m], | |
232 z1 = m_z_coarse[m + 1], | |
233 z2 = m_z_coarse[m + 2]; | |
234 m_constants[3 * m + 0] = 1.0 / (z1 - z0); | |
235 m_constants[3 * m + 1] = 1.0 / (z2 - z0); | |
236 m_constants[3 * m + 2] = 1.0 / (z2 - z1); | |
237 } | |
238 } | |
239 | |
240 } | |
241 | |
242 } // end of namespace pism |