1 #ifndef ERF_TERRAINPOISSON_3D_K_H_
2 #define ERF_TERRAINPOISSON_3D_K_H_
4 #include <AMReX_FArrayBox.H>
7 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
9 amrex::Array4<T const>
const& sol,
10 amrex::Array4<T const>
const& zp,
18 Real px_lo = (sol(i,j,k) - sol(i-1,j,k)) * dxinv;
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;
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) );
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;
40 px_lo -= Real(0.5) * ( pz_lo_md_hi + pz_lo_md_lo );
45 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
47 amrex::Array4<T const>
const& sol,
48 amrex::Array4<T const>
const& zp,
55 Real py_lo = (sol(i,j,k) - sol(i,j-1,k)) * dyinv;
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;
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) );
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;
77 py_lo -= Real(0.5) * ( pz_md_lo_hi + pz_md_lo_lo );
83 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
85 amrex::Array4<T const>
const& sol,
86 amrex::Array4<T const>
const& zp,
87 T dxinv,
T dyinv) noexcept
91 Real h_xi, h_eta, h_zeta;
94 Real pz_lo = (sol(i,j,k ) - sol(i,j,k-1));
95 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))
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;
100 Real px_hi_md_lo = Real(0.5) * ( 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 = Real(0.5) * ( 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 = Real(0.5) * ( 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 = Real(0.5) * ( sol(i ,j ,k ) - sol(i ,j-1,k )
107 +sol(i ,j ,k-1) - sol(i ,j-1,k-1)) * dyinv;
110 Real pz_hi_md_lo = Real(0.5) * ( sol(i+1,j,k ) + sol(i,j,k )
111 -sol(i+1,j,k-1) - sol(i,j,k-1) );
112 h_xi = Real(0.25) * ( 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 = Real(0.25) * ( 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;
119 Real pz_lo_md_lo = Real(0.5) * ( sol(i,j,k ) + sol(i-1,j,k )
120 -sol(i,j,k-1) - sol(i-1,j,k-1) );
121 h_xi = Real(0.25) * ( 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 = Real(0.25) * ( 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;
128 Real pz_md_hi_lo = Real(0.5) * ( sol(i,j+1,k ) + sol(i,j,k )
129 -sol(i,j+1,k-1) - sol(i,j,k-1) );
130 h_eta = Real(0.25) * ( 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 = Real(0.25) * ( 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;
137 Real pz_md_lo_lo = Real(0.5) * ( sol(i,j,k ) + sol(i,j-1,k )
138 -sol(i,j,k-1) - sol(i,j-1,k-1) );
139 h_eta = Real(0.25) * ( 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 = Real(0.25) * ( 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;
146 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;
147 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;
149 pz_lo -= 0.5 * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
150 pz_lo -= 0.5 * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
155 template <
typename T>
156 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
158 amrex::Array4<T const>
const& x,
159 amrex::Array4<T const>
const& zp,
160 T dxinv,
T dyinv) noexcept
163 Real h_xi, h_eta, h_zeta;
169 Real px_hi = (
x(i+1,j,k) -
x(i,j,k)) * dxinv;
172 Real pz_hi_md_hi = Real(0.5) * (
x(i+1,j ,k+1) +
x(i ,j ,k+1)
173 -
x(i+1,j ,k ) -
x(i ,j ,k ) );
174 h_xi = Real(0.25) * ( zp(i+1,j ,k+1) - zp(i-1,j ,k+1)
175 +zp(i+1,j+1,k+1) - zp(i-1,j+1,k+1) ) * dxinv;
176 h_zeta = Real(0.25) * ( zp(i+1,j ,k+2) - zp(i+1,j ,k )
177 +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
178 pz_hi_md_hi *= h_xi / h_zeta;
181 Real pz_hi_md_lo = Real(0.5) * (
x(i+1,j ,k ) +
x(i ,j ,k )
182 -
x(i+1,j ,k-1) -
x(i ,j ,k-1) );
183 h_xi = Real(0.25) * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
184 +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
185 h_zeta = Real(0.25) * ( zp(i+1,j ,k+1) - zp(i+1,j ,k-1)
186 +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
187 pz_hi_md_lo *= h_xi / h_zeta;
190 px_hi -= Real(0.5) * ( pz_hi_md_hi + pz_hi_md_lo );
195 Real px_lo = (
x(i,j,k) -
x(i-1,j,k)) * dxinv;
198 Real pz_lo_md_hi = Real(0.5) * (
x(i,j,k+1) +
x(i-1,j,k+1)
199 -
x(i,j,k ) -
x(i-1,j,k ) );
200 h_xi = Real(0.25) * ( zp(i,j ,k+1) - zp(i-2,j ,k+1)
201 +zp(i,j+1,k+1) - zp(i-2,j+1,k+1) ) * dxinv;
202 h_zeta = Real(0.25) * ( zp(i,j ,k+2) - zp(i ,j ,k )
203 +zp(i,j+1,k+2) - zp(i ,j+1,k ) );
204 pz_lo_md_hi *= h_xi / h_zeta;
207 Real pz_lo_md_lo = Real(0.5) * (
x(i,j,k ) +
x(i-1,j,k )
208 -
x(i,j,k-1) -
x(i-1,j,k-1) );
209 h_xi = Real(0.25) * ( zp(i,j ,k ) - zp(i-2,j ,k )
210 +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
211 h_zeta = Real(0.25) * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
212 +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
213 pz_lo_md_lo *= h_xi / h_zeta;
216 px_lo -= Real(0.5) * ( pz_lo_md_hi + pz_lo_md_lo );
222 Real py_hi = (
x(i,j+1,k) -
x(i,j,k)) * dyinv;
225 Real pz_md_hi_hi = Real(0.5) * (
x(i,j+1,k+1) +
x(i,j,k+1)
226 -
x(i,j+1,k ) -
x(i,j,k ) );
227 h_eta = Real(0.25) * ( zp(i ,j+1,k+1) - zp(i ,j-1,k+1)
228 +zp(i+1,j+1,k+1) - zp(i+1,j-1,k+1) ) * dyinv;
229 h_zeta = Real(0.25) * ( zp(i ,j+1,k+2) - zp(i ,j+1,k )
230 +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
231 pz_md_hi_hi *= h_eta/ h_zeta;
234 Real pz_md_hi_lo = Real(0.5) * (
x(i,j+1,k ) +
x(i,j,k )
235 -
x(i,j+1,k-1) -
x(i,j,k-1) );
236 h_eta = Real(0.25) * ( zp(i ,j+1,k ) - zp(i ,j-1,k)
237 +zp(i+1,j+1,k ) - zp(i+1,j-1,k) ) * dyinv;
238 h_zeta = Real(0.25) * ( zp(i ,j+1,k+1) - zp(i ,j+1,k-1)
239 +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
240 pz_md_hi_lo *= h_eta/ h_zeta;
243 py_hi -= Real(0.5) * ( pz_md_hi_hi + pz_md_hi_lo );
249 Real py_lo = (
x(i,j,k) -
x(i,j-1,k)) * dyinv;
252 Real pz_md_lo_hi = Real(0.5) * (
x(i ,j,k+1) +
x(i ,j-1,k+1)
253 -
x(i ,j,k ) -
x(i ,j-1,k ) );
254 h_eta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j-2,k+1)
255 +zp(i+1,j,k+1) - zp(i+1,j-2,k+1) ) * dyinv;
256 h_zeta = Real(0.25) * ( zp(i ,j,k+2) - zp(i ,j ,k )
257 +zp(i+1,j,k+2) - zp(i+1,j ,k ) );
258 pz_md_lo_hi *= h_eta/ h_zeta;
261 Real pz_md_lo_lo = Real(0.5) * (
x(i ,j,k ) +
x(i ,j-1,k )
262 -
x(i ,j,k-1) -
x(i ,j-1,k-1) );
263 h_eta = Real(0.25) * ( zp(i ,j,k ) - zp(i ,j-2,k )
264 +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
265 h_zeta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
266 +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
267 pz_md_lo_lo *= h_eta/ h_zeta;
270 py_lo -= Real(0.5) * ( pz_md_lo_hi + pz_md_lo_lo );
276 Real pz_hi =
x(i,j,k+1) -
x(i,j,k );
277 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))
278 -(zp(i,j,k ) + zp(i+1,j,k ) + zp(i,j+1,k ) + zp(i+1,j+1,k )) );
279 pz_hi *= hzeta_inv_on_zhi;
282 Real px_hi_md_hi = Real(0.5) * (
x(i+1,j,k+1) -
x(i ,j ,k+1)
283 +
x(i+1,j,k ) -
x(i ,j ,k )) * dxinv;
284 Real px_lo_md_hi = Real(0.5) * (
x(i ,j,k+1) -
x(i-1,j ,k+1)
285 +
x(i ,j,k ) -
x(i-1,j ,k )) * dxinv;
286 Real py_md_hi_hi = Real(0.5) * (
x(i,j+1,k+1) -
x(i ,j ,k+1)
287 +
x(i,j+1,k ) -
x(i ,j ,k )) * dyinv;
288 Real py_md_lo_hi = Real(0.5) * (
x(i,j ,k+1) -
x(i ,j-1,k+1)
289 +
x(i,j ,k ) -
x(i ,j-1,k )) * dyinv;
292 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;
293 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;
297 pz_hi -= 0.5 * h_xi_on_zhi * ( (px_hi_md_hi + px_lo_md_hi) - (pz_hi_md_hi + pz_lo_md_hi) );
298 pz_hi -= 0.5 * h_eta_on_zhi * ( (py_md_hi_hi + py_md_lo_hi) - (pz_md_hi_hi + pz_md_lo_hi) );
304 Real pz_lo =
x(i,j,k ) -
x(i,j,k-1);
305 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))
306 -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
307 pz_lo *= hzeta_inv_on_zlo;
310 Real px_hi_md_lo = Real(0.5) * (
x(i+1,j ,k ) -
x(i ,j ,k )
311 +
x(i+1,j ,k-1) -
x(i ,j ,k-1)) * dxinv;
312 Real px_lo_md_lo = Real(0.5) * (
x(i ,j ,k ) -
x(i-1,j ,k )
313 +
x(i ,j ,k-1) -
x(i-1,j ,k-1)) * dxinv;
314 Real py_md_hi_lo = Real(0.5) * (
x(i ,j+1,k ) -
x(i ,j ,k )
315 +
x(i ,j+1,k-1) -
x(i ,j ,k-1)) * dyinv;
316 Real py_md_lo_lo = Real(0.5) * (
x(i ,j ,k ) -
x(i ,j-1,k )
317 +
x(i ,j ,k-1) -
x(i ,j-1,k-1)) * dyinv;
320 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;
321 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;
325 pz_lo -= 0.5 * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
326 pz_lo -= 0.5 * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
346 Real invdJ = 4.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)
347 -zp(i,j,k ) - zp(i+1,j,k ) - zp(i,j+1,k ) - zp(i+1,j+1,k ) );
349 Real ax_lo = .5 * (zp(i ,j,k+1) + zp(i ,j+1,k+1) - zp(i ,j,k) - zp(i ,j+1,k));
350 Real ax_hi = .5 * (zp(i+1,j,k+1) + zp(i+1,j+1,k+1) - zp(i+1,j,k) - zp(i+1,j+1,k));
351 Real ay_lo = .5 * (zp(i,j ,k+1) + zp(i+1,j ,k+1) - zp(i,j ,k) - zp(i+1,j ,k));
352 Real ay_hi = .5 * (zp(i,j+1,k+1) + zp(i+1,j+1,k+1) - zp(i,j+1,k) - zp(i+1,j+1,k));
354 y(i,j,k) = ( (ax_hi*px_hi - ax_lo*px_lo) * dxinv
355 +(ay_hi*py_hi - ay_lo*py_lo) * dyinv
356 +( pz_hi - pz_lo) ) * invdJ;
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 &zp, T dxinv, T dyinv) noexcept
Definition: ERF_TerrainPoisson_3D_K.H:157
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
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
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
@ T
Definition: ERF_IndexDefines.H:99