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