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 template <typename T>
157 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
158 void terrpoisson_adotx (int i, int j, int k, amrex::Array4<T> const& y,
159  amrex::Array4<T const> const& x,
160  amrex::Array4<T const> const& ax,
161  amrex::Array4<T const> const& ay,
162  amrex::Array4<T const> const& az,
163  amrex::Array4<T const> const& dJ,
164  amrex::Array4<T const> const& zp,
165  T dxinv, T dyinv, T dzinv) noexcept
166 {
167  using amrex::Real;
168  Real h_xi, h_eta, h_zeta;
169 
170 #if 1
171  // *********************************************************
172  // Hi x-face
173  // *********************************************************
174  // On x-face
175  Real px_hi = (x(i+1,j,k) - x(i,j,k)) * dxinv;
176 
177  // On y-edges
178  Real pz_hi_md_hi = Real(0.5) * ( x(i+1,j ,k+1) + x(i ,j ,k+1)
179  -x(i+1,j ,k ) - x(i ,j ,k ) );
180  h_xi = Real(0.25) * ( zp(i+1,j ,k+1) - zp(i-1,j ,k+1)
181  +zp(i+1,j+1,k+1) - zp(i-1,j+1,k+1) ) * dxinv;
182  h_zeta = Real(0.25) * ( zp(i+1,j ,k+2) - zp(i+1,j ,k )
183  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
184  pz_hi_md_hi *= h_xi / h_zeta;
185 
186  // On y-edges
187  Real pz_hi_md_lo = Real(0.5) * ( x(i+1,j ,k ) + x(i ,j ,k )
188  -x(i+1,j ,k-1) - x(i ,j ,k-1) );
189  h_xi = Real(0.25) * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
190  +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
191  h_zeta = Real(0.25) * ( zp(i+1,j ,k+1) - zp(i+1,j ,k-1)
192  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
193  pz_hi_md_lo *= h_xi / h_zeta;
194 
195  // On x-face
196  px_hi -= Real(0.5) * ( pz_hi_md_hi + pz_hi_md_lo );
197 
198  // *********************************************************
199  // Lo x-face
200  // ********************************************************* // On x-face
201  Real px_lo = (x(i,j,k) - x(i-1,j,k)) * dxinv;
202 
203  // On y-edges
204  Real pz_lo_md_hi = Real(0.5) * ( x(i,j,k+1) + x(i-1,j,k+1)
205  -x(i,j,k ) - x(i-1,j,k ) );
206  h_xi = Real(0.25) * ( zp(i,j ,k+1) - zp(i-2,j ,k+1)
207  +zp(i,j+1,k+1) - zp(i-2,j+1,k+1) ) * dxinv;
208  h_zeta = Real(0.25) * ( zp(i,j ,k+2) - zp(i ,j ,k )
209  +zp(i,j+1,k+2) - zp(i ,j+1,k ) );
210  pz_lo_md_hi *= h_xi / h_zeta;
211 
212  // On y-edges
213  Real pz_lo_md_lo = Real(0.5) * ( x(i,j,k ) + x(i-1,j,k )
214  -x(i,j,k-1) - x(i-1,j,k-1) );
215  h_xi = Real(0.25) * ( zp(i,j ,k ) - zp(i-2,j ,k )
216  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
217  h_zeta = Real(0.25) * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
218  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
219  pz_lo_md_lo *= h_xi / h_zeta;
220 
221  // On x-face
222  px_lo -= Real(0.5) * ( pz_lo_md_hi + pz_lo_md_lo );
223 
224  // *********************************************************
225  // Hi y-face
226  // *********************************************************
227  // On y-face
228  Real py_hi = (x(i,j+1,k) - x(i,j,k)) * dyinv;
229 
230  // On x-edges
231  Real pz_md_hi_hi = Real(0.5) * ( x(i,j+1,k+1) + x(i,j,k+1)
232  -x(i,j+1,k ) - x(i,j,k ) );
233  h_eta = Real(0.25) * ( zp(i ,j+1,k+1) - zp(i ,j-1,k+1)
234  +zp(i+1,j+1,k+1) - zp(i+1,j-1,k+1) ) * dyinv;
235  h_zeta = Real(0.25) * ( zp(i ,j+1,k+2) - zp(i ,j+1,k )
236  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
237  pz_md_hi_hi *= h_eta/ h_zeta;
238 
239  // On x-edges
240  Real pz_md_hi_lo = Real(0.5) * ( x(i,j+1,k ) + x(i,j,k )
241  -x(i,j+1,k-1) - x(i,j,k-1) );
242  h_eta = Real(0.25) * ( zp(i ,j+1,k ) - zp(i ,j-1,k)
243  +zp(i+1,j+1,k ) - zp(i+1,j-1,k) ) * dyinv;
244  h_zeta = Real(0.25) * ( zp(i ,j+1,k+1) - zp(i ,j+1,k-1)
245  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
246  pz_md_hi_lo *= h_eta/ h_zeta;
247 
248  // On y-face
249  py_hi -= Real(0.5) * ( pz_md_hi_hi + pz_md_hi_lo );
250 
251  // *********************************************************
252  // Lo y-face
253  // *********************************************************
254  // On y-face
255  Real py_lo = (x(i,j,k) - x(i,j-1,k)) * dyinv;
256 
257  // On x-edges
258  Real pz_md_lo_hi = Real(0.5) * ( x(i ,j,k+1) + x(i ,j-1,k+1)
259  -x(i ,j,k ) - x(i ,j-1,k ) );
260  h_eta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j-2,k+1)
261  +zp(i+1,j,k+1) - zp(i+1,j-2,k+1) ) * dyinv;
262  h_zeta = Real(0.25) * ( zp(i ,j,k+2) - zp(i ,j ,k )
263  +zp(i+1,j,k+2) - zp(i+1,j ,k ) );
264  pz_md_lo_hi *= h_eta/ h_zeta;
265 
266  // On x-edges
267  Real pz_md_lo_lo = Real(0.5) * ( x(i ,j,k ) + x(i ,j-1,k )
268  -x(i ,j,k-1) - x(i ,j-1,k-1) );
269  h_eta = Real(0.25) * ( zp(i ,j,k ) - zp(i ,j-2,k )
270  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
271  h_zeta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
272  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
273  pz_md_lo_lo *= h_eta/ h_zeta;
274 
275  // On y-face
276  py_lo -= Real(0.5) * ( pz_md_lo_hi + pz_md_lo_lo );
277 
278  // *********************************************************
279  // Hi z-face
280  // *********************************************************
281  // On z-face
282  Real pz_hi = x(i,j,k+1) - x(i,j,k );
283  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))
284  -(zp(i,j,k ) + zp(i+1,j,k ) + zp(i,j+1,k ) + zp(i+1,j+1,k )) );
285  pz_hi *= hzeta_inv_on_zhi;
286 
287  // On corners
288  Real px_hi_md_hi = Real(0.5) * ( x(i+1,j,k+1) - x(i ,j ,k+1)
289  +x(i+1,j,k ) - x(i ,j ,k )) * dxinv;
290  Real px_lo_md_hi = Real(0.5) * ( x(i ,j,k+1) - x(i-1,j ,k+1)
291  +x(i ,j,k ) - x(i-1,j ,k )) * dxinv;
292  Real py_md_hi_hi = Real(0.5) * ( x(i,j+1,k+1) - x(i ,j ,k+1)
293  +x(i,j+1,k ) - x(i ,j ,k )) * dyinv;
294  Real py_md_lo_hi = Real(0.5) * ( x(i,j ,k+1) - x(i ,j-1,k+1)
295  +x(i,j ,k ) - x(i ,j-1,k )) * dyinv;
296 
297  // On z-face
298  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;
299  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;
300  //
301  // Note we do not need to recalculate pz_...hi here
302  //
303  pz_hi -= 0.5 * h_xi_on_zhi * ( (px_hi_md_hi + px_lo_md_hi) - (pz_hi_md_hi + pz_lo_md_hi) );
304  pz_hi -= 0.5 * h_eta_on_zhi * ( (py_md_hi_hi + py_md_lo_hi) - (pz_md_hi_hi + pz_md_lo_hi) );
305 
306  // *********************************************************
307  // Lo z-face
308  // *********************************************************
309  // On z-face
310  Real pz_lo = x(i,j,k ) - x(i,j,k-1);
311  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))
312  -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
313  pz_lo *= hzeta_inv_on_zlo;
314 
315  // On corners
316  Real px_hi_md_lo = Real(0.5) * ( x(i+1,j ,k ) - x(i ,j ,k )
317  +x(i+1,j ,k-1) - x(i ,j ,k-1)) * dxinv;
318  Real px_lo_md_lo = Real(0.5) * ( x(i ,j ,k ) - x(i-1,j ,k )
319  +x(i ,j ,k-1) - x(i-1,j ,k-1)) * dxinv;
320  Real py_md_hi_lo = Real(0.5) * ( x(i ,j+1,k ) - x(i ,j ,k )
321  +x(i ,j+1,k-1) - x(i ,j ,k-1)) * dyinv;
322  Real py_md_lo_lo = Real(0.5) * ( x(i ,j ,k ) - x(i ,j-1,k )
323  +x(i ,j ,k-1) - x(i ,j-1,k-1)) * dyinv;
324 
325  // On z-face
326  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;
327  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;
328  //
329  // Note we do not need to recalculate pz_...lo here
330  //
331  pz_lo -= 0.5 * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
332  pz_lo -= 0.5 * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
333 
334 #else
335 
336  // *********************************************************
337  // Version which calls flux routines
338  // *********************************************************
339  //
340  // This option uses calls to the flux routines so there is
341  // some duplicated computation
342  // This option should give the same answer as above
343  //
344  Real px_lo = -terrpoisson_flux_x(i ,j,k,x,zp,dxinv);
345  Real px_hi = -terrpoisson_flux_x(i+1,j,k,x,zp,dxinv);
346  Real py_lo = -terrpoisson_flux_y(i,j ,k,x,zp,dxinv);
347  Real py_hi = -terrpoisson_flux_y(i,j+1,k,x,zp,dyinv);
348  Real pz_lo = -terrpoisson_flux_z(i,j,k ,x,zp,dxinv,dyinv);
349  Real pz_hi = -terrpoisson_flux_z(i,j,k+1,x,zp,dxinv,dyinv);
350 #endif
351  //
352  // *********************************************************
353  // Adotx
354  // *********************************************************
355  if (k == 0) pz_lo = 0.0;
356 
357  y(i,j,k) = (ax(i+1,j,k)*px_hi - ax(i,j,k)*px_lo) * dxinv
358  +(ay(i,j+1,k)*py_hi - ay(i,j,k)*py_lo) * dyinv
359  +(az(i,j,k+1)*pz_hi - az(i,j,k)*pz_lo) * dzinv;
360  y(i,j,k) /= dJ(i,j,k);
361 }
362 #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_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:158
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