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