238 Real h_xi, h_eta, h_zeta;
245 Real px_hi = (
x(i+1,j,k) -
x(i,j,k)) * dxinv;
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;
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;
266 px_hi -=
myhalf * ( pz_hi_md_hi + pz_hi_md_lo );
271 Real px_lo = (
x(i,j,k) -
x(i-1,j,k)) * dxinv;
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;
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;
292 px_lo -=
myhalf * ( pz_lo_md_hi + pz_lo_md_lo );
298 Real py_hi = (
x(i,j+1,k) -
x(i,j,k)) * dyinv;
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;
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;
319 py_hi -=
myhalf * ( pz_md_hi_hi + pz_md_hi_lo );
325 Real py_lo = (
x(i,j,k) -
x(i,j-1,k)) * dyinv;
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;
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;
346 py_lo -=
myhalf * ( pz_md_lo_hi + pz_md_lo_lo );
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;
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;
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;
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) );
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;
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;
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;
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) );
425 if (k == 0) pz_lo =
zero;
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);
constexpr amrex::Real fourth
Definition: ERF_Constants.H:12
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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
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:47