ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_TerrainPoisson_3D_K.H
Go to the documentation of this file.
1 /**
2  * \file ERF_TerrainPoisson_3D_K.H
3  */
4 #ifndef ERF_TERRAINPOISSON_3D_K_H_
5 #define ERF_TERRAINPOISSON_3D_K_H_
6 
7 #include <AMReX_FArrayBox.H>
8 #include "ERF_Constants.H"
9 
10 /**
11  * Compute the x-face terrain Poisson flux.
12  *
13  * @tparam T Real-valued array and spacing type
14  * @param i x-index of the x-face
15  * @param j y-index of the x-face
16  * @param k z-index of the x-face
17  * @param sol Cell-centered solution field
18  * @param zp Node-centered physical height field
19  * @param dxinv Inverse x grid spacing
20  * @return x-face flux for the terrain Poisson operator
21  */
22 template <typename T>
23 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
24 T terrpoisson_flux_x (int i, int j, int k,
25  amrex::Array4<T const> const& sol,
26  amrex::Array4<T const> const& zp,
27  T dxinv) noexcept
28 {
29  using amrex::Real;
30 
31  Real h_xi, h_zeta;
32 
33  // On x-face
34  Real px_lo = (sol(i,j,k) - sol(i-1,j,k)) * dxinv;
35 
36  // On y-edges
37  Real pz_lo_md_hi = myhalf * ( sol(i,j,k+1) + sol(i-1,j,k+1)
38  -sol(i,j,k ) - sol(i-1,j,k ) );
39  h_xi = fourth * ( zp(i,j ,k+1) - zp(i-2,j ,k+1)
40  +zp(i,j+1,k+1) - zp(i-2,j+1,k+1) ) * dxinv;
41  h_zeta = fourth * ( zp(i,j ,k+2) - zp(i ,j ,k )
42  +zp(i,j+1,k+2) - zp(i ,j+1,k ) );
43  pz_lo_md_hi *= h_xi / h_zeta;
44 
45  // On y-edges
46  Real pz_lo_md_lo = myhalf * ( sol(i,j,k ) + sol(i-1,j,k )
47  -sol(i,j,k-1) - sol(i-1,j,k-1) );
48 
49  h_xi = fourth * ( zp(i,j ,k ) - zp(i-2,j ,k )
50  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
51  h_zeta = fourth * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
52  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
53  pz_lo_md_lo *= h_xi / h_zeta;
54 
55  // On x-face
56  px_lo -= myhalf * ( pz_lo_md_hi + pz_lo_md_lo );
57 
58  return -px_lo;
59 }
60 
61 /**
62  * Compute the y-face terrain Poisson flux.
63  *
64  * @tparam T Real-valued array and spacing type
65  * @param i x-index of the y-face
66  * @param j y-index of the y-face
67  * @param k z-index of the y-face
68  * @param sol Cell-centered solution field
69  * @param zp Node-centered physical height field
70  * @param dyinv Inverse y grid spacing
71  * @return y-face flux for the terrain Poisson operator
72  */
73 template <typename T>
74 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
75 T terrpoisson_flux_y (int i, int j, int k,
76  amrex::Array4<T const> const& sol,
77  amrex::Array4<T const> const& zp,
78  T dyinv) noexcept
79 {
80  using amrex::Real;
81 
82  Real h_eta, h_zeta;
83 
84  Real py_lo = (sol(i,j,k) - sol(i,j-1,k)) * dyinv;
85 
86  // On x-edges
87  Real pz_md_lo_hi = myhalf * ( sol(i,j,k+1) + sol(i,j-1,k+1)
88  -sol(i,j,k ) - sol(i,j-1,k ) );
89  h_eta = fourth * ( zp(i ,j,k+1) - zp(i ,j-2,k+1)
90  +zp(i+1,j,k+1) - zp(i+1,j-2,k+1) ) * dyinv;
91  h_zeta = fourth * ( zp(i ,j,k+2) - zp(i ,j ,k )
92  +zp(i+1,j,k+2) - zp(i+1,j ,k ) );
93  pz_md_lo_hi *= h_eta/ h_zeta;
94 
95  // On x-edges
96  Real pz_md_lo_lo = myhalf * ( sol(i,j,k ) + sol(i,j-1,k )
97  -sol(i,j,k-1) - sol(i,j-1,k-1) );
98 
99  h_eta = fourth * ( zp(i ,j,k ) - zp(i ,j-2,k )
100  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
101  h_zeta = fourth * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
102  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
103  pz_md_lo_lo *= h_eta/ h_zeta;
104 
105  // On y-face
106  py_lo -= myhalf * ( pz_md_lo_hi + pz_md_lo_lo );
107  return -py_lo;
108 }
109 
110 /**
111  * Compute the z-face terrain Poisson flux.
112  *
113  * @tparam T Real-valued array and spacing type
114  * @param i x-index of the z-face
115  * @param j y-index of the z-face
116  * @param k z-index of the z-face
117  * @param sol Cell-centered solution field
118  * @param zp Node-centered physical height field
119  * @param dxinv Inverse x grid spacing
120  * @param dyinv Inverse y grid spacing
121  * @return z-face flux for the terrain Poisson operator
122  */
123 template <typename T>
124 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
125 T terrpoisson_flux_z (int i, int j, int k,
126  amrex::Array4<T const> const& sol,
127  amrex::Array4<T const> const& zp,
128  T dxinv, T dyinv) noexcept
129 {
130  using amrex::Real;
131 
132  Real h_xi, h_eta, h_zeta;
133 
134  // On z-face
135  Real pz_lo = (sol(i,j,k ) - sol(i,j,k-1));
136  Real hzeta_inv_on_zlo = amrex::Real(8.0) / ( (zp(i,j,k+1) + zp(i+1,j,k+1) + zp(i,j+1,k+1) + zp(i+1,j+1,k+1))
137  -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
138  pz_lo *= hzeta_inv_on_zlo;
139 
140  // On corners
141  Real px_hi_md_lo = myhalf * ( sol(i+1,j ,k ) - sol(i ,j ,k )
142  +sol(i+1,j ,k-1) - sol(i ,j ,k-1)) * dxinv;
143  Real px_lo_md_lo = myhalf * ( sol(i ,j ,k ) - sol(i-1,j ,k )
144  +sol(i ,j ,k-1) - sol(i-1,j ,k-1)) * dxinv;
145  Real py_md_hi_lo = myhalf * ( sol(i ,j+1,k ) - sol(i ,j ,k )
146  +sol(i ,j+1,k-1) - sol(i ,j ,k-1)) * dyinv;
147  Real py_md_lo_lo = myhalf * ( sol(i ,j ,k ) - sol(i ,j-1,k )
148  +sol(i ,j ,k-1) - sol(i ,j-1,k-1)) * dyinv;
149 
150  // On y-edges
151  Real pz_hi_md_lo = myhalf * ( sol(i+1,j,k ) + sol(i,j,k )
152  -sol(i+1,j,k-1) - sol(i,j,k-1) );
153  h_xi = fourth * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
154  +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
155  h_zeta = fourth * ( zp(i+1,j ,k+1) - zp(i+1,j ,k-1)
156  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
157  pz_hi_md_lo *= h_xi / h_zeta;
158 
159  // On y-edges
160  Real pz_lo_md_lo = myhalf * ( sol(i,j,k ) + sol(i-1,j,k )
161  -sol(i,j,k-1) - sol(i-1,j,k-1) );
162  h_xi = fourth * ( zp(i,j ,k ) - zp(i-2,j ,k )
163  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
164  h_zeta = fourth * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
165  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
166  pz_lo_md_lo *= h_xi / h_zeta;
167 
168  // On x-edges
169  Real pz_md_hi_lo = myhalf * ( sol(i,j+1,k ) + sol(i,j,k )
170  -sol(i,j+1,k-1) - sol(i,j,k-1) );
171  h_eta = fourth * ( zp(i ,j+1,k ) - zp(i ,j-1,k)
172  +zp(i+1,j+1,k ) - zp(i+1,j-1,k) ) * dyinv;
173  h_zeta = fourth * ( zp(i ,j+1,k+1) - zp(i ,j+1,k-1)
174  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
175  pz_md_hi_lo *= h_eta/ h_zeta;
176 
177  // On x-edges
178  Real pz_md_lo_lo = myhalf * ( sol(i,j,k ) + sol(i,j-1,k )
179  -sol(i,j,k-1) - sol(i,j-1,k-1) );
180  h_eta = fourth * ( zp(i ,j,k ) - zp(i ,j-2,k )
181  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
182  h_zeta = fourth * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
183  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
184  pz_md_lo_lo *= h_eta/ h_zeta;
185 
186  // On z-face
187  Real h_xi_on_zlo = myhalf * (zp(i+1,j+1,k) + zp(i+1,j,k) - zp(i,j+1,k) - zp(i,j,k)) * dxinv;
188  Real h_eta_on_zlo = myhalf * (zp(i+1,j+1,k) + zp(i,j+1,k) - zp(i+1,j,k) - zp(i,j,k)) * dyinv;
189 
190  pz_lo -= myhalf * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
191  pz_lo -= myhalf * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
192 
193  if (k == 0) pz_lo = zero;
194 
195  return -pz_lo;
196 }
197 
198 /**
199  * Simplified gradient calc with dirichlet BC (phi_zlo==0)
200  * - compare with terrpoisson_flux_z
201  *
202  * @tparam T Real-valued array and spacing type
203  * @param i x-index of the z-face
204  * @param j y-index of the z-face
205  * @param k z-index of the z-face
206  * @param sol Cell-centered solution field
207  * @param zp Node-centered physical height field
208  * @param dxinv Inverse x grid spacing
209  * @param dyinv Inverse y grid spacing
210  * @return z-face flux with a low-z Dirichlet value
211  */
212 template <typename T>
213 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
214 T terrpoisson_flux_zlo_dir (int i, int j, int k,
215  amrex::Array4<T const> const& sol,
216  amrex::Array4<T const> const& zp,
217  T dxinv, T dyinv) noexcept
218 {
219  using amrex::Real;
220 
221  Real h_xi, h_eta, h_zeta;
222 
223  // On z-face, one-sided
224  Real pz_lo = sol(i,j,k);
225  Real hzeta_inv_on_zlo = amrex::Real(8.0) / ( (zp(i,j,k+1) + zp(i+1,j,k+1) + zp(i,j+1,k+1) + zp(i+1,j+1,k+1))
226  -(zp(i,j,k ) + zp(i+1,j,k ) + zp(i,j+1,k ) + zp(i+1,j+1,k )) );
227  pz_lo *= hzeta_inv_on_zlo;
228 
229  // On corners
230  constexpr Real px_hi_md_lo = zero;
231  constexpr Real px_lo_md_lo = zero;
232  constexpr Real py_md_hi_lo = zero;
233  constexpr Real py_md_lo_lo = zero;
234 
235  // On y-edges, one-sided
236  Real pz_hi_md_lo = myhalf * ( sol(i+1,j,k ) + sol(i,j,k ));
237  h_xi = fourth * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
238  +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
239  h_zeta = fourth * ( zp(i+1,j ,k+1) - zp(i+1,j ,k )
240  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k ) );
241  pz_hi_md_lo *= h_xi / h_zeta;
242 
243  // On y-edges, one-sided
244  Real pz_lo_md_lo = myhalf * ( sol(i,j,k ) + sol(i-1,j,k ));
245  h_xi = fourth * ( zp(i,j ,k ) - zp(i-2,j ,k )
246  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
247  h_zeta = fourth * ( zp(i,j ,k+1) - zp(i ,j ,k )
248  +zp(i,j+1,k+1) - zp(i ,j+1,k ) );
249  pz_lo_md_lo *= h_xi / h_zeta;
250 
251  // On x-edges, one-sided
252  Real pz_md_hi_lo = myhalf * ( sol(i,j+1,k ) + sol(i,j,k ));
253  h_eta = fourth * ( zp(i ,j+1,k ) - zp(i ,j-1,k)
254  +zp(i+1,j+1,k ) - zp(i+1,j-1,k) ) * dyinv;
255  h_zeta = fourth * ( zp(i ,j+1,k+1) - zp(i ,j+1,k )
256  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k ) );
257  pz_md_hi_lo *= h_eta/ h_zeta;
258 
259  // On x-edges, one-sided
260  Real pz_md_lo_lo = myhalf * ( sol(i,j,k ) + sol(i,j-1,k ));
261  h_eta = fourth * ( zp(i ,j,k ) - zp(i ,j-2,k )
262  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
263  h_zeta = fourth * ( zp(i ,j,k+1) - zp(i ,j ,k )
264  +zp(i+1,j,k+1) - zp(i+1,j ,k ) );
265  pz_md_lo_lo *= h_eta/ h_zeta;
266 
267  // On z-face
268  Real h_xi_on_zlo = myhalf * (zp(i+1,j+1,k) + zp(i+1,j,k) - zp(i,j+1,k) - zp(i,j,k)) * dxinv;
269  Real h_eta_on_zlo = myhalf * (zp(i+1,j+1,k) + zp(i,j+1,k) - zp(i+1,j,k) - zp(i,j,k)) * dyinv;
270 
271  pz_lo -= myhalf * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
272  pz_lo -= myhalf * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
273 
274  return -pz_lo;
275 }
276 
277 /**
278  * Apply the terrain Poisson matrix to a cell-centered field.
279  *
280  * @tparam T Real-valued array and spacing type
281  * @param i x-index of the cell
282  * @param j y-index of the cell
283  * @param k z-index of the cell
284  * @param y Destination operator result
285  * @param x Source solution field
286  * @param ax Terrain metric coefficient on x-faces
287  * @param ay Terrain metric coefficient on y-faces
288  * @param az Terrain metric coefficient on z-faces
289  * @param dJ Cell-centered Jacobian determinant
290  * @param zp Node-centered physical height field
291  * @param dxinv Inverse x grid spacing
292  * @param dyinv Inverse y grid spacing
293  * @param dzinv Inverse z grid spacing
294  */
295 template <typename T>
296 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
297 void terrpoisson_adotx (int i, int j, int k, amrex::Array4<T> const& y,
298  amrex::Array4<T const> const& x,
299  amrex::Array4<T const> const& ax,
300  amrex::Array4<T const> const& ay,
301  amrex::Array4<T const> const& az,
302  amrex::Array4<T const> const& dJ,
303  amrex::Array4<T const> const& zp,
304  T dxinv, T dyinv, T dzinv) noexcept
305 {
306  using amrex::Real;
307  Real h_xi, h_eta, h_zeta;
308 
309 #if 1
310  // *********************************************************
311  // Hi x-face
312  // *********************************************************
313  // On x-face
314  Real px_hi = (x(i+1,j,k) - x(i,j,k)) * dxinv;
315 
316  // On y-edges
317  Real pz_hi_md_hi = myhalf * ( x(i+1,j ,k+1) + x(i ,j ,k+1)
318  -x(i+1,j ,k ) - x(i ,j ,k ) );
319  h_xi = fourth * ( zp(i+1,j ,k+1) - zp(i-1,j ,k+1)
320  +zp(i+1,j+1,k+1) - zp(i-1,j+1,k+1) ) * dxinv;
321  h_zeta = fourth * ( zp(i+1,j ,k+2) - zp(i+1,j ,k )
322  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
323  pz_hi_md_hi *= h_xi / h_zeta;
324 
325  // On y-edges
326  Real pz_hi_md_lo = myhalf * ( x(i+1,j ,k ) + x(i ,j ,k )
327  -x(i+1,j ,k-1) - x(i ,j ,k-1) );
328  h_xi = fourth * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
329  +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
330  h_zeta = fourth * ( zp(i+1,j ,k+1) - zp(i+1,j ,k-1)
331  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
332  pz_hi_md_lo *= h_xi / h_zeta;
333 
334  // On x-face
335  px_hi -= myhalf * ( pz_hi_md_hi + pz_hi_md_lo );
336 
337  // *********************************************************
338  // Lo x-face
339  // ********************************************************* // On x-face
340  Real px_lo = (x(i,j,k) - x(i-1,j,k)) * dxinv;
341 
342  // On y-edges
343  Real pz_lo_md_hi = myhalf * ( x(i,j,k+1) + x(i-1,j,k+1)
344  -x(i,j,k ) - x(i-1,j,k ) );
345  h_xi = fourth * ( zp(i,j ,k+1) - zp(i-2,j ,k+1)
346  +zp(i,j+1,k+1) - zp(i-2,j+1,k+1) ) * dxinv;
347  h_zeta = fourth * ( zp(i,j ,k+2) - zp(i ,j ,k )
348  +zp(i,j+1,k+2) - zp(i ,j+1,k ) );
349  pz_lo_md_hi *= h_xi / h_zeta;
350 
351  // On y-edges
352  Real pz_lo_md_lo = myhalf * ( x(i,j,k ) + x(i-1,j,k )
353  -x(i,j,k-1) - x(i-1,j,k-1) );
354  h_xi = fourth * ( zp(i,j ,k ) - zp(i-2,j ,k )
355  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
356  h_zeta = fourth * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
357  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
358  pz_lo_md_lo *= h_xi / h_zeta;
359 
360  // On x-face
361  px_lo -= myhalf * ( pz_lo_md_hi + pz_lo_md_lo );
362 
363  // *********************************************************
364  // Hi y-face
365  // *********************************************************
366  // On y-face
367  Real py_hi = (x(i,j+1,k) - x(i,j,k)) * dyinv;
368 
369  // On x-edges
370  Real pz_md_hi_hi = myhalf * ( x(i,j+1,k+1) + x(i,j,k+1)
371  -x(i,j+1,k ) - x(i,j,k ) );
372  h_eta = fourth * ( zp(i ,j+1,k+1) - zp(i ,j-1,k+1)
373  +zp(i+1,j+1,k+1) - zp(i+1,j-1,k+1) ) * dyinv;
374  h_zeta = fourth * ( zp(i ,j+1,k+2) - zp(i ,j+1,k )
375  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
376  pz_md_hi_hi *= h_eta/ h_zeta;
377 
378  // On x-edges
379  Real pz_md_hi_lo = myhalf * ( x(i,j+1,k ) + x(i,j,k )
380  -x(i,j+1,k-1) - x(i,j,k-1) );
381  h_eta = fourth * ( zp(i ,j+1,k ) - zp(i ,j-1,k)
382  +zp(i+1,j+1,k ) - zp(i+1,j-1,k) ) * dyinv;
383  h_zeta = fourth * ( zp(i ,j+1,k+1) - zp(i ,j+1,k-1)
384  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
385  pz_md_hi_lo *= h_eta/ h_zeta;
386 
387  // On y-face
388  py_hi -= myhalf * ( pz_md_hi_hi + pz_md_hi_lo );
389 
390  // *********************************************************
391  // Lo y-face
392  // *********************************************************
393  // On y-face
394  Real py_lo = (x(i,j,k) - x(i,j-1,k)) * dyinv;
395 
396  // On x-edges
397  Real pz_md_lo_hi = myhalf * ( x(i ,j,k+1) + x(i ,j-1,k+1)
398  -x(i ,j,k ) - x(i ,j-1,k ) );
399  h_eta = fourth * ( zp(i ,j,k+1) - zp(i ,j-2,k+1)
400  +zp(i+1,j,k+1) - zp(i+1,j-2,k+1) ) * dyinv;
401  h_zeta = fourth * ( zp(i ,j,k+2) - zp(i ,j ,k )
402  +zp(i+1,j,k+2) - zp(i+1,j ,k ) );
403  pz_md_lo_hi *= h_eta/ h_zeta;
404 
405  // On x-edges
406  Real pz_md_lo_lo = myhalf * ( x(i ,j,k ) + x(i ,j-1,k )
407  -x(i ,j,k-1) - x(i ,j-1,k-1) );
408  h_eta = fourth * ( zp(i ,j,k ) - zp(i ,j-2,k )
409  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
410  h_zeta = fourth * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
411  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
412  pz_md_lo_lo *= h_eta/ h_zeta;
413 
414  // On y-face
415  py_lo -= myhalf * ( pz_md_lo_hi + pz_md_lo_lo );
416 
417  // *********************************************************
418  // Hi z-face
419  // *********************************************************
420  // On z-face
421  Real pz_hi = x(i,j,k+1) - x(i,j,k );
422  Real hzeta_inv_on_zhi = amrex::Real(8.0) / ( (zp(i,j,k+2) + zp(i+1,j,k+2) + zp(i,j+1,k+2) + zp(i+1,j+1,k+2))
423  -(zp(i,j,k ) + zp(i+1,j,k ) + zp(i,j+1,k ) + zp(i+1,j+1,k )) );
424  pz_hi *= hzeta_inv_on_zhi;
425 
426  // On corners
427  Real px_hi_md_hi = myhalf * ( x(i+1,j,k+1) - x(i ,j ,k+1)
428  +x(i+1,j,k ) - x(i ,j ,k )) * dxinv;
429  Real px_lo_md_hi = myhalf * ( x(i ,j,k+1) - x(i-1,j ,k+1)
430  +x(i ,j,k ) - x(i-1,j ,k )) * dxinv;
431  Real py_md_hi_hi = myhalf * ( x(i,j+1,k+1) - x(i ,j ,k+1)
432  +x(i,j+1,k ) - x(i ,j ,k )) * dyinv;
433  Real py_md_lo_hi = myhalf * ( x(i,j ,k+1) - x(i ,j-1,k+1)
434  +x(i,j ,k ) - x(i ,j-1,k )) * dyinv;
435 
436  // On z-face
437  Real h_xi_on_zhi = myhalf * ( zp(i+1,j+1,k+1) + zp(i+1,j,k+1) - zp(i,j+1,k+1) - zp(i,j,k+1) ) * dxinv;
438  Real h_eta_on_zhi = myhalf * ( zp(i+1,j+1,k+1) + zp(i,j+1,k+1) - zp(i+1,j,k+1) - zp(i,j,k+1) ) * dyinv;
439  //
440  // Note we do not need to recalculate pz_...hi here
441  //
442  pz_hi -= myhalf * h_xi_on_zhi * ( (px_hi_md_hi + px_lo_md_hi) - (pz_hi_md_hi + pz_lo_md_hi) );
443  pz_hi -= myhalf * h_eta_on_zhi * ( (py_md_hi_hi + py_md_lo_hi) - (pz_md_hi_hi + pz_md_lo_hi) );
444 
445  // *********************************************************
446  // Lo z-face
447  // *********************************************************
448  // On z-face
449  Real pz_lo = x(i,j,k ) - x(i,j,k-1);
450  Real hzeta_inv_on_zlo = amrex::Real(8.0) / ( (zp(i,j,k+1) + zp(i+1,j,k+1) + zp(i,j+1,k+1) + zp(i+1,j+1,k+1))
451  -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
452  pz_lo *= hzeta_inv_on_zlo;
453 
454  // On corners
455  Real px_hi_md_lo = myhalf * ( x(i+1,j ,k ) - x(i ,j ,k )
456  +x(i+1,j ,k-1) - x(i ,j ,k-1)) * dxinv;
457  Real px_lo_md_lo = myhalf * ( x(i ,j ,k ) - x(i-1,j ,k )
458  +x(i ,j ,k-1) - x(i-1,j ,k-1)) * dxinv;
459  Real py_md_hi_lo = myhalf * ( x(i ,j+1,k ) - x(i ,j ,k )
460  +x(i ,j+1,k-1) - x(i ,j ,k-1)) * dyinv;
461  Real py_md_lo_lo = myhalf * ( x(i ,j ,k ) - x(i ,j-1,k )
462  +x(i ,j ,k-1) - x(i ,j-1,k-1)) * dyinv;
463 
464  // On z-face
465  Real h_xi_on_zlo = myhalf * (zp(i+1,j+1,k) + zp(i+1,j,k) - zp(i,j+1,k) - zp(i,j,k)) * dxinv;
466  Real h_eta_on_zlo = myhalf * (zp(i+1,j+1,k) + zp(i,j+1,k) - zp(i+1,j,k) - zp(i,j,k)) * dyinv;
467  //
468  // Note we do not need to recalculate pz_...lo here
469  //
470  pz_lo -= myhalf * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
471  pz_lo -= myhalf * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
472 
473 #else
474 
475  // *********************************************************
476  // Version which calls flux routines
477  // *********************************************************
478  //
479  // This option uses calls to the flux routines so there is
480  // some duplicated computation
481  // This option should give the same answer as above
482  //
483  Real px_lo = -terrpoisson_flux_x(i ,j,k,x,zp,dxinv);
484  Real px_hi = -terrpoisson_flux_x(i+1,j,k,x,zp,dxinv);
485  Real py_lo = -terrpoisson_flux_y(i,j ,k,x,zp,dyinv);
486  Real py_hi = -terrpoisson_flux_y(i,j+1,k,x,zp,dyinv);
487  Real pz_lo = -terrpoisson_flux_z(i,j,k ,x,zp,dxinv,dyinv);
488  Real pz_hi = -terrpoisson_flux_z(i,j,k+1,x,zp,dxinv,dyinv);
489 #endif
490  //
491  // *********************************************************
492  // Adotx
493  // *********************************************************
494  if (k == 0) pz_lo = zero;
495 
496  y(i,j,k) = (ax(i+1,j,k)*px_hi - ax(i,j,k)*px_lo) * dxinv
497  +(ay(i,j+1,k)*py_hi - ay(i,j,k)*py_lo) * dyinv
498  +(az(i,j,k+1)*pz_hi - az(i,j,k)*pz_lo) * dzinv;
499  y(i,j,k) /= dJ(i,j,k);
500 }
501 #endif
ERF_Constants.H
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
Coord::y
@ y
Coord::x
@ x
T
Real T
Definition: ERF_InitCustomPert_Bubble.H:106
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
terrpoisson_flux_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_x(int i, int j, int k, amrex::Array4< T const > const &sol, amrex::Array4< T const > const &zp, T dxinv) noexcept
Definition: ERF_TerrainPoisson_3D_K.H:24
terrpoisson_flux_zlo_dir
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_zlo_dir(int i, int j, int k, amrex::Array4< T const > const &sol, amrex::Array4< T const > const &zp, T dxinv, T dyinv) noexcept
Definition: ERF_TerrainPoisson_3D_K.H:214
terrpoisson_flux_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_z(int i, int j, int k, amrex::Array4< T const > const &sol, amrex::Array4< T const > const &zp, T dxinv, T dyinv) noexcept
Definition: ERF_TerrainPoisson_3D_K.H:125
terrpoisson_adotx
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void terrpoisson_adotx(int i, int j, int k, amrex::Array4< T > const &y, amrex::Array4< T const > const &x, amrex::Array4< T const > const &ax, amrex::Array4< T const > const &ay, amrex::Array4< T const > const &az, amrex::Array4< T const > const &dJ, amrex::Array4< T const > const &zp, T dxinv, T dyinv, T dzinv) noexcept
Definition: ERF_TerrainPoisson_3D_K.H:297
terrpoisson_flux_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_y(int i, int j, int k, amrex::Array4< T const > const &sol, amrex::Array4< T const > const &zp, T dyinv) noexcept
Definition: ERF_TerrainPoisson_3D_K.H:75