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