ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_Linear_interpolate.H
Go to the documentation of this file.
1 //
2 // interpolation of data in latitude, longitude, and time
3 //
4 #ifndef ERF_LINEAR_INTERP_H_
5 #define ERF_LINEAR_INTERP_H_
6 
7 #include "rrtmgp_const.h"
8 using yakl::fortran::parallel_for;
9 using yakl::fortran::SimpleBounds;
10 
11 class LinInterp
12 {
13  public:
18  };
19 
20  struct InterpType {
21  real1d wgts;
22  real1d wgtn;
23  int1d jjm;
24  int1d jjp;
25  };
26 
27  // Initialize a variable of type(interp_type) with weights for linear interpolation.
28  // this variable can then be used in calls to lininterp1d and lininterp2d.
29  // yin is a 1d array of length n_in of locations to interpolate from - this array must
30  // be monotonic but can be increasing or decreasing
31  // yout is a 1d array of length nout of locations to interpolate to, this array need
32  // not be ordered
33  // extrap_method determines how to handle yout points beyond the bounds of yin
34  // if 0 set values outside output grid to 0
35  // if 1 set to boundary value
36  // if 2 set to cyclic boundaries
37  // optional values cyclicmin and cyclicmax can be used to set the bounds of the
38  // cyclic mapping - these default to 0 and 360.
39  YAKL_INLINE
40  static void init (const real1d& yin,
41  const int& n_in,
42  const real1d& yout,
43  const int& nout,
44  const InterpMethod& extrap_method,
45  InterpType& interp_wgts,
46  real cyclicmin = 0.,
47  real cyclicmax= 0.)
48  {
49  real cmin, cmax;
50  real dyinwrap;
51  real avgdyin;
52  //
53  // Check validity of input coordinate arrays: must be monotonically increasing,
54  // and have a total of at least 2 elements
55  //
56  if(cyclicmin != 0.) {
57  cmin=cyclicmin;
58  } else {
59  cmin=0;
60  }
61  if(cyclicmax != 0.) {
62  cmax=cyclicmax;
63  } else {
64  cmax=360.0;
65  }
66 
67  bool increasing = true;
68 
69  interp_wgts.jjm = int1d("jjm", nout);
70  interp_wgts.jjp = int1d("jjp", nout);
71  interp_wgts.wgts = real1d("wgts", nout);
72  interp_wgts.wgtn = real1d("wgtn", nout);
73 
74  parallel_for(SimpleBounds<1>(n_in-1), YAKL_LAMBDA (int j)
75  {
76  if(yin(j) > yin(j+1)) {
77  printf("init: inputs are not monotonic!\n");
78  return;
79  }
80  });
81 
82  //
83  // Initialize index arrays for later checking
84  //
85  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
86  {
87  interp_wgts.jjm(j) = 0;
88  interp_wgts.jjp(j) = 0;
89  });
90 
91  switch (extrap_method) {
92  case extrap_method_zero:
93  //
94  // For values which extend beyond boundaries, set weights
95  // such that values will be 0.
96  //
97  if(increasing) {
98  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
99  {
100  if (yout(j) < yin(1)) {
101  interp_wgts.jjm(j) = 1;
102  interp_wgts.jjp(j) = 1;
103  interp_wgts.wgts(j) = 0.;
104  interp_wgts.wgtn(j) = 0.;
105  }
106  else if (yout(j) > yin(n_in)) {
107  interp_wgts.jjm(j) = n_in;
108  interp_wgts.jjp(j) = n_in;
109  interp_wgts.wgts(j) = 0.;
110  interp_wgts.wgtn(j) = 0.;
111  }
112  });
113  }
114  else {
115  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
116  {
117  if (yout(j) > yin(1)) {
118  interp_wgts.jjm(j) = 1;
119  interp_wgts.jjp(j) = 1;
120  interp_wgts.wgts(j) = 0.;
121  interp_wgts.wgtn(j) = 0.;
122  }
123  else if (yout(j) < yin(n_in)) {
124  interp_wgts.jjm(j) = n_in;
125  interp_wgts.jjp(j) = n_in;
126  interp_wgts.wgts(j) = 0.;
127  interp_wgts.wgtn(j) = 0.;
128  }
129  });
130  }
131  break;
132  case extrap_method_bndry:
133  //
134  // For values which extend beyond boundaries, set weights
135  // such that values will just be copied.
136  //
137  if(increasing) {
138  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
139  {
140  if (yout(j) <= yin(1)) {
141  interp_wgts.jjm(j) = 1;
142  interp_wgts.jjp(j) = 1;
143  interp_wgts.wgts(j) = 1.;
144  interp_wgts.wgtn(j) = 0.;
145  }
146  else if (yout(j) > yin(n_in)) {
147  interp_wgts.jjm(j) = n_in;
148  interp_wgts.jjp(j) = n_in;
149  interp_wgts.wgts(j) = 1.;
150  interp_wgts.wgtn(j) = 0.;
151  }
152  });
153  }
154  else {
155  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
156  {
157  if (yout(j) > yin(1)) {
158  interp_wgts.jjm(j) = 1;
159  interp_wgts.jjp(j) = 1;
160  interp_wgts.wgts(j) = 1.;
161  interp_wgts.wgtn(j) = 0.;
162  }
163  else if (yout(j) <= yin(n_in)) {
164  interp_wgts.jjm(j) = n_in;
165  interp_wgts.jjp(j) = n_in;
166  interp_wgts.wgts(j) = 1.;
167  interp_wgts.wgtn(j) = 0.;
168  }
169  });
170  }
171  break;
172  case extrap_method_cycle:
173  //
174  // For values which extend beyond boundaries, set weights
175  // for circular boundaries
176  //
177  dyinwrap = yin(1) + (cmax-cmin) - yin(n_in);
178  avgdyin = abs(yin(n_in)-yin(1))/(n_in-1.);
179  auto ratio = dyinwrap/avgdyin;
180  if (ratio < 0.9 || ratio > 1.1) {
181  printf("ratio is too large, dyinwrap= %13.6e, avgdyin= %13.6e, yin(1) = %13.6e, yin(n_in)= %13.6e\n",
182  dyinwrap, avgdyin, yin(1), yin(n_in));
183  }
184 
185  if(increasing) {
186  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
187  {
188  if (yout(j) <= yin(1)) {
189  interp_wgts.jjm(j) = n_in;
190  interp_wgts.jjp(j) = 1;
191  interp_wgts.wgts(j) = (yin(1)-yout(j))/dyinwrap;
192  interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin) - yin(n_in))/dyinwrap;
193  }
194  else if (yout(j) > yin(n_in)) {
195  interp_wgts.jjm(j) = n_in;
196  interp_wgts.jjp(j) = 1;
197  interp_wgts.wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap;
198  interp_wgts.wgtn(j) = (yout(j)-yin(n_in))/dyinwrap;
199  }
200  });
201  }
202  else {
203  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
204  {
205  if (yout(j) > yin(1)) {
206  interp_wgts.jjm(j) = n_in;
207  interp_wgts.jjp(j) = 1;
208  interp_wgts.wgts(j) = (yin(1)-yout(j))/dyinwrap;
209  interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin) - yin(n_in))/dyinwrap;
210  }
211  else if (yout(j) <= yin(n_in)) {
212  interp_wgts.jjm(j) = n_in;
213  interp_wgts.jjp(j) = 1;
214  interp_wgts.wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap;
215  interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin)-yin(n_in))/dyinwrap;
216  }
217  });
218  }
219  break;
220  }
221  //
222  // Loop though output indices finding input indices and weights
223  //
224  if(increasing) {
225  parallel_for(SimpleBounds<2>(nout, n_in-1), YAKL_LAMBDA (int j, int jj)
226  {
227  if (yout(j) > yin(jj) && yout(j) <= yin(jj+1)) {
228  interp_wgts.jjm(j) = jj;
229  interp_wgts.jjp(j) = jj + 1;
230  interp_wgts.wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj));
231  interp_wgts.wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj));
232  return;
233  }
234  });
235  }
236  else {
237  parallel_for(SimpleBounds<2>(nout, n_in-1), YAKL_LAMBDA (int j, int jj)
238  {
239  if (yout(j) <= yin(jj) && yout(j) > yin(jj+1)) {
240  interp_wgts.jjm(j) = jj;
241  interp_wgts.jjp(j) = jj + 1;
242  interp_wgts.wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj));
243  interp_wgts.wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj));
244  return;
245  }
246  });
247  }
248 
249  //
250  // Check that interp/extrap points have been found for all outputs
251  //
252  parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j)
253  {
254  int count = 0;
255  amrex::ignore_unused(count);
256  if (interp_wgts.jjm(j) == 0 || interp_wgts.jjp(j) == 0) count += 1;
257  auto frac = interp_wgts.wgts(j) + interp_wgts.wgtn(j);
258  if ((frac < 0.9 || frac > 1.1) && extrap_method != 0)
259  printf("interpolation error!\n");
260  });
261 
262  }
263 
264  //
265  // Purpose: Do a linear interpolation from input mesh to output
266  // mesh with weights as set in lininterp_init.
267  YAKL_INLINE
268  static void interp1d (const real1d& arrin,
269  const int& n1,
270  const real1d& arrout,
271  const int& m1,
272  const InterpType& interp_wgts)
273  {
274  //
275  // Do the interpolation
276  //
277  parallel_for(SimpleBounds<1>(m1), YAKL_LAMBDA (int j)
278  {
279  arrout(j) = arrin(interp_wgts.jjm(j))*interp_wgts.wgts(j) +
280  arrin(interp_wgts.jjp(j))*interp_wgts.wgtn(j);
281  });
282  }
283 
284  YAKL_INLINE
285  static void interp2d2d (const real2d& arrin,
286  const int& n1,
287  const int& n2,
288  const real2d& arrout,
289  const int& m1,
290  const int& m2,
291  const InterpType& wgt1,
292  const InterpType& wgt2)
293  {
294 
295  real2d arrtmp("arrtmp",n1,m2);
296 
297  parallel_for(SimpleBounds<2>(n1, m2), YAKL_LAMBDA (int i, int j)
298  {
299  arrtmp(i,j) = arrin(i,wgt2.jjm(j))*wgt2.wgts(j) + arrin(i,wgt2.jjp(j))*wgt2.wgtn(j);
300  });
301 
302  parallel_for(SimpleBounds<2>(n1, m2), YAKL_LAMBDA (int i, int j)
303  {
304  arrout(i,j) = arrtmp(wgt1.jjm(i),j)*wgt1.wgts(i) + arrtmp(wgt1.jjp(i),j)*wgt1.wgtn(i);
305  });
306  }
307 
308  YAKL_INLINE
309  static void interp2d1d (const real2d& arrin,
310  const int& n1,
311  const int& n2,
312  const real1d& arrout,
313  const int& m1,
314  const InterpType& wgt1,
315  const InterpType& wgt2)
316  {
317  parallel_for(SimpleBounds<1>(m1), YAKL_LAMBDA (int i)
318  {
319  arrout(i) = arrin(wgt1.jjm(i),wgt2.jjm(i))*wgt1.wgts(i)*wgt2.wgts(i)+arrin(wgt1.jjp(i),wgt2.jjm(i))*wgt1.wgtn(i)*wgt2.wgts(i)
320  + arrin(wgt1.jjm(i),wgt2.jjp(i))*wgt1.wgtn(i)*wgt2.wgtn(i)+arrin(wgt1.jjp(i),wgt2.jjp(i))*wgt1.wgtn(i)*wgt2.wgtn(i);
321  });
322  }
323 
324  YAKL_INLINE
325  static void interp3d2d (const real3d& arrin,
326  const int& n1,
327  const int& n2,
328  const int& n3,
329  const real2d& arrout,
330  const int& m1,
331  const int& len1,
332  const InterpType& wgt1,
333  const InterpType& wgt2)
334  {
335  parallel_for(SimpleBounds<2>(m1, n3), YAKL_LAMBDA (int i, int k)
336  {
337  arrout(i,k) = arrin(wgt1.jjm(i),wgt2.jjm(i),k)*wgt1.wgts(i)*wgt2.wgts(i)+arrin(wgt1.jjp(i),wgt2.jjm(i),k)*wgt1.wgtn(i)*wgt2.wgts(i) +
338  arrin(wgt1.jjm(i),wgt2.jjp(i),k)*wgt1.wgts(i)*wgt2.wgtn(i)+arrin(wgt1.jjp(i),wgt2.jjp(i),k)*wgt1.wgtn(i)*wgt2.wgtn(i);
339  });
340  }
341 };
342 #endif
Definition: ERF_Linear_interpolate.H:12
static YAKL_INLINE void init(const real1d &yin, const int &n_in, const real1d &yout, const int &nout, const InterpMethod &extrap_method, InterpType &interp_wgts, real cyclicmin=0., real cyclicmax=0.)
Definition: ERF_Linear_interpolate.H:40
InterpMethod
Definition: ERF_Linear_interpolate.H:14
@ extrap_method_zero
Definition: ERF_Linear_interpolate.H:15
@ extrap_method_bndry
Definition: ERF_Linear_interpolate.H:16
@ extrap_method_cycle
Definition: ERF_Linear_interpolate.H:17
static YAKL_INLINE void interp3d2d(const real3d &arrin, const int &n1, const int &n2, const int &n3, const real2d &arrout, const int &m1, const int &len1, const InterpType &wgt1, const InterpType &wgt2)
Definition: ERF_Linear_interpolate.H:325
static YAKL_INLINE void interp2d2d(const real2d &arrin, const int &n1, const int &n2, const real2d &arrout, const int &m1, const int &m2, const InterpType &wgt1, const InterpType &wgt2)
Definition: ERF_Linear_interpolate.H:285
static YAKL_INLINE void interp1d(const real1d &arrin, const int &n1, const real1d &arrout, const int &m1, const InterpType &interp_wgts)
Definition: ERF_Linear_interpolate.H:268
static YAKL_INLINE void interp2d1d(const real2d &arrin, const int &n1, const int &n2, const real1d &arrout, const int &m1, const InterpType &wgt1, const InterpType &wgt2)
Definition: ERF_Linear_interpolate.H:309
Definition: ERF.H:92
real1d wgts
Definition: ERF_Linear_interpolate.H:21
int1d jjp
Definition: ERF_Linear_interpolate.H:24
real1d wgtn
Definition: ERF_Linear_interpolate.H:22
int1d jjm
Definition: ERF_Linear_interpolate.H:23