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