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_NumericalConstants.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+1,j ,k+1) - zp(i-1,j ,k+1)
40  +zp(i+1,j+1,k+1) - zp(i-1,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+1,j ,k ) - zp(i-1,j ,k )
50  +zp(i+1,j+1,k ) - zp(i-1,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+1,k+1) - zp(i ,j-1,k+1)
90  +zp(i+1,j+1,k+1) - zp(i+1,j-1,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+1,k ) - zp(i ,j-1,k )
100  +zp(i+1,j+1,k ) - zp(i+1,j-1,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+2,j ,k ) - zp(i ,j ,k )
154  +zp(i+2,j+1,k ) - zp(i ,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+1,j ,k ) - zp(i-1,j ,k )
163  +zp(i+1,j+1,k ) - zp(i-1,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+2,k ) - zp(i ,j ,k )
172  +zp(i+1,j+2,k ) - zp(i+1,j ,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+1,k ) - zp(i ,j-1,k )
181  +zp(i+1,j+1,k ) - zp(i+1,j-1,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+2,j ,k ) - zp(i ,j ,k )
238  +zp(i+2,j+1,k ) - zp(i ,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+1,j ,k ) - zp(i-1,j ,k )
246  +zp(i+1,j+1,k ) - zp(i-1,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+2,k ) - zp(i ,j ,k )
254  +zp(i+1,j+2,k ) - zp(i+1,j ,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+1,k ) - zp(i ,j-1,k )
262  +zp(i+1,j+1,k ) - zp(i+1,j-1,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  * Bits of the face mask taken by terrpoisson_adotx: one per face of the solve region whose
279  * flux is set to zero (TerrainPoisson masks the Neumann faces of its region).
280  */
281 enum TerrainPoissonFaceBit : int {
282  tp_face_xlo = 1, tp_face_xhi = 2,
283  tp_face_ylo = 4, tp_face_yhi = 8,
284  tp_face_zlo = 16, tp_face_zhi = 32
285 };
286 
287 /**
288  * Apply the terrain Poisson matrix to a cell-centered field.
289  *
290  * @tparam T Real-valued array and spacing type
291  * @param i x-index of the cell
292  * @param j y-index of the cell
293  * @param k z-index of the cell
294  * @param y Destination operator result
295  * @param x Source solution field
296  * @param ax Terrain metric coefficient on x-faces
297  * @param ay Terrain metric coefficient on y-faces
298  * @param az Terrain metric coefficient on z-faces
299  * @param dJ Cell-centered Jacobian determinant
300  * @param zp Node-centered physical height field
301  * @param dxinv Inverse x grid spacing
302  * @param dyinv Inverse y grid spacing
303  * @param dzinv Inverse z grid spacing
304  * @param dlo Low corner of the solve region (used with zf)
305  * @param dhi High corner of the solve region (used with zf)
306  * @param zf Faces of the region whose flux is zero (sum of TerrainPoissonFaceBit values;
307  * 0, the default, keeps every face flux)
308  *
309  * A masked face is a Neumann face of the solve region: a face of the region inside the domain
310  * (a coarse/fine face of a refined level) or a Neumann (wall, inflow) face of the domain. The
311  * even mirror that TerrainPoisson::apply_bcs puts in the ghost cell zeroes the normal-gradient
312  * part of the flux through such a face but not the terrain cross term, so on sloping terrain
313  * the face would let correction through; with the face masked the whole flux through it is
314  * zero, as in AMReX's MLTerrainPoisson. Faces outside the mask are untouched, bit for bit.
315  */
316 template <typename T>
317 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
318 void terrpoisson_adotx (int i, int j, int k, amrex::Array4<T> const& y,
319  amrex::Array4<T const> const& x,
320  amrex::Array4<T const> const& ax,
321  amrex::Array4<T const> const& ay,
322  amrex::Array4<T const> const& az,
323  amrex::Array4<T const> const& dJ,
324  amrex::Array4<T const> const& zp,
325  T dxinv, T dyinv, T dzinv,
326  amrex::Dim3 const& dlo = amrex::Dim3{0,0,0},
327  amrex::Dim3 const& dhi = amrex::Dim3{0,0,0},
328  int zf = 0) noexcept
329 {
330  using amrex::Real;
331  Real h_xi, h_eta, h_zeta;
332 
333 #if 1
334  // *********************************************************
335  // Hi x-face
336  // *********************************************************
337  // On x-face
338  Real px_hi = (x(i+1,j,k) - x(i,j,k)) * dxinv;
339 
340  // On y-edges
341  Real pz_hi_md_hi = myhalf * ( x(i+1,j ,k+1) + x(i ,j ,k+1)
342  -x(i+1,j ,k ) - x(i ,j ,k ) );
343  h_xi = fourth * ( zp(i+2,j ,k+1) - zp(i ,j ,k+1)
344  +zp(i+2,j+1,k+1) - zp(i ,j+1,k+1) ) * dxinv;
345  h_zeta = fourth * ( zp(i+1,j ,k+2) - zp(i+1,j ,k )
346  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
347  pz_hi_md_hi *= h_xi / h_zeta;
348 
349  // On y-edges
350  Real pz_hi_md_lo = myhalf * ( x(i+1,j ,k ) + x(i ,j ,k )
351  -x(i+1,j ,k-1) - x(i ,j ,k-1) );
352  h_xi = fourth * ( zp(i+2,j ,k ) - zp(i ,j ,k )
353  +zp(i+2,j+1,k ) - zp(i ,j+1,k ) ) * dxinv;
354  h_zeta = fourth * ( zp(i+1,j ,k+1) - zp(i+1,j ,k-1)
355  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
356  pz_hi_md_lo *= h_xi / h_zeta;
357 
358  // On x-face
359  px_hi -= myhalf * ( pz_hi_md_hi + pz_hi_md_lo );
360 
361  // *********************************************************
362  // Lo x-face
363  // ********************************************************* // On x-face
364  Real px_lo = (x(i,j,k) - x(i-1,j,k)) * dxinv;
365 
366  // On y-edges
367  Real pz_lo_md_hi = myhalf * ( x(i,j,k+1) + x(i-1,j,k+1)
368  -x(i,j,k ) - x(i-1,j,k ) );
369  h_xi = fourth * ( zp(i+1,j ,k+1) - zp(i-1,j ,k+1)
370  +zp(i+1,j+1,k+1) - zp(i-1,j+1,k+1) ) * dxinv;
371  h_zeta = fourth * ( zp(i,j ,k+2) - zp(i ,j ,k )
372  +zp(i,j+1,k+2) - zp(i ,j+1,k ) );
373  pz_lo_md_hi *= h_xi / h_zeta;
374 
375  // On y-edges
376  Real pz_lo_md_lo = myhalf * ( x(i,j,k ) + x(i-1,j,k )
377  -x(i,j,k-1) - x(i-1,j,k-1) );
378  h_xi = fourth * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
379  +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
380  h_zeta = fourth * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
381  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
382  pz_lo_md_lo *= h_xi / h_zeta;
383 
384  // On x-face
385  px_lo -= myhalf * ( pz_lo_md_hi + pz_lo_md_lo );
386 
387  // *********************************************************
388  // Hi y-face
389  // *********************************************************
390  // On y-face
391  Real py_hi = (x(i,j+1,k) - x(i,j,k)) * dyinv;
392 
393  // On x-edges
394  Real pz_md_hi_hi = myhalf * ( x(i,j+1,k+1) + x(i,j,k+1)
395  -x(i,j+1,k ) - x(i,j,k ) );
396  h_eta = fourth * ( zp(i ,j+2,k+1) - zp(i ,j ,k+1)
397  +zp(i+1,j+2,k+1) - zp(i+1,j ,k+1) ) * dyinv;
398  h_zeta = fourth * ( zp(i ,j+1,k+2) - zp(i ,j+1,k )
399  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
400  pz_md_hi_hi *= h_eta/ h_zeta;
401 
402  // On x-edges
403  Real pz_md_hi_lo = myhalf * ( x(i,j+1,k ) + x(i,j,k )
404  -x(i,j+1,k-1) - x(i,j,k-1) );
405  h_eta = fourth * ( zp(i ,j+2,k ) - zp(i ,j ,k )
406  +zp(i+1,j+2,k ) - zp(i+1,j ,k ) ) * dyinv;
407  h_zeta = fourth * ( zp(i ,j+1,k+1) - zp(i ,j+1,k-1)
408  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
409  pz_md_hi_lo *= h_eta/ h_zeta;
410 
411  // On y-face
412  py_hi -= myhalf * ( pz_md_hi_hi + pz_md_hi_lo );
413 
414  // *********************************************************
415  // Lo y-face
416  // *********************************************************
417  // On y-face
418  Real py_lo = (x(i,j,k) - x(i,j-1,k)) * dyinv;
419 
420  // On x-edges
421  Real pz_md_lo_hi = myhalf * ( x(i ,j,k+1) + x(i ,j-1,k+1)
422  -x(i ,j,k ) - x(i ,j-1,k ) );
423  h_eta = fourth * ( zp(i ,j+1,k+1) - zp(i ,j-1,k+1)
424  +zp(i+1,j+1,k+1) - zp(i+1,j-1,k+1) ) * dyinv;
425  h_zeta = fourth * ( zp(i ,j,k+2) - zp(i ,j ,k )
426  +zp(i+1,j,k+2) - zp(i+1,j ,k ) );
427  pz_md_lo_hi *= h_eta/ h_zeta;
428 
429  // On x-edges
430  Real pz_md_lo_lo = myhalf * ( x(i ,j,k ) + x(i ,j-1,k )
431  -x(i ,j,k-1) - x(i ,j-1,k-1) );
432  h_eta = fourth * ( zp(i ,j+1,k ) - zp(i ,j-1,k )
433  +zp(i+1,j+1,k ) - zp(i+1,j-1,k ) ) * dyinv;
434  h_zeta = fourth * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
435  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
436  pz_md_lo_lo *= h_eta/ h_zeta;
437 
438  // On y-face
439  py_lo -= myhalf * ( pz_md_lo_hi + pz_md_lo_lo );
440 
441  // *********************************************************
442  // Hi z-face
443  // *********************************************************
444  // On z-face
445  Real pz_hi = x(i,j,k+1) - x(i,j,k );
446  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))
447  -(zp(i,j,k ) + zp(i+1,j,k ) + zp(i,j+1,k ) + zp(i+1,j+1,k )) );
448  pz_hi *= hzeta_inv_on_zhi;
449 
450  // On corners
451  Real px_hi_md_hi = myhalf * ( x(i+1,j,k+1) - x(i ,j ,k+1)
452  +x(i+1,j,k ) - x(i ,j ,k )) * dxinv;
453  Real px_lo_md_hi = myhalf * ( x(i ,j,k+1) - x(i-1,j ,k+1)
454  +x(i ,j,k ) - x(i-1,j ,k )) * dxinv;
455  Real py_md_hi_hi = myhalf * ( x(i,j+1,k+1) - x(i ,j ,k+1)
456  +x(i,j+1,k ) - x(i ,j ,k )) * dyinv;
457  Real py_md_lo_hi = myhalf * ( x(i,j ,k+1) - x(i ,j-1,k+1)
458  +x(i,j ,k ) - x(i ,j-1,k )) * dyinv;
459 
460  // On z-face
461  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;
462  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;
463  //
464  // Note we do not need to recalculate pz_...hi here
465  //
466  pz_hi -= myhalf * h_xi_on_zhi * ( (px_hi_md_hi + px_lo_md_hi) - (pz_hi_md_hi + pz_lo_md_hi) );
467  pz_hi -= myhalf * h_eta_on_zhi * ( (py_md_hi_hi + py_md_lo_hi) - (pz_md_hi_hi + pz_md_lo_hi) );
468 
469  // *********************************************************
470  // Lo z-face
471  // *********************************************************
472  // On z-face
473  Real pz_lo = x(i,j,k ) - x(i,j,k-1);
474  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))
475  -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
476  pz_lo *= hzeta_inv_on_zlo;
477 
478  // On corners
479  Real px_hi_md_lo = myhalf * ( x(i+1,j ,k ) - x(i ,j ,k )
480  +x(i+1,j ,k-1) - x(i ,j ,k-1)) * dxinv;
481  Real px_lo_md_lo = myhalf * ( x(i ,j ,k ) - x(i-1,j ,k )
482  +x(i ,j ,k-1) - x(i-1,j ,k-1)) * dxinv;
483  Real py_md_hi_lo = myhalf * ( x(i ,j+1,k ) - x(i ,j ,k )
484  +x(i ,j+1,k-1) - x(i ,j ,k-1)) * dyinv;
485  Real py_md_lo_lo = myhalf * ( x(i ,j ,k ) - x(i ,j-1,k )
486  +x(i ,j ,k-1) - x(i ,j-1,k-1)) * dyinv;
487 
488  // On z-face
489  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;
490  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;
491  //
492  // Note we do not need to recalculate pz_...lo here
493  //
494  pz_lo -= myhalf * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
495  pz_lo -= myhalf * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
496 
497 #else
498 
499  // *********************************************************
500  // Version which calls flux routines
501  // *********************************************************
502  //
503  // This option uses calls to the flux routines so there is
504  // some duplicated computation
505  // This option should give the same answer as above
506  //
507  Real px_lo = -terrpoisson_flux_x(i ,j,k,x,zp,dxinv);
508  Real px_hi = -terrpoisson_flux_x(i+1,j,k,x,zp,dxinv);
509  Real py_lo = -terrpoisson_flux_y(i,j ,k,x,zp,dyinv);
510  Real py_hi = -terrpoisson_flux_y(i,j+1,k,x,zp,dyinv);
511  Real pz_lo = -terrpoisson_flux_z(i,j,k ,x,zp,dxinv,dyinv);
512  Real pz_hi = -terrpoisson_flux_z(i,j,k+1,x,zp,dxinv,dyinv);
513 #endif
514  //
515  // *********************************************************
516  // Adotx
517  // *********************************************************
518  if (k == 0) pz_lo = zero;
519  if ((zf & tp_face_xlo) && i == dlo.x) px_lo = zero;
520  if ((zf & tp_face_xhi) && i == dhi.x) px_hi = zero;
521  if ((zf & tp_face_ylo) && j == dlo.y) py_lo = zero;
522  if ((zf & tp_face_yhi) && j == dhi.y) py_hi = zero;
523  if ((zf & tp_face_zlo) && k == dlo.z) pz_lo = zero;
524  if ((zf & tp_face_zhi) && k == dhi.z) pz_hi = zero;
525 
526  y(i,j,k) = (ax(i+1,j,k)*px_hi - ax(i,j,k)*px_lo) * dxinv
527  +(ay(i,j+1,k)*py_hi - ay(i,j,k)*py_lo) * dyinv
528  +(az(i,j,k+1)*pz_hi - az(i,j,k)*pz_lo) * dzinv;
529  y(i,j,k) /= dJ(i,j,k);
530 }
531 
532 
533 #endif
Coord::y
@ y
Coord::x
@ x
ERF_NumericalConstants.H
Dimensionless numeric literals and pure mathematical constants.
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TerrainPoissonFaceBit
TerrainPoissonFaceBit
Definition: ERF_TerrainPoisson_3D_K.H:281
tp_face_xlo
@ tp_face_xlo
Definition: ERF_TerrainPoisson_3D_K.H:282
tp_face_zhi
@ tp_face_zhi
Definition: ERF_TerrainPoisson_3D_K.H:284
tp_face_yhi
@ tp_face_yhi
Definition: ERF_TerrainPoisson_3D_K.H:283
tp_face_zlo
@ tp_face_zlo
Definition: ERF_TerrainPoisson_3D_K.H:284
tp_face_ylo
@ tp_face_ylo
Definition: ERF_TerrainPoisson_3D_K.H:283
tp_face_xhi
@ tp_face_xhi
Definition: ERF_TerrainPoisson_3D_K.H:282
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_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
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, amrex::Dim3 const &dlo=amrex::Dim3{0, 0, 0}, amrex::Dim3 const &dhi=amrex::Dim3{0, 0, 0}, int zf=0) noexcept
Definition: ERF_TerrainPoisson_3D_K.H:318
RealBdyVars::T
@ T
Definition: ERF_IndexDefines.H:131