ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_Interpolation.H
Go to the documentation of this file.
1 #ifndef ERF_INTERPOLATE_H_
2 #define ERF_INTERPOLATE_H_
3 
4 #include "ERF_DataStruct.H"
5 #include "ERF_Interpolation_UPW.H"
6 #include "ERF_Interpolation_WENO.H"
7 #include "ERF_Interpolation_WENO_Z.H"
8 
9 /**
10  * Interpolation operators used in construction of advective fluxes using
11  * non-WENO centered/upwind schemes.
12  */
13 
14 AMREX_GPU_HOST_DEVICE
15 AMREX_FORCE_INLINE
16 bool
17 isSupportedInterpolationAdvType (const AdvType adv_type) noexcept
18 {
19  return adv_type == AdvType::Centered_2nd ||
20  adv_type == AdvType::Upwind_3rd ||
21  adv_type == AdvType::Centered_4th ||
22  adv_type == AdvType::Upwind_5th ||
23  adv_type == AdvType::Centered_6th;
24 }
25 
26 /**
27  * Compute interpolated value based on averages, differences, and advection type.
28  *
29  * @param avg1 First average
30  * @param avg2 Second average
31  * @param avg3 Third average
32  * @param diff1 First difference
33  * @param diff2 Second difference
34  * @param diff3 Third difference
35  * @param scaled_upw Scaled upwind factor
36  * @param adv_type Advection scheme type
37  * @return Interpolated value
38  */
39 AMREX_GPU_DEVICE
40 AMREX_FORCE_INLINE
41 amrex::Real
42 interpolatedVal (amrex::Real avg1, amrex::Real avg2, amrex::Real avg3,
43  amrex::Real diff1, amrex::Real diff2, amrex::Real diff3,
44  amrex::Real scaled_upw, const AdvType adv_type)
45 {
46  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
47  isSupportedInterpolationAdvType(adv_type),
48  "interpolatedVal only supports Centered_2nd, Upwind_3rd, Centered_4th, Upwind_5th, and Centered_6th");
49 
50  amrex::Real myInterpolatedVal(amrex::Real(0));
51  if (adv_type == AdvType::Centered_2nd) {
52  myInterpolatedVal = myhalf * avg1;
53  } else if (adv_type == AdvType::Upwind_3rd) {
54  myInterpolatedVal = (amrex::Real(7.0)/amrex::Real(12.0))*avg1 -(amrex::Real(1)/amrex::Real(12.0))*avg2 + (scaled_upw/amrex::Real(12.0))*(diff2 - amrex::Real(3)*diff1);
55  } else if (adv_type == AdvType::Centered_4th) {
56  myInterpolatedVal = (amrex::Real(7.0)/amrex::Real(12.0))*avg1 -(amrex::Real(1)/amrex::Real(12.0))*avg2;
57  } else if (adv_type == AdvType::Upwind_5th) {
58  myInterpolatedVal = (amrex::Real(37.0)/amrex::Real(60.0))*avg1 -(amrex::Real(2)/amrex::Real(15.0))*avg2 +(amrex::Real(1)/amrex::Real(60.0))*avg3
59  -(scaled_upw/amrex::Real(60.0))*(diff3 - amrex::Real(5.0)*diff2 + amrex::Real(10.0)*diff1);
60  } else if (adv_type == AdvType::Centered_6th) {
61  myInterpolatedVal = (amrex::Real(37.0)/amrex::Real(60.0))*avg1 -(amrex::Real(2)/amrex::Real(15.0))*avg2 +(amrex::Real(1)/amrex::Real(60.0))*avg3;
62  } else {
63  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
64  false,
65  "interpolatedVal only supports Centered_2nd, Upwind_3rd, Centered_4th, Upwind_5th, and Centered_6th");
66  }
67  return myInterpolatedVal;
68 }
69 
70 /**
71  * Interpolate a quantity in the x-direction.
72  *
73  * @param i x-index
74  * @param j y-index
75  * @param k z-index
76  * @param qty Quantity field
77  * @param qty_index Component index
78  * @param upw Upwind direction
79  * @param adv_type Advection scheme type
80  * @return Interpolated value in x
81  */
82 AMREX_GPU_DEVICE
83 AMREX_FORCE_INLINE
84 amrex::Real
85 InterpolateInX (int i, int j, int k, const amrex::Array4<const amrex::Real>& qty,
86  int qty_index, amrex::Real upw, const AdvType adv_type)
87 {
88 
89  if (adv_type == AdvType::Centered_2nd) {
90  return myhalf * (qty(i,j,k,qty_index) + qty(i-1,j,k,qty_index));
91  } else {
92 
93  amrex::Real avg1 = amrex::Real(0); amrex::Real avg2 = amrex::Real(0); amrex::Real avg3 = amrex::Real(0);
94  amrex::Real diff1 = amrex::Real(0); amrex::Real diff2 = amrex::Real(0); amrex::Real diff3 = amrex::Real(0);
95  amrex::Real scaled_upw = amrex::Real(0);
96  //
97  // The value that comes in has not been normalized so we do that here
98  if (upw != amrex::Real(0)) { scaled_upw = (upw > 0) ? amrex::Real(1) : -amrex::Real(1); }
99 
100  avg1 = (qty(i, j, k, qty_index) + qty(i-1, j, k, qty_index));
101  diff1 = (qty(i, j, k, qty_index) - qty(i-1, j, k, qty_index));
102  avg2 = (qty(i+1, j, k, qty_index) + qty(i-2, j, k, qty_index));
103  diff2 = (qty(i+1, j, k, qty_index) - qty(i-2, j, k, qty_index));
104  if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th)
105  {
106  avg3 = (qty(i+2, j, k, qty_index) + qty(i-3, j, k, qty_index));
107  diff3 = (qty(i+2, j, k, qty_index) - qty(i-3, j, k, qty_index));
108  }
109  return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type);
110  }
111 }
112 
113 /**
114  * Interpolate a quantity in the y-direction.
115  *
116  * @param i x-index
117  * @param j y-index
118  * @param k z-index
119  * @param qty Quantity field
120  * @param qty_index Component index
121  * @param upw Upwind direction
122  * @param adv_type Advection scheme type
123  * @return Interpolated value in y
124  */
125 AMREX_GPU_DEVICE
126 AMREX_FORCE_INLINE
127 amrex::Real
128 InterpolateInY (int i, int j, int k, const amrex::Array4<const amrex::Real>& qty,
129  int qty_index, amrex::Real upw, const AdvType adv_type)
130 {
131  if (adv_type == AdvType::Centered_2nd) {
132  return myhalf * (qty(i,j,k,qty_index) + qty(i,j-1,k,qty_index));
133  } else {
134 
135  amrex::Real avg1; amrex::Real avg2; amrex::Real avg3 = amrex::Real(0);
136  amrex::Real diff1; amrex::Real diff2; amrex::Real diff3 = amrex::Real(0);
137  amrex::Real scaled_upw = amrex::Real(0);
138 
139  // The value that comes in has not been normalized so we do that here
140  if (upw != amrex::Real(0)) { scaled_upw = (upw > 0) ? amrex::Real(1) : -amrex::Real(1); }
141 
142  avg1 = (qty(i, j , k, qty_index) + qty(i, j-1, k, qty_index));
143  diff1 = (qty(i, j , k, qty_index) - qty(i, j-1, k, qty_index));
144  avg2 = (qty(i, j+1, k, qty_index) + qty(i, j-2, k, qty_index));
145  diff2 = (qty(i, j+1, k, qty_index) - qty(i, j-2, k, qty_index));
146  if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th)
147  {
148  avg3 = (qty(i, j+2, k, qty_index) + qty(i, j-3, k, qty_index));
149  diff3 = (qty(i, j+2, k, qty_index) - qty(i, j-3, k, qty_index));
150  }
151  return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type);
152  }
153 }
154 
155 /**
156  * Interpolate a quantity in the z-direction.
157  *
158  * @param i x-index
159  * @param j y-index
160  * @param k z-index
161  * @param qty Quantity field
162  * @param qty_index Component index
163  * @param upw Upwind direction
164  * @param adv_type Advection scheme type
165  * @return Interpolated value in z
166  */
167 AMREX_GPU_DEVICE
168 AMREX_FORCE_INLINE
169 amrex::Real
170 InterpolateInZ (int i, int j, int k, const amrex::Array4<const amrex::Real>& qty,
171  int qty_index, amrex::Real upw, const AdvType adv_type)
172 {
173  if (adv_type == AdvType::Centered_2nd) {
174  return myhalf * (qty(i,j,k,qty_index) + qty(i,j,k-1,qty_index));
175  } else {
176 
177  amrex::Real avg1 = amrex::Real(0); amrex::Real avg2 = amrex::Real(0); amrex::Real avg3 = amrex::Real(0);
178  amrex::Real diff1 = amrex::Real(0); amrex::Real diff2 = amrex::Real(0); amrex::Real diff3 = amrex::Real(0);
179  amrex::Real scaled_upw = amrex::Real(0);
180  // The value that comes in has not been normalized so we do that here
181  if (upw != amrex::Real(0)) { scaled_upw = (upw > 0) ? amrex::Real(1) : -amrex::Real(1); }
182 
183  avg1 = (qty(i, j, k , qty_index) + qty(i, j, k-1, qty_index));
184  diff1 = (qty(i, j, k , qty_index) - qty(i, j, k-1, qty_index));
185  avg2 = (qty(i, j, k+1, qty_index) + qty(i, j, k-2, qty_index));
186  diff2 = (qty(i, j, k+1, qty_index) - qty(i, j, k-2, qty_index));
187  if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th)
188  {
189  avg3 = (qty(i, j, k+2, qty_index) + qty(i, j, k-3, qty_index));
190  diff3 = (qty(i, j, k+2, qty_index) - qty(i, j, k-3, qty_index));
191  }
192  return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type);
193  }
194 }
195 
196 /**
197  * Interpolate the perturbation of a quantity from cell center to face.
198  *
199  * @param i x-index
200  * @param j y-index
201  * @param k z-index
202  * @param qty Quantity field
203  * @param qty_index Component index
204  * @param upw Upwind direction
205  * @param coordDir Coordinate direction
206  * @param adv_type Advection scheme type
207  * @param r0_arr Base state field
208  * @return Interpolated perturbation value
209  */
210 AMREX_GPU_DEVICE
211 AMREX_FORCE_INLINE
212 amrex::Real
213 InterpolatePertFromCell (int i, int j, int k,
214  const amrex::Array4<const amrex::Real>& qty,
215  int qty_index, amrex::Real upw, Coord coordDir,
216  const AdvType adv_type, const amrex::Array4<const amrex::Real>& r0_arr)
217 {
218  // Reconstruct the perturbation qty - r0_arr, so both averages and
219  // differences must be formed from perturbation quantities.
220  amrex::Real avg1 = amrex::Real(0); amrex::Real avg2 = amrex::Real(0); amrex::Real avg3 = amrex::Real(0);
221  amrex::Real diff1 = amrex::Real(0); amrex::Real diff2 = amrex::Real(0); amrex::Real diff3 = amrex::Real(0);
222  amrex::Real scaled_upw = amrex::Real(0);
223 
224  // The value that comes in has not been normalized so we do that here
225  if (upw != amrex::Real(0)) { scaled_upw = (upw > 0) ? amrex::Real(1) : -amrex::Real(1); }
226 
227  if (coordDir == Coord::x) {
228  avg1 = (qty(i , j, k, qty_index) + qty(i-1, j, k, qty_index));
229  avg1 -= (r0_arr(i,j,k) + r0_arr(i-1,j,k));
230  diff1 = (qty(i , j, k, qty_index) - qty(i-1, j, k, qty_index));
231  diff1 -= (r0_arr(i,j,k) - r0_arr(i-1,j,k));
232  if (adv_type != AdvType::Centered_2nd)
233  {
234  avg2 = (qty(i+1, j, k, qty_index) + qty(i-2, j, k, qty_index));
235  avg2 -= (r0_arr(i+1,j,k) + r0_arr(i-2,j,k));
236  diff2 = (qty(i+1, j, k, qty_index) - qty(i-2, j, k, qty_index));
237  diff2 -= (r0_arr(i+1,j,k) - r0_arr(i-2,j,k));
238  }
239  if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th)
240  {
241  avg3 = (qty(i+2, j, k, qty_index) + qty(i-3, j, k, qty_index));
242  avg3 -= (r0_arr(i+2,j,k) + r0_arr(i-3,j,k));
243  diff3 = (qty(i+2, j, k, qty_index) - qty(i-3, j, k, qty_index));
244  diff3 -= (r0_arr(i+2,j,k) - r0_arr(i-3,j,k));
245  }
246  } else if (coordDir == Coord::y) {
247  avg1 = (qty(i, j , k, qty_index) + qty(i, j-1, k, qty_index));
248  avg1 -= (r0_arr(i,j,k) + r0_arr(i,j-1,k));
249  diff1 = (qty(i, j , k, qty_index) - qty(i, j-1, k, qty_index));
250  diff1 -= (r0_arr(i,j,k) - r0_arr(i,j-1,k));
251  if (adv_type != AdvType::Centered_2nd)
252  {
253  avg2 = (qty(i, j+1, k, qty_index) + qty(i, j-2, k, qty_index));
254  avg2 -= (r0_arr(i,j+1,k) + r0_arr(i,j-2,k));
255  diff2 = (qty(i, j+1, k, qty_index) - qty(i, j-2, k, qty_index));
256  diff2 -= (r0_arr(i,j+1,k) - r0_arr(i,j-2,k));
257  }
258  if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th)
259  {
260  avg3 = (qty(i, j+2, k, qty_index) + qty(i, j-3, k, qty_index));
261  avg3 -= (r0_arr(i,j+2,k) + r0_arr(i,j-3,k));
262  diff3 = (qty(i, j+2, k, qty_index) - qty(i, j-3, k, qty_index));
263  diff3 -= (r0_arr(i,j+2,k) - r0_arr(i,j-3,k));
264  }
265  } else {
266  avg1 = (qty(i, j, k , qty_index) + qty(i, j, k-1, qty_index));
267  diff1 = (qty(i, j, k , qty_index) - qty(i, j, k-1, qty_index));
268  avg1 -= (r0_arr(i,j,k) + r0_arr(i,j,k-1));
269  diff1 -= (r0_arr(i,j,k) - r0_arr(i,j,k-1));
270 
271  if (adv_type != AdvType::Centered_2nd)
272  {
273  avg2 = (qty(i, j, k+1, qty_index) + qty(i, j, k-2, qty_index));
274  diff2 = (qty(i, j, k+1, qty_index) - qty(i, j, k-2, qty_index));
275  avg2 -= (r0_arr(i,j,k+1) + r0_arr(i,j,k-2));
276  diff2 -= (r0_arr(i,j,k+1) - r0_arr(i,j,k-2));
277  }
278  if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th)
279  {
280  avg3 = (qty(i, j, k+2, qty_index) + qty(i, j, k-3, qty_index));
281  diff3 = (qty(i, j, k+2, qty_index) - qty(i, j, k-3, qty_index));
282  avg3 -= (r0_arr(i,j,k+2) + r0_arr(i,j,k-3));
283  diff3 -= (r0_arr(i,j,k+2) - r0_arr(i,j,k-3));
284  }
285  }
286 
287  return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type);
288 }
289 
290 /**
291  * Interpolate the density perturbation from cell center to face.
292  *
293  * @param i x-index
294  * @param j y-index
295  * @param k z-index
296  * @param cons_in Conservative variable field
297  * @param upw Upwind direction
298  * @param coordDir Coordinate direction
299  * @param adv_type Advection scheme type
300  * @param r0_arr Base state field
301  * @return Interpolated density perturbation
302  */
303 AMREX_GPU_DEVICE
304 AMREX_FORCE_INLINE
305 amrex::Real
306 InterpolateDensityPertFromCellToFace (int i, int j, int k, const amrex::Array4<const amrex::Real>& cons_in,
307  amrex::Real upw, Coord coordDir, const AdvType adv_type,
308  const amrex::Array4<const amrex::Real>& r0_arr)
309 {
310  return InterpolatePertFromCell(i, j, k, cons_in, Rho_comp, upw, coordDir, adv_type, r0_arr);
311 }
312 #endif
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
ERF_DataStruct.H
Coord
Coord
Coordinate-axis selector.
Definition: ERF_DataStruct.H:144
Coord::y
@ y
Coord::x
@ x
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
AdvType
AdvType
Definition: ERF_IndexDefines.H:304
AdvType::Upwind_3rd
@ Upwind_3rd
AdvType::Centered_4th
@ Centered_4th
AdvType::Upwind_5th
@ Upwind_5th
AdvType::Centered_6th
@ Centered_6th
AdvType::Centered_2nd
@ Centered_2nd
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
InterpolateInZ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateInZ(int i, int j, int k, const amrex::Array4< const amrex::Real > &qty, int qty_index, amrex::Real upw, const AdvType adv_type)
Definition: ERF_Interpolation.H:170
isSupportedInterpolationAdvType
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isSupportedInterpolationAdvType(const AdvType adv_type) noexcept
Definition: ERF_Interpolation.H:17
InterpolatePertFromCell
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolatePertFromCell(int i, int j, int k, const amrex::Array4< const amrex::Real > &qty, int qty_index, amrex::Real upw, Coord coordDir, const AdvType adv_type, const amrex::Array4< const amrex::Real > &r0_arr)
Definition: ERF_Interpolation.H:213
InterpolateDensityPertFromCellToFace
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateDensityPertFromCellToFace(int i, int j, int k, const amrex::Array4< const amrex::Real > &cons_in, amrex::Real upw, Coord coordDir, const AdvType adv_type, const amrex::Array4< const amrex::Real > &r0_arr)
Definition: ERF_Interpolation.H:306
interpolatedVal
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real interpolatedVal(amrex::Real avg1, amrex::Real avg2, amrex::Real avg3, amrex::Real diff1, amrex::Real diff2, amrex::Real diff3, amrex::Real scaled_upw, const AdvType adv_type)
Definition: ERF_Interpolation.H:42
InterpolateInY
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateInY(int i, int j, int k, const amrex::Array4< const amrex::Real > &qty, int qty_index, amrex::Real upw, const AdvType adv_type)
Definition: ERF_Interpolation.H:128
InterpolateInX
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateInX(int i, int j, int k, const amrex::Array4< const amrex::Real > &qty, int qty_index, amrex::Real upw, const AdvType adv_type)
Definition: ERF_Interpolation.H:85
ERF_Interpolation_UPW.H
ERF_Interpolation_WENO.H
ERF_Interpolation_WENO_Z.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19