ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Classes | Public Types | Static Public Member Functions | List of all members
LinInterp Class Reference

#include <ERF_LinearInterpolate.H>

Classes

struct  InterpType
 

Public Types

enum  InterpMethod { extrap_method_zero = 0 , extrap_method_bndry = 1 , extrap_method_cycle = 2 }
 

Static Public Member Functions

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.)
 
static YAKL_INLINE void interp1d (const real1d &arrin, const int &n1, const real1d &arrout, const int &m1, const InterpType &interp_wgts)
 
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)
 
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)
 
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)
 

Member Enumeration Documentation

◆ InterpMethod

enum LinInterp::InterpMethod
Enumerator
extrap_method_zero 
extrap_method_bndry 
extrap_method_cycle 
14  {
15  extrap_method_zero = 0,
16  extrap_method_bndry = 1,
17  extrap_method_cycle = 2,
18  };
LinInterp::extrap_method_zero
@ extrap_method_zero
Definition: ERF_LinearInterpolate.H:15
LinInterp::extrap_method_bndry
@ extrap_method_bndry
Definition: ERF_LinearInterpolate.H:16
LinInterp::extrap_method_cycle
@ extrap_method_cycle
Definition: ERF_LinearInterpolate.H:17

Member Function Documentation

◆ init()

static YAKL_INLINE void LinInterp::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. 
)
inlinestatic
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  }

Referenced by CloudRadProps::get_mu_lambda_weights(), CloudRadProps::mitchell_ice_optics_lw(), and CloudRadProps::mitchell_ice_optics_sw().

Here is the caller graph for this function:

◆ interp1d()

static YAKL_INLINE void LinInterp::interp1d ( const real1d &  arrin,
const int &  n1,
const real1d &  arrout,
const int &  m1,
const InterpType &  interp_wgts 
)
inlinestatic
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  }

Referenced by CloudRadProps::get_mu_lambda_weights(), CloudRadProps::mitchell_ice_optics_lw(), and CloudRadProps::mitchell_ice_optics_sw().

Here is the caller graph for this function:

◆ interp2d1d()

static YAKL_INLINE void LinInterp::interp2d1d ( const real2d &  arrin,
const int &  n1,
const int &  n2,
const real1d &  arrout,
const int &  m1,
const InterpType &  wgt1,
const InterpType &  wgt2 
)
inlinestatic
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  }

Referenced by CloudRadProps::gam_liquid_lw(), and CloudRadProps::gam_liquid_sw().

Here is the caller graph for this function:

◆ interp2d2d()

static YAKL_INLINE void LinInterp::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 
)
inlinestatic
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  }

◆ interp3d2d()

static YAKL_INLINE void LinInterp::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 
)
inlinestatic
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  }
The documentation for this class was generated from the following file: