ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_Interpolation_WENO.H
Go to the documentation of this file.
1 #ifndef ERF_INTERPOLATE_WENO_H_
2 #define ERF_INTERPOLATE_WENO_H_
3 
4 #include "ERF_DataStruct.H"
5 
6 /**
7  * Interpolation operators used for WENO-5 scheme
8  */
9 struct WENO3
10 {
11  WENO3 (const amrex::Array4<const amrex::Real>& phi)
12  : m_phi(phi) {}
13 
14  AMREX_GPU_DEVICE
15  AMREX_FORCE_INLINE
16  void
17  InterpolateInX (const int& i,
18  const int& j,
19  const int& k,
20  const int& qty_index,
21  amrex::Real& val_lo,
22  amrex::Real upw_lo,
23  const amrex::Real /*upw_frac*/) const
24  {
25  // Data to interpolate on
26  amrex::Real sp1 = m_phi(i+1, j , k , qty_index);
27  amrex::Real s = m_phi(i , j , k , qty_index);
28  amrex::Real sm1 = m_phi(i-1, j , k , qty_index);
29  amrex::Real sm2 = m_phi(i-2, j , k , qty_index);
30 
31  if (upw_lo > tol) {
32  val_lo = Evaluate(sm2,sm1,s );
33  } else if (upw_lo < -tol) {
34  val_lo = Evaluate(sp1,s ,sm1);
35  } else {
36  val_lo = 0.5 * (s + sm1);
37  }
38  }
39 
40  AMREX_GPU_DEVICE
41  AMREX_FORCE_INLINE
42  void
43  InterpolateInY (const int& i,
44  const int& j,
45  const int& k,
46  const int& qty_index,
47  amrex::Real& val_lo,
48  amrex::Real upw_lo,
49  const amrex::Real /*upw_frac*/) const
50  {
51  // Data to interpolate on
52  amrex::Real sp1 = m_phi(i , j+1, k , qty_index);
53  amrex::Real s = m_phi(i , j , k , qty_index);
54  amrex::Real sm1 = m_phi(i , j-1, k , qty_index);
55  amrex::Real sm2 = m_phi(i , j-2, k , qty_index);
56 
57  if (upw_lo > tol) {
58  val_lo = Evaluate(sm2,sm1,s );
59  } else if (upw_lo < -tol) {
60  val_lo = Evaluate(sp1,s ,sm1);
61  } else {
62  val_lo = 0.5 * (s + sm1);
63  }
64  }
65 
66  AMREX_GPU_DEVICE
67  AMREX_FORCE_INLINE
68  void
69  InterpolateInZ (const int& i,
70  const int& j,
71  const int& k,
72  const int& qty_index,
73  amrex::Real& val_lo,
74  amrex::Real upw_lo,
75  const amrex::Real /*upw_frac*/) const
76  {
77  // Data to interpolate on
78  amrex::Real sp1 = m_phi(i , j , k+1, qty_index);
79  amrex::Real s = m_phi(i , j , k , qty_index);
80  amrex::Real sm1 = m_phi(i , j , k-1, qty_index);
81  amrex::Real sm2 = m_phi(i , j , k-2, qty_index);
82 
83  if (upw_lo > tol) {
84  val_lo = Evaluate(sm2,sm1,s );
85  } else if (upw_lo < -tol) {
86  val_lo = Evaluate(sp1,s ,sm1);
87  } else {
88  val_lo = 0.5 * (s + sm1);
89  }
90  }
91 
92  AMREX_GPU_DEVICE
93  AMREX_FORCE_INLINE
94  amrex::Real
95  Evaluate (const amrex::Real& sm1,
96  const amrex::Real& s ,
97  const amrex::Real& sp1) const
98  {
99  // Smoothing factors
100  amrex::Real b1 = (s - sm1) * (s - sm1);
101  amrex::Real b2 = (sp1 - s) * (sp1 - s);
102 
103  // Weight factors
104  amrex::Real w1 = g1 / ( (eps + b1) * (eps + b1) );
105  amrex::Real w2 = g2 / ( (eps + b2) * (eps + b2) );
106 
107  // Weight factor norm
108  amrex::Real wsum = w1 + w2;
109 
110  // Taylor expansions
111  amrex::Real v1 = -sm1 + 3.0 * s;
112  amrex::Real v2 = s + sp1;
113 
114  // Interpolated value
115  return ( (w1 * v1 + w2 * v2) / (2.0 * wsum) );
116  }
117 
118 private:
119  amrex::Array4<const amrex::Real> m_phi; // Quantity to interpolate
120  const amrex::Real eps=1.0e-6;
121  const amrex::Real tol=1.0e-12;
122  static constexpr amrex::Real g1=(1.0/3.0);
123  static constexpr amrex::Real g2=(2.0/3.0);
124 };
125 
126 /**
127  * Interpolation operators used for WENO-5 scheme
128  */
129 struct WENO5
130 {
131  WENO5 (const amrex::Array4<const amrex::Real>& phi)
132  : m_phi(phi) {}
133 
134  AMREX_GPU_DEVICE
135  AMREX_FORCE_INLINE
136  void
137  InterpolateInX (const int& i,
138  const int& j,
139  const int& k,
140  const int& qty_index,
141  amrex::Real& val_lo,
142  amrex::Real upw_lo,
143  const amrex::Real /*upw_frac*/) const
144  {
145  // Data to interpolate on
146  amrex::Real sp2 = m_phi(i+2, j , k , qty_index);
147  amrex::Real sp1 = m_phi(i+1, j , k , qty_index);
148  amrex::Real s = m_phi(i , j , k , qty_index);
149  amrex::Real sm1 = m_phi(i-1, j , k , qty_index);
150  amrex::Real sm2 = m_phi(i-2, j , k , qty_index);
151  amrex::Real sm3 = m_phi(i-3, j , k , qty_index);
152 
153  if (upw_lo > tol) {
154  val_lo = Evaluate(sm3,sm2,sm1,s ,sp1);
155  } else if (upw_lo < -tol) {
156  val_lo = Evaluate(sp2,sp1,s ,sm1,sm2);
157  } else {
158  val_lo = 0.5 * (s + sm1);
159  }
160  }
161 
162  AMREX_GPU_DEVICE
163  AMREX_FORCE_INLINE
164  void
165  InterpolateInY (const int& i,
166  const int& j,
167  const int& k,
168  const int& qty_index,
169  amrex::Real& val_lo,
170  amrex::Real upw_lo,
171  const amrex::Real /*upw_frac*/) const
172  {
173  // Data to interpolate on
174  amrex::Real sp2 = m_phi(i , j+2, k , qty_index);
175  amrex::Real sp1 = m_phi(i , j+1, k , qty_index);
176  amrex::Real s = m_phi(i , j , k , qty_index);
177  amrex::Real sm1 = m_phi(i , j-1, k , qty_index);
178  amrex::Real sm2 = m_phi(i , j-2, k , qty_index);
179  amrex::Real sm3 = m_phi(i , j-3, k , qty_index);
180 
181  if (upw_lo > tol) {
182  val_lo = Evaluate(sm3,sm2,sm1,s ,sp1);
183  } else if (upw_lo < -tol) {
184  val_lo = Evaluate(sp2,sp1,s ,sm1,sm2);
185  } else {
186  val_lo = 0.5 * (s + sm1);
187  }
188  }
189 
190  AMREX_GPU_DEVICE
191  AMREX_FORCE_INLINE
192  void
193  InterpolateInZ (const int& i,
194  const int& j,
195  const int& k,
196  const int& qty_index,
197  amrex::Real& val_lo,
198  amrex::Real upw_lo,
199  const amrex::Real /*upw_frac*/) const
200  {
201  // Data to interpolate on
202  amrex::Real sp2 = m_phi(i , j , k+2, qty_index);
203  amrex::Real sp1 = m_phi(i , j , k+1, qty_index);
204  amrex::Real s = m_phi(i , j , k , qty_index);
205  amrex::Real sm1 = m_phi(i , j , k-1, qty_index);
206  amrex::Real sm2 = m_phi(i , j , k-2, qty_index);
207  amrex::Real sm3 = m_phi(i , j , k-3, qty_index);
208 
209  if (upw_lo > tol) {
210  val_lo = Evaluate(sm3,sm2,sm1,s ,sp1);
211  } else if (upw_lo < -tol) {
212  val_lo = Evaluate(sp2,sp1,s ,sm1,sm2);
213  } else {
214  val_lo = 0.5 * (s + sm1);
215  }
216  }
217 
218  AMREX_GPU_DEVICE
219  AMREX_FORCE_INLINE
220  amrex::Real
221  Evaluate (const amrex::Real& sm2,
222  const amrex::Real& sm1,
223  const amrex::Real& s ,
224  const amrex::Real& sp1,
225  const amrex::Real& sp2) const
226  {
227  // Smoothing factors
228  amrex::Real b1 = c1 * (sm2 - 2.0 * sm1 + s) * (sm2 - 2.0 * sm1 + s) +
229  0.25 * (sm2 - 4.0 * sm1 + 3.0 * s) * (sm2 - 4.0 * sm1 + 3.0 * s);
230  amrex::Real b2 = c1 * (sm1 - 2.0 * s + sp1) * (sm1 - 2.0 * s + sp1) +
231  0.25 * (sm1 - sp1) * (sm1 - sp1);
232  amrex::Real b3 = c1 * (s - 2.0 * sp1 + sp2) * (s - 2.0 * sp1 + sp2) +
233  0.25 * (3.0 * s - 4.0 * sp1 + sp2) * (3.0 * s - 4.0 * sp1 + sp2);
234 
235  // Weight factors
236  amrex::Real w1 = g1 / ( (eps + b1) * (eps + b1) );
237  amrex::Real w2 = g2 / ( (eps + b2) * (eps + b2) );
238  amrex::Real w3 = g3 / ( (eps + b3) * (eps + b3) );
239 
240  // Weight factor norm
241  amrex::Real wsum = w1 + w2 + w3;
242 
243  // Taylor expansions
244  amrex::Real v1 = 2.0 * sm2 - 7.0 * sm1 + 11.0 * s;
245  amrex::Real v2 = -sm1 + 5.0 * s + 2.0 * sp1;
246  amrex::Real v3 = 2.0 * s + 5.0 * sp1 - sp2;
247 
248  // Interpolated value
249  return ( (w1 * v1 + w2 * v2 + w3 * v3) / (6.0 * wsum) );
250  }
251 
252 private:
253  amrex::Array4<const amrex::Real> m_phi; // Quantity to interpolate
254  const amrex::Real eps=1.0e-6;
255  const amrex::Real tol=1.0e-12;
256  static constexpr amrex::Real c1=(13.0/12.0);
257  static constexpr amrex::Real g1=(1.0/10.0);
258  static constexpr amrex::Real g2=(3.0/5.0);
259  static constexpr amrex::Real g3=(3.0/10.0);
260 };
261 
262 /**
263  * Interpolation operators used for WENO-7 scheme
264  */
265 struct WENO7
266 {
267  WENO7 (const amrex::Array4<const amrex::Real>& phi)
268  : m_phi(phi) {}
269 
270  AMREX_GPU_DEVICE
271  AMREX_FORCE_INLINE
272  void
273  InterpolateInX (const int& i,
274  const int& j,
275  const int& k,
276  const int& qty_index,
277  amrex::Real& val_lo,
278  amrex::Real upw_lo,
279  const amrex::Real /*upw_frac*/) const
280  {
281  // Data to interpolate on
282  amrex::Real sp3 = m_phi(i+3, j , k , qty_index);
283  amrex::Real sp2 = m_phi(i+2, j , k , qty_index);
284  amrex::Real sp1 = m_phi(i+1, j , k , qty_index);
285  amrex::Real s = m_phi(i , j , k , qty_index);
286  amrex::Real sm1 = m_phi(i-1, j , k , qty_index);
287  amrex::Real sm2 = m_phi(i-2, j , k , qty_index);
288  amrex::Real sm3 = m_phi(i-3, j , k , qty_index);
289  amrex::Real sm4 = m_phi(i-4, j , k , qty_index);
290 
291  if (upw_lo > tol) {
292  val_lo = Evaluate(sm4,sm3,sm2,sm1,s ,sp1,sp2);
293  } else if (upw_lo < -tol) {
294  val_lo = Evaluate(sp3,sp2,sp1,s ,sm1,sm2,sm3);
295  } else {
296  val_lo = 0.5 * (s + sm1);
297  }
298  }
299 
300  AMREX_GPU_DEVICE
301  AMREX_FORCE_INLINE
302  void
303  InterpolateInY (const int& i,
304  const int& j,
305  const int& k,
306  const int& qty_index,
307  amrex::Real& val_lo,
308  amrex::Real upw_lo,
309  const amrex::Real /*upw_frac*/) const
310  {
311  // Data to interpolate on
312  amrex::Real sp3 = m_phi(i , j+3, k , qty_index);
313  amrex::Real sp2 = m_phi(i , j+2, k , qty_index);
314  amrex::Real sp1 = m_phi(i , j+1, k , qty_index);
315  amrex::Real s = m_phi(i , j , k , qty_index);
316  amrex::Real sm1 = m_phi(i , j-1, k , qty_index);
317  amrex::Real sm2 = m_phi(i , j-2, k , qty_index);
318  amrex::Real sm3 = m_phi(i , j-3, k , qty_index);
319  amrex::Real sm4 = m_phi(i , j-4, k , qty_index);
320 
321  if (upw_lo > tol) {
322  val_lo = Evaluate(sm4,sm3,sm2,sm1,s ,sp1,sp2);
323  } else if (upw_lo < -tol) {
324  val_lo = Evaluate(sp3,sp2,sp1,s ,sm1,sm2,sm3);
325  } else {
326  val_lo = 0.5 * (s + sm1);
327  }
328  }
329 
330  AMREX_GPU_DEVICE
331  AMREX_FORCE_INLINE
332  void
333  InterpolateInZ (const int& i,
334  const int& j,
335  const int& k,
336  const int& qty_index,
337  amrex::Real& val_lo,
338  amrex::Real upw_lo,
339  const amrex::Real /*upw_frac*/) const
340  {
341  // Data to interpolate on
342  amrex::Real sp3 = m_phi(i , j , k+2, qty_index);
343  amrex::Real sp2 = m_phi(i , j , k+2, qty_index);
344  amrex::Real sp1 = m_phi(i , j , k+1, qty_index);
345  amrex::Real s = m_phi(i , j , k , qty_index);
346  amrex::Real sm1 = m_phi(i , j , k-1, qty_index);
347  amrex::Real sm2 = m_phi(i , j , k-2, qty_index);
348  amrex::Real sm3 = m_phi(i , j , k-3, qty_index);
349  amrex::Real sm4 = m_phi(i , j , k-3, qty_index);
350 
351  if (upw_lo > tol) {
352  val_lo = Evaluate(sm4,sm3,sm2,sm1,s ,sp1,sp2);
353  } else if (upw_lo < -tol) {
354  val_lo = Evaluate(sp3,sp2,sp1,s ,sm1,sm2,sm3);
355  } else {
356  val_lo = 0.5 * (s + sm1);
357  }
358  }
359 
360  AMREX_GPU_DEVICE
361  AMREX_FORCE_INLINE
362  amrex::Real
363  Evaluate (const amrex::Real& sm3,
364  const amrex::Real& sm2,
365  const amrex::Real& sm1,
366  const amrex::Real& s ,
367  const amrex::Real& sp1,
368  const amrex::Real& sp2,
369  const amrex::Real& sp3) const
370  {
371  // Smoothing factors
372  amrex::Real b1 = ( sm3 * sm3 * 6649./2880.0
373  - sm3 * sm2 * 2623./160.0
374  + sm3 * sm1 * 9449./480.0
375  - sm3 * s * 11389./1440.0
376  + sm2 * sm2 * 28547./960.0
377  - sm2 * sm1 * 35047./480.0
378  + sm2 * s * 14369./480.0
379  + sm1 * sm1 * 44747./960.0
380  - sm1 * s * 6383./160.0
381  + s * s * 25729./2880.0 );
382  amrex::Real b2 = ( sm2 * sm2 * 3169/2880.0
383  - sm2 * sm1 * 3229/480.0
384  + sm2 * s * 3169/480.0
385  - sm2 * sp1 * 2989/1440.0
386  + sm1 * sm1 * 11147/960.0
387  - sm1 * s * 11767/480.0
388  + sm1 * sp1 * 1283/160.0
389  + s * s * 13667/960.0
390  - s * sp1 * 5069/480.0
391  + sp1 * sp1 * 6649/2880.0 );
392  amrex::Real b3 = ( sm1 * sm1 * 6649./2880.0
393  - sm1 * s * 5069./480.0
394  + sm1 * sp1 * 1283./160.0
395  - sm1 * sp2 * 2989./1440.0
396  + s * s * 13667./960.0
397  - s * sp1 * 11767./480.0
398  + s * sp2 * 3169./480.0
399  + sp1 * sp1 * 11147./960.0
400  - sp1 * sp2 * 3229./480.0
401  + sp2 * sp2 * 3169./2880.0 );
402  amrex::Real b4 = ( s * s * 25729./2880.
403  - s * sp1 * 6383./160.
404  + s * sp2 * 14369./480.
405  - s * sp3 * 11389./1440.
406  + sp1 * sp1 * 44747./960.
407  - sp1 * sp2 * 35047./480.
408  + sp1 * sp3 * 9449./480.
409  + sp2 * sp2 * 28547./960.
410  - sp2 * sp3 * 2623./160.
411  + sp3 * sp3 * 6649./2880. );
412 
413  // Weight factors
414  amrex::Real w1 = g1 / ( (eps + b1) * (eps + b1) );
415  amrex::Real w2 = g2 / ( (eps + b2) * (eps + b2) );
416  amrex::Real w3 = g3 / ( (eps + b3) * (eps + b3) );
417  amrex::Real w4 = g4 / ( (eps + b4) * (eps + b4) );
418 
419  // Weight factor norm
420  amrex::Real wsum = w1 + w2 + w3 + w4;
421 
422  // Taylor expansions
423  amrex::Real v1 = (-0.3125)*sm3 + ( 1.3125)*sm2 + (-2.1875)*sm1 + ( 2.1875)*s;
424  amrex::Real v2 = ( 0.0625)*sm2 + (-0.3125)*sm1 + ( 0.9375)*s + ( 0.3125)*sp1;
425  amrex::Real v3 = (-0.0625)*sm1 + ( 0.5625)*s + ( 0.5625)*sp1 + (-0.0625)*sp2;
426  amrex::Real v4 = ( 0.3125)*s + ( 0.9375)*sp1 + (-0.3125)*sp2 + ( 0.0625)*sp3;
427 
428  // Interpolated value
429  return ( (w1 * v1 + w2 * v2 + w3 * v3 + w4 * v4) / (wsum) );
430  }
431 
432 private:
433  amrex::Array4<const amrex::Real> m_phi; // Quantity to interpolate
434  const amrex::Real eps=1.0e-8;
435  const amrex::Real tol=1.0e-12;
436  static constexpr amrex::Real g1=( 1.0/64.0);
437  static constexpr amrex::Real g2=(21.0/64.0);
438  static constexpr amrex::Real g3=(35.0/64.0);
439  static constexpr amrex::Real g4=( 7.0/64.0);
440 };
441 #endif
Definition: ERF_Interpolation_WENO.H:10
static constexpr amrex::Real g1
Definition: ERF_Interpolation_WENO.H:122
amrex::Array4< const amrex::Real > m_phi
Definition: ERF_Interpolation_WENO.H:119
const amrex::Real eps
Definition: ERF_Interpolation_WENO.H:120
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Evaluate(const amrex::Real &sm1, const amrex::Real &s, const amrex::Real &sp1) const
Definition: ERF_Interpolation_WENO.H:95
static constexpr amrex::Real g2
Definition: ERF_Interpolation_WENO.H:123
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInY(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:43
WENO3(const amrex::Array4< const amrex::Real > &phi)
Definition: ERF_Interpolation_WENO.H:11
const amrex::Real tol
Definition: ERF_Interpolation_WENO.H:121
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInX(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:17
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInZ(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:69
Definition: ERF_Interpolation_WENO.H:130
const amrex::Real eps
Definition: ERF_Interpolation_WENO.H:254
amrex::Array4< const amrex::Real > m_phi
Definition: ERF_Interpolation_WENO.H:253
static constexpr amrex::Real c1
Definition: ERF_Interpolation_WENO.H:256
static constexpr amrex::Real g2
Definition: ERF_Interpolation_WENO.H:258
WENO5(const amrex::Array4< const amrex::Real > &phi)
Definition: ERF_Interpolation_WENO.H:131
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInX(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:137
const amrex::Real tol
Definition: ERF_Interpolation_WENO.H:255
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInZ(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:193
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInY(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:165
static constexpr amrex::Real g3
Definition: ERF_Interpolation_WENO.H:259
static constexpr amrex::Real g1
Definition: ERF_Interpolation_WENO.H:257
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Evaluate(const amrex::Real &sm2, const amrex::Real &sm1, const amrex::Real &s, const amrex::Real &sp1, const amrex::Real &sp2) const
Definition: ERF_Interpolation_WENO.H:221
Definition: ERF_Interpolation_WENO.H:266
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInY(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:303
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Evaluate(const amrex::Real &sm3, const amrex::Real &sm2, const amrex::Real &sm1, const amrex::Real &s, const amrex::Real &sp1, const amrex::Real &sp2, const amrex::Real &sp3) const
Definition: ERF_Interpolation_WENO.H:363
WENO7(const amrex::Array4< const amrex::Real > &phi)
Definition: ERF_Interpolation_WENO.H:267
const amrex::Real eps
Definition: ERF_Interpolation_WENO.H:434
static constexpr amrex::Real g4
Definition: ERF_Interpolation_WENO.H:439
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInX(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:273
static constexpr amrex::Real g3
Definition: ERF_Interpolation_WENO.H:438
static constexpr amrex::Real g2
Definition: ERF_Interpolation_WENO.H:437
static constexpr amrex::Real g1
Definition: ERF_Interpolation_WENO.H:436
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void InterpolateInZ(const int &i, const int &j, const int &k, const int &qty_index, amrex::Real &val_lo, amrex::Real upw_lo, const amrex::Real) const
Definition: ERF_Interpolation_WENO.H:333
amrex::Array4< const amrex::Real > m_phi
Definition: ERF_Interpolation_WENO.H:433
const amrex::Real tol
Definition: ERF_Interpolation_WENO.H:435