PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
ColumnInterpolation.cc
Go to the documentation of this file.
1 /* Copyright (C) 2014, 2015, 2021, 2022, 2023 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 the
6  * terms of the GNU General Public License as published by the Free Software
7  * Foundation; either version 3 of the License, or (at your option) any later
8  * version.
9  *
10  * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
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 USA
18  */
19 
20 #include "pism/util/ColumnInterpolation.hh"
21 #include <algorithm> // for max, min
22 #include <cmath> // for fabs
23 
24 namespace pism {
25 
26 ColumnInterpolation::ColumnInterpolation(const std::vector<double> &new_z_coarse,
27  const std::vector<double> &new_z_fine)
28  : m_z_fine(new_z_fine), m_z_coarse(new_z_coarse) {
30 }
31 
32 std::vector<double> ColumnInterpolation::coarse_to_fine(const std::vector<double> &input,
33  unsigned int k_max_result) const {
34  std::vector<double> result(Mz_fine());
35  coarse_to_fine(input.data(), k_max_result, result.data());
36  return result;
37 }
38 
39 void ColumnInterpolation::coarse_to_fine(const double *input,
40  unsigned int k_max_result,
41  double *result) const {
43  coarse_to_fine_linear(input, k_max_result, result);
44  } else {
45  coarse_to_fine_quadratic(input, k_max_result, result);
46  }
47 }
48 
49 void ColumnInterpolation::coarse_to_fine_linear(const double *input, unsigned int k_max_result,
50  double *result) const {
51  const unsigned int Mzfine = Mz_fine();
52  const unsigned int Mzcoarse = Mz_coarse();
53 
54  for (unsigned int k = 0; k < Mzfine; ++k) {
55  if (k > k_max_result) {
56  result[k] = input[m_coarse2fine[k]];
57  continue;
58  }
59 
60  unsigned int m = m_coarse2fine[k];
61 
62  // extrapolate (if necessary):
63  if (m == Mzcoarse - 1) {
64  result[k] = input[Mzcoarse - 1];
65  continue;
66  }
67 
68  const double incr = (m_z_fine[k] - m_z_coarse[m]) / (m_z_coarse[m + 1] - m_z_coarse[m]);
69  result[k] = input[m] + incr * (input[m + 1] - input[m]);
70  }
71 }
72 
73 void ColumnInterpolation::coarse_to_fine_quadratic(const double *input, unsigned int k_max_result,
74  double *result) const {
75  unsigned int k = 0, m = 0;
76  const unsigned int Mz = Mz_coarse();
77  for (m = 0; m < Mz - 2 and k <= k_max_result; ++m) {
78 
79  const double
80  z0 = m_z_coarse[m],
81  z1 = m_z_coarse[m + 1],
82  dz_inv = m_constants[3 * m + 0], // = 1.0 / (z1 - z0)
83  dz1_inv = m_constants[3 * m + 1], // = 1.0 / (z2 - z0)
84  dz2_inv = m_constants[3 * m + 2], // = 1.0 / (z2 - z1)
85  f0 = input[m],
86  f1 = input[m + 1],
87  f2 = input[m + 2];
88 
89  const double
90  d1 = (f1 - f0) * dz_inv,
91  d2 = (f2 - f0) * dz1_inv,
92  b = (d2 - d1) * dz2_inv,
93  a = d1 - b * (z1 - z0),
94  c = f0;
95 
96  for (; m_z_fine[k] < z1 and k <= k_max_result; ++k) {
97  const double s = m_z_fine[k] - z0;
98 
99  result[k] = s * (a + b * s) + c;
100  }
101  } // m-loop
102 
103  // check if we got to the end of the m-loop and use linear
104  // interpolation between the remaining 2 coarse levels
105  if (m == Mz - 2) {
106  const double
107  z0 = m_z_coarse[m],
108  z1 = m_z_coarse[m + 1],
109  f0 = input[m],
110  f1 = input[m + 1],
111  lambda = (f1 - f0) / (z1 - z0);
112 
113  for (; m_z_fine[k] < z1; ++k) {
114  result[k] = f0 + lambda * (m_z_fine[k] - z0);
115  }
116  }
117 
118  // fill the rest using constant extrapolation
119  const double f0 = input[Mz - 1];
120  for (; k <= k_max_result; ++k) {
121  result[k] = f0;
122  }
123 }
124 
125 std::vector<double> ColumnInterpolation::fine_to_coarse(const std::vector<double> &input) const {
126  std::vector<double> result(Mz_coarse());
127  fine_to_coarse(input.data(), result.data());
128  return result;
129 }
130 
131 void ColumnInterpolation::fine_to_coarse(const double *input, double *result) const {
132  const unsigned int N = Mz_coarse();
133 
134  for (unsigned int k = 0; k < N - 1; ++k) {
135  const int m = m_fine2coarse[k];
136 
137  const double increment = (m_z_coarse[k] - m_z_fine[m]) / (m_z_fine[m + 1] - m_z_fine[m]);
138  result[k] = input[m] + increment * (input[m + 1] - input[m]);
139  }
140 
141  result[N - 1] = input[m_fine2coarse[N - 1]];
142 }
143 
144 unsigned int ColumnInterpolation::Mz_coarse() const {
145  return m_z_coarse.size();
146 }
147 
148 unsigned int ColumnInterpolation::Mz_fine() const {
149  return m_z_fine.size();
150 }
151 
153  return m_z_fine[1] - m_z_fine[0];
154 }
155 
156 const std::vector<double>& ColumnInterpolation::z_fine() const {
157  return m_z_fine;
158 }
159 
160 const std::vector<double>& ColumnInterpolation::z_coarse() const {
161  return m_z_coarse;
162 }
163 
164 /*!
165  * Given two 1D grids, `z_input` and `z_output`, for each `z_output`
166  * index `k` we find an index `m` so that
167  *
168  * `z_input[m] < z_output[k] <= z_input[m+1]`
169  *
170  * In other words, we look for two consecutive points in the input
171  * grid that bracket a point in the output grid.
172  *
173  * This function sets `result[k] = m`. This information is then used
174  * to interpolate from the grid defined by `z_input` to the one
175  * defined by `z_output`.
176  *
177  * We use constant extrapolation outside the range defined by `z_input`.
178  */
179 static std::vector<unsigned int> init_interpolation_indexes(const std::vector<double>& z_input,
180  const std::vector<double>& z_output) {
181  std::vector<unsigned int> result(z_output.size());
182 
183  unsigned int m = 0;
184  for (unsigned int k = 0; k < z_output.size(); ++k) {
185 
186  if (z_output[k] <= z_input.front()) {
187  result[k] = 0;
188  continue;
189  }
190 
191  if (z_output[k] >= z_input.back()) {
192  result[k] = z_input.size() - 1;
193  continue;
194  }
195 
196  while (z_input[m + 1] < z_output[k]) {
197  ++m;
198  }
199 
200  result[k] = m;
201  }
202 
203  return result;
204 }
205 
207 
208  // coarse -> fine
210 
211  // fine -> coarse
213 
214  // decide if we're going to use linear or quadratic interpolation
215  double dz_min = m_z_coarse.back();
216  double dz_max = 0.0;
217  for (unsigned int k = 0; k < Mz_coarse() - 1; ++k) {
218  const double dz = m_z_coarse[k + 1] - m_z_coarse[k];
219  dz_min = std::min(dz, dz_min);
220  dz_max = std::max(dz, dz_max);
221  }
222 
223  const double eps = 1.0e-8;
224  m_use_linear_interpolation = (fabs(dz_max - dz_min) <= eps);
225 
226  // initialize quadratic interpolation constants
227  if (not m_use_linear_interpolation) {
228  const unsigned int N = Mz_coarse() - 2;
229  m_constants.resize(3 * N);
230  for (unsigned int m = 0; m < N; ++m) {
231  const double
232  z0 = m_z_coarse[m],
233  z1 = m_z_coarse[m + 1],
234  z2 = m_z_coarse[m + 2];
235  m_constants[3 * m + 0] = 1.0 / (z1 - z0);
236  m_constants[3 * m + 1] = 1.0 / (z2 - z0);
237  m_constants[3 * m + 2] = 1.0 / (z2 - z1);
238  }
239  }
240 
241 }
242 
243 } // end of namespace pism
void coarse_to_fine_quadratic(const double *input, unsigned int ks, double *result) const
std::vector< double > m_z_fine
std::vector< double > m_z_coarse
unsigned int Mz_fine() const
const std::vector< double > & z_coarse() const
void coarse_to_fine_linear(const double *input, unsigned int ks, double *result) const
std::vector< unsigned int > m_coarse2fine
std::vector< double > m_constants
unsigned int Mz_coarse() const
std::vector< unsigned int > m_fine2coarse
void coarse_to_fine(const double *input, unsigned int ks, double *result) const
ColumnInterpolation(const std::vector< double > &z_coarse, const std::vector< double > &z_fine)
const std::vector< double > & z_fine() const
void fine_to_coarse(const double *input, double *result) const
static const double z0
Definition: exactTestL.cc:42
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:193
static std::vector< unsigned int > init_interpolation_indexes(const std::vector< double > &z_input, const std::vector< double > &z_output)
static const double k
Definition: exactTestP.cc:42
const double d2
Definition: ssafd_code.cc:1