#include <AMReX_FArrayBox.H>
#include "ERF_NumericalConstants.H"

Include dependency graph for ERF_TerrainPoisson_3D_K.H:

This graph shows which files directly or indirectly include this file:

Enumerations
enum	TerrainPoissonFaceBit : int { tp_face_xlo = 1 , tp_face_xhi = 2 , tp_face_ylo = 4 , tp_face_yhi = 8 , tp_face_zlo = 16 , tp_face_zhi = 32 }

Functions
template<typename T >
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

template<typename T >
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

template<typename T >
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

template<typename T >
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

template<typename T >
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, amrex::Dim3 const &dlo=amrex::Dim3{0, 0, 0}, amrex::Dim3 const &dhi=amrex::Dim3{0, 0, 0}, int zf=0) noexcept

Enumeration Type Documentation

◆ TerrainPoissonFaceBit

enum TerrainPoissonFaceBit : int

Bits of the face mask taken by terrpoisson_adotx: one per face of the solve region whose flux is set to zero (TerrainPoisson masks the Neumann faces of its region).

Enumerator
tp_face_xlo
tp_face_xhi
tp_face_ylo
tp_face_yhi
tp_face_zlo
tp_face_zhi

                            : int {
     tp_face_xlo = 1, tp_face_xhi = 2,
     tp_face_ylo = 4, tp_face_yhi = 8,
     tp_face_zlo = 16, tp_face_zhi = 32
 };

Function Documentation

◆ terrpoisson_adotx()

template<typename T >

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,
		amrex::Dim3 const &	dlo = `amrex::Dim3{0,0,0}`,
		amrex::Dim3 const &	dhi = `amrex::Dim3{0,0,0}`,
		int	zf = `0`
	)

noexcept

Apply the terrain Poisson matrix to a cell-centered field.

Template Parameters

T	Real-valued array and spacing type

Parameters

i	x-index of the cell
j	y-index of the cell
k	z-index of the cell
y	Destination operator result
x	Source solution field
ax	Terrain metric coefficient on x-faces
ay	Terrain metric coefficient on y-faces
az	Terrain metric coefficient on z-faces
dJ	Cell-centered Jacobian determinant
zp	Node-centered physical height field
dxinv	Inverse x grid spacing
dyinv	Inverse y grid spacing
dzinv	Inverse z grid spacing
dlo	Low corner of the solve region (used with zf)
dhi	High corner of the solve region (used with zf)
zf	Faces of the region whose flux is zero (sum of TerrainPoissonFaceBit values; 0, the default, keeps every face flux)

A masked face is a Neumann face of the solve region: a face of the region inside the domain (a coarse/fine face of a refined level) or a Neumann (wall, inflow) face of the domain. The even mirror that TerrainPoisson::apply_bcs puts in the ghost cell zeroes the normal-gradient part of the flux through such a face but not the terrain cross term, so on sloping terrain the face would let correction through; with the face masked the whole flux through it is zero, as in AMReX's MLTerrainPoisson. Faces outside the mask are untouched, bit for bit.

                                                         {0,0,0},
                         amrex::Dim3 const& dhi = amrex::Dim3{0,0,0},
                         int zf = 0) noexcept
 {
     using amrex::Real;
     Real h_xi, h_eta, h_zeta;
  
 #if 1
     // *********************************************************
     // Hi x-face
     // *********************************************************
     // On x-face
     Real px_hi = (x(i+1,j,k) - x(i,j,k)) * dxinv;
  
     // On y-edges
     Real pz_hi_md_hi  = myhalf  * (  x(i+1,j  ,k+1) +  x(i  ,j  ,k+1)
                                     -x(i+1,j  ,k  ) -  x(i  ,j  ,k  ) );
          h_xi         = fourth * ( zp(i+2,j  ,k+1) - zp(i  ,j  ,k+1)
                                   +zp(i+2,j+1,k+1) - zp(i  ,j+1,k+1) ) * dxinv;
          h_zeta       = fourth * ( zp(i+1,j  ,k+2) - zp(i+1,j  ,k  )
                                   +zp(i+1,j+1,k+2) - zp(i+1,j+1,k  ) );
          pz_hi_md_hi *= h_xi / h_zeta;
  
     // On y-edges
     Real pz_hi_md_lo  = myhalf  *  ( x(i+1,j  ,k  ) +  x(i  ,j  ,k  )
                                     -x(i+1,j  ,k-1) -  x(i  ,j  ,k-1) );
          h_xi         = fourth * ( zp(i+2,j  ,k  ) - zp(i  ,j  ,k  )
                                   +zp(i+2,j+1,k  ) - zp(i  ,j+1,k  ) ) * dxinv;
          h_zeta       = fourth * ( zp(i+1,j  ,k+1) - zp(i+1,j  ,k-1)
                                   +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
          pz_hi_md_lo *= h_xi / h_zeta;
  
     // On x-face
     px_hi -= myhalf * ( pz_hi_md_hi + pz_hi_md_lo );
  
     // *********************************************************
     // Lo x-face
     // ********************************************************* // On x-face
     Real px_lo = (x(i,j,k) - x(i-1,j,k)) * dxinv;
  
     // On y-edges
     Real pz_lo_md_hi  = myhalf  * (  x(i,j,k+1) + x(i-1,j,k+1)
                                     -x(i,j,k  ) - x(i-1,j,k  ) );
          h_xi         = fourth * ( zp(i+1,j  ,k+1) - zp(i-1,j  ,k+1)
                                   +zp(i+1,j+1,k+1) - zp(i-1,j+1,k+1) ) * dxinv;
          h_zeta       = fourth * ( zp(i,j  ,k+2) - zp(i  ,j  ,k  )
                                   +zp(i,j+1,k+2) - zp(i  ,j+1,k  ) );
          pz_lo_md_hi *= h_xi / h_zeta;
  
     // On y-edges
     Real pz_lo_md_lo  = myhalf  * (  x(i,j,k  ) + x(i-1,j,k  )
                                     -x(i,j,k-1) - x(i-1,j,k-1) );
          h_xi         = fourth * ( zp(i+1,j  ,k  ) - zp(i-1,j  ,k  )
                                   +zp(i+1,j+1,k  ) - zp(i-1,j+1,k  ) ) * dxinv;
          h_zeta       = fourth * ( zp(i,j  ,k+1) - zp(i  ,j  ,k-1)
                                   +zp(i,j+1,k+1) - zp(i  ,j+1,k-1) );
          pz_lo_md_lo *= h_xi / h_zeta;
  
     // On x-face
     px_lo -= myhalf * ( pz_lo_md_hi + pz_lo_md_lo );
  
     // *********************************************************
     // Hi y-face
     // *********************************************************
     // On y-face
     Real py_hi = (x(i,j+1,k) - x(i,j,k)) * dyinv;
  
     // On x-edges
     Real pz_md_hi_hi  = myhalf  * (  x(i,j+1,k+1) + x(i,j,k+1)
                                        -x(i,j+1,k  ) - x(i,j,k  ) );
          h_eta        = fourth * ( zp(i  ,j+2,k+1) - zp(i  ,j  ,k+1)
                                       +zp(i+1,j+2,k+1) - zp(i+1,j  ,k+1) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j+1,k+2) - zp(i  ,j+1,k  )
                                       +zp(i+1,j+1,k+2) - zp(i+1,j+1,k  ) );
          pz_md_hi_hi *= h_eta/ h_zeta;
  
     // On x-edges
     Real pz_md_hi_lo  = myhalf  * (  x(i,j+1,k  ) + x(i,j,k  )
                                        -x(i,j+1,k-1) - x(i,j,k-1) );
          h_eta        = fourth * ( zp(i  ,j+2,k  ) - zp(i  ,j  ,k  )
                                       +zp(i+1,j+2,k  ) - zp(i+1,j  ,k  ) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j+1,k+1) - zp(i  ,j+1,k-1)
                                       +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
          pz_md_hi_lo *= h_eta/ h_zeta;
  
     // On y-face
     py_hi -= myhalf * ( pz_md_hi_hi + pz_md_hi_lo );
  
     // *********************************************************
     // Lo y-face
     // *********************************************************
     // On y-face
     Real py_lo = (x(i,j,k) - x(i,j-1,k)) * dyinv;
  
     // On x-edges
     Real pz_md_lo_hi  = myhalf  * (  x(i  ,j,k+1) +  x(i  ,j-1,k+1)
                                        -x(i  ,j,k  ) -  x(i  ,j-1,k  ) );
          h_eta        = fourth * ( zp(i  ,j+1,k+1) - zp(i  ,j-1,k+1)
                                       +zp(i+1,j+1,k+1) - zp(i+1,j-1,k+1) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j,k+2) - zp(i  ,j  ,k  )
                                       +zp(i+1,j,k+2) - zp(i+1,j  ,k  ) );
          pz_md_lo_hi *= h_eta/ h_zeta;
  
     // On x-edges
     Real pz_md_lo_lo  = myhalf  * (  x(i  ,j,k  ) +  x(i  ,j-1,k  )
                                        -x(i  ,j,k-1) -  x(i  ,j-1,k-1) );
          h_eta        = fourth * ( zp(i  ,j+1,k  ) - zp(i  ,j-1,k  )
                                       +zp(i+1,j+1,k  ) - zp(i+1,j-1,k  ) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j,k+1) - zp(i  ,j  ,k-1)
                                       +zp(i+1,j,k+1) - zp(i+1,j  ,k-1) );
          pz_md_lo_lo *= h_eta/ h_zeta;
  
     // On y-face
     py_lo -= myhalf * ( pz_md_lo_hi + pz_md_lo_lo );
  
     // *********************************************************
     // Hi z-face
     // *********************************************************
     // On z-face
     Real pz_hi    = x(i,j,k+1) - x(i,j,k  );
     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))
                                    -(zp(i,j,k  ) + zp(i+1,j,k  ) + zp(i,j+1,k  ) + zp(i+1,j+1,k  )) );
     pz_hi *= hzeta_inv_on_zhi;
  
     // On corners
     Real px_hi_md_hi = myhalf * ( x(i+1,j,k+1) - x(i  ,j  ,k+1)
                                     +x(i+1,j,k  ) - x(i  ,j  ,k  )) * dxinv;
     Real px_lo_md_hi = myhalf * ( x(i  ,j,k+1) - x(i-1,j  ,k+1)
                                     +x(i  ,j,k  ) - x(i-1,j  ,k  )) * dxinv;
     Real py_md_hi_hi = myhalf * ( x(i,j+1,k+1) - x(i  ,j  ,k+1)
                                     +x(i,j+1,k  ) - x(i  ,j  ,k  )) * dyinv;
     Real py_md_lo_hi = myhalf * ( x(i,j  ,k+1) - x(i  ,j-1,k+1)
                                     +x(i,j  ,k  ) - x(i  ,j-1,k  )) * dyinv;
  
     // On z-face
     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;
     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;
     //
     // Note we do not need to recalculate pz_...hi here
     //
     pz_hi -= myhalf * h_xi_on_zhi  * ( (px_hi_md_hi + px_lo_md_hi) - (pz_hi_md_hi + pz_lo_md_hi) );
     pz_hi -= myhalf * h_eta_on_zhi * ( (py_md_hi_hi + py_md_lo_hi) - (pz_md_hi_hi + pz_md_lo_hi) );
  
     // *********************************************************
     // Lo z-face
     // *********************************************************
     // On z-face
     Real pz_lo = x(i,j,k  ) - x(i,j,k-1);
     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))
                                    -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
     pz_lo *= hzeta_inv_on_zlo;
  
     // On corners
     Real px_hi_md_lo = myhalf * ( x(i+1,j  ,k  ) - x(i  ,j  ,k  )
                                     +x(i+1,j  ,k-1) - x(i  ,j  ,k-1)) * dxinv;
     Real px_lo_md_lo = myhalf * ( x(i  ,j  ,k  ) - x(i-1,j  ,k  )
                                     +x(i  ,j  ,k-1) - x(i-1,j  ,k-1)) * dxinv;
     Real py_md_hi_lo = myhalf * ( x(i  ,j+1,k  ) - x(i  ,j  ,k  )
                                     +x(i  ,j+1,k-1) - x(i  ,j  ,k-1)) * dyinv;
     Real py_md_lo_lo = myhalf * ( x(i  ,j  ,k  ) - x(i  ,j-1,k  )
                                     +x(i  ,j  ,k-1) - x(i  ,j-1,k-1)) * dyinv;
  
     // On z-face
     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;
     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;
     //
     // Note we do not need to recalculate pz_...lo here
     //
     pz_lo -= myhalf * h_xi_on_zlo  * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
     pz_lo -= myhalf * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
  
 #else
  
     // *********************************************************
     // Version which calls flux routines
     // *********************************************************
     //
     // This option uses calls to the flux routines so there is
     //      some duplicated computation
     // This option should give the same answer as above
     //
         Real px_lo = -terrpoisson_flux_x(i  ,j,k,x,zp,dxinv);
         Real px_hi = -terrpoisson_flux_x(i+1,j,k,x,zp,dxinv);
         Real py_lo = -terrpoisson_flux_y(i,j  ,k,x,zp,dyinv);
         Real py_hi = -terrpoisson_flux_y(i,j+1,k,x,zp,dyinv);
         Real pz_lo = -terrpoisson_flux_z(i,j,k  ,x,zp,dxinv,dyinv);
         Real pz_hi = -terrpoisson_flux_z(i,j,k+1,x,zp,dxinv,dyinv);
 #endif
     //
     // *********************************************************
     // Adotx
     // *********************************************************
     if (k == 0) pz_lo = zero;
     if ((zf & tp_face_xlo) && i == dlo.x) px_lo = zero;
     if ((zf & tp_face_xhi) && i == dhi.x) px_hi = zero;
     if ((zf & tp_face_ylo) && j == dlo.y) py_lo = zero;
     if ((zf & tp_face_yhi) && j == dhi.y) py_hi = zero;
     if ((zf & tp_face_zlo) && k == dlo.z) pz_lo = zero;
     if ((zf & tp_face_zhi) && k == dhi.z) pz_hi = zero;
  
     y(i,j,k) =  (ax(i+1,j,k)*px_hi - ax(i,j,k)*px_lo) * dxinv
                +(ay(i,j+1,k)*py_hi - ay(i,j,k)*py_lo) * dyinv
                +(az(i,j,k+1)*pz_hi - az(i,j,k)*pz_lo) * dzinv;
     y(i,j,k) /=  dJ(i,j,k);
 }

◆ terrpoisson_flux_x()

template<typename T >

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

Compute the x-face terrain Poisson flux.

Template Parameters

T	Real-valued array and spacing type

Parameters

i	x-index of the x-face
j	y-index of the x-face
k	z-index of the x-face
sol	Cell-centered solution field
zp	Node-centered physical height field
dxinv	Inverse x grid spacing

Returns: x-face flux for the terrain Poisson operator

 {
     using amrex::Real;
  
     Real h_xi, h_zeta;
  
     // On x-face
     Real px_lo = (sol(i,j,k) - sol(i-1,j,k)) * dxinv;
  
     // On y-edges
     Real pz_lo_md_hi  = myhalf * ( sol(i,j,k+1) + sol(i-1,j,k+1)
                                      -sol(i,j,k  ) - sol(i-1,j,k  ) );
          h_xi         = fourth * ( zp(i+1,j  ,k+1) - zp(i-1,j  ,k+1)
                                       +zp(i+1,j+1,k+1) - zp(i-1,j+1,k+1) ) * dxinv;
          h_zeta       = fourth * ( zp(i,j  ,k+2) - zp(i  ,j  ,k  )
                                       +zp(i,j+1,k+2) - zp(i  ,j+1,k  ) );
          pz_lo_md_hi *= h_xi / h_zeta;
  
     // On y-edges
     Real pz_lo_md_lo  = myhalf * ( sol(i,j,k  ) + sol(i-1,j,k  )
                                      -sol(i,j,k-1) - sol(i-1,j,k-1) );
  
          h_xi         = fourth * ( zp(i+1,j  ,k  ) - zp(i-1,j  ,k  )
                                       +zp(i+1,j+1,k  ) - zp(i-1,j+1,k  ) ) * dxinv;
          h_zeta       = fourth * ( zp(i,j  ,k+1) - zp(i  ,j  ,k-1)
                                       +zp(i,j+1,k+1) - zp(i  ,j+1,k-1) );
          pz_lo_md_lo *= h_xi / h_zeta;
  
     // On x-face
     px_lo -= myhalf * ( pz_lo_md_hi + pz_lo_md_lo );
  
     return -px_lo;
 }

◆ terrpoisson_flux_y()

template<typename T >

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

Compute the y-face terrain Poisson flux.

Template Parameters

T	Real-valued array and spacing type

Parameters

i	x-index of the y-face
j	y-index of the y-face
k	z-index of the y-face
sol	Cell-centered solution field
zp	Node-centered physical height field
dyinv	Inverse y grid spacing

Returns: y-face flux for the terrain Poisson operator

 {
     using amrex::Real;
  
     Real h_eta, h_zeta;
  
     Real py_lo = (sol(i,j,k) - sol(i,j-1,k)) * dyinv;
  
     // On x-edges
     Real pz_md_lo_hi  = myhalf * ( sol(i,j,k+1) + sol(i,j-1,k+1)
                                      -sol(i,j,k  ) - sol(i,j-1,k  ) );
          h_eta        = fourth * ( zp(i  ,j+1,k+1) - zp(i  ,j-1,k+1)
                                       +zp(i+1,j+1,k+1) - zp(i+1,j-1,k+1) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j,k+2) - zp(i  ,j  ,k  )
                                       +zp(i+1,j,k+2) - zp(i+1,j  ,k  ) );
          pz_md_lo_hi *= h_eta/ h_zeta;
  
     // On x-edges
     Real pz_md_lo_lo  = myhalf * ( sol(i,j,k  ) + sol(i,j-1,k  )
                                      -sol(i,j,k-1) - sol(i,j-1,k-1) );
  
          h_eta        = fourth * ( zp(i  ,j+1,k  ) - zp(i  ,j-1,k  )
                                       +zp(i+1,j+1,k  ) - zp(i+1,j-1,k  ) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j,k+1) - zp(i  ,j  ,k-1)
                                       +zp(i+1,j,k+1) - zp(i+1,j  ,k-1) );
          pz_md_lo_lo *= h_eta/ h_zeta;
  
     // On y-face
     py_lo -= myhalf * ( pz_md_lo_hi + pz_md_lo_lo );
     return -py_lo;
 }

◆ terrpoisson_flux_z()

template<typename T >

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

Compute the z-face terrain Poisson flux.

Template Parameters

T	Real-valued array and spacing type

Parameters

i	x-index of the z-face
j	y-index of the z-face
k	z-index of the z-face
sol	Cell-centered solution field
zp	Node-centered physical height field
dxinv	Inverse x grid spacing
dyinv	Inverse y grid spacing

Returns: z-face flux for the terrain Poisson operator

 {
     using amrex::Real;
  
     Real h_xi, h_eta, h_zeta;
  
     // On z-face
     Real pz_lo = (sol(i,j,k  ) - sol(i,j,k-1));
     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))
                                    -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
     pz_lo *= hzeta_inv_on_zlo;
  
     // On corners
     Real px_hi_md_lo = myhalf * ( sol(i+1,j  ,k  ) - sol(i  ,j  ,k  )
                                     +sol(i+1,j  ,k-1) - sol(i  ,j  ,k-1)) * dxinv;
     Real px_lo_md_lo = myhalf * ( sol(i  ,j  ,k  ) - sol(i-1,j  ,k  )
                                     +sol(i  ,j  ,k-1) - sol(i-1,j  ,k-1)) * dxinv;
     Real py_md_hi_lo = myhalf * ( sol(i  ,j+1,k  ) - sol(i  ,j  ,k  )
                                     +sol(i  ,j+1,k-1) - sol(i  ,j  ,k-1)) * dyinv;
     Real py_md_lo_lo = myhalf * ( sol(i  ,j  ,k  ) - sol(i  ,j-1,k  )
                                     +sol(i  ,j  ,k-1) - sol(i  ,j-1,k-1)) * dyinv;
  
     // On y-edges
     Real pz_hi_md_lo  = myhalf * ( sol(i+1,j,k  ) + sol(i,j,k  )
                                      -sol(i+1,j,k-1) - sol(i,j,k-1) );
          h_xi         = fourth * ( zp(i+2,j  ,k  ) - zp(i  ,j  ,k  )
                                       +zp(i+2,j+1,k  ) - zp(i  ,j+1,k  ) ) * dxinv;
          h_zeta       = fourth * ( zp(i+1,j  ,k+1) - zp(i+1,j  ,k-1)
                                       +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
          pz_hi_md_lo *= h_xi / h_zeta;
  
     // On y-edges
     Real pz_lo_md_lo  = myhalf * ( sol(i,j,k  ) + sol(i-1,j,k  )
                                      -sol(i,j,k-1) - sol(i-1,j,k-1) );
          h_xi         = fourth * ( zp(i+1,j  ,k  ) - zp(i-1,j  ,k  )
                                       +zp(i+1,j+1,k  ) - zp(i-1,j+1,k  ) ) * dxinv;
          h_zeta       = fourth * ( zp(i,j  ,k+1) - zp(i  ,j  ,k-1)
                                       +zp(i,j+1,k+1) - zp(i  ,j+1,k-1) );
          pz_lo_md_lo *= h_xi / h_zeta;
  
     // On x-edges
     Real pz_md_hi_lo  = myhalf * ( sol(i,j+1,k  ) + sol(i,j,k  )
                                      -sol(i,j+1,k-1) - sol(i,j,k-1) );
          h_eta        = fourth * ( zp(i  ,j+2,k  ) - zp(i  ,j  ,k  )
                                       +zp(i+1,j+2,k  ) - zp(i+1,j  ,k  ) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j+1,k+1) - zp(i  ,j+1,k-1)
                                       +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
          pz_md_hi_lo *= h_eta/ h_zeta;
  
     // On x-edges
     Real pz_md_lo_lo  = myhalf * ( sol(i,j,k  ) + sol(i,j-1,k  )
                                      -sol(i,j,k-1) - sol(i,j-1,k-1) );
          h_eta        = fourth * ( zp(i  ,j+1,k  ) - zp(i  ,j-1,k  )
                                       +zp(i+1,j+1,k  ) - zp(i+1,j-1,k  ) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j,k+1) - zp(i  ,j  ,k-1)
                                       +zp(i+1,j,k+1) - zp(i+1,j  ,k-1) );
          pz_md_lo_lo *= h_eta/ h_zeta;
  
     // On z-face
     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;
     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;
  
     pz_lo -= myhalf * h_xi_on_zlo  * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
     pz_lo -= myhalf * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
  
     if (k == 0) pz_lo = zero;
  
     return -pz_lo;
 }

◆ terrpoisson_flux_zlo_dir()

template<typename T >

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

Simplified gradient calc with dirichlet BC (phi_zlo==0)

compare with terrpoisson_flux_z

Template Parameters

T	Real-valued array and spacing type

Parameters

i	x-index of the z-face
j	y-index of the z-face
k	z-index of the z-face
sol	Cell-centered solution field
zp	Node-centered physical height field
dxinv	Inverse x grid spacing
dyinv	Inverse y grid spacing

Returns: z-face flux with a low-z Dirichlet value

 {
     using amrex::Real;
  
     Real h_xi, h_eta, h_zeta;
  
     // On z-face, one-sided
     Real pz_lo = sol(i,j,k);
     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))
                                    -(zp(i,j,k  ) + zp(i+1,j,k  ) + zp(i,j+1,k  ) + zp(i+1,j+1,k  )) );
     pz_lo *= hzeta_inv_on_zlo;
  
     // On corners
     constexpr Real px_hi_md_lo = zero;
     constexpr Real px_lo_md_lo = zero;
     constexpr Real py_md_hi_lo = zero;
     constexpr Real py_md_lo_lo = zero;
  
     // On y-edges, one-sided
     Real pz_hi_md_lo  = myhalf * ( sol(i+1,j,k  ) + sol(i,j,k  ));
          h_xi         = fourth * ( zp(i+2,j  ,k  ) - zp(i  ,j  ,k  )
                                       +zp(i+2,j+1,k  ) - zp(i  ,j+1,k  ) ) * dxinv;
          h_zeta       = fourth * ( zp(i+1,j  ,k+1) - zp(i+1,j  ,k  )
                                       +zp(i+1,j+1,k+1) - zp(i+1,j+1,k  ) );
          pz_hi_md_lo *= h_xi / h_zeta;
  
     // On y-edges, one-sided
     Real pz_lo_md_lo  = myhalf * ( sol(i,j,k  ) + sol(i-1,j,k  ));
          h_xi         = fourth * ( zp(i+1,j  ,k  ) - zp(i-1,j  ,k  )
                                       +zp(i+1,j+1,k  ) - zp(i-1,j+1,k  ) ) * dxinv;
          h_zeta       = fourth * ( zp(i,j  ,k+1) - zp(i  ,j  ,k  )
                                       +zp(i,j+1,k+1) - zp(i  ,j+1,k  ) );
          pz_lo_md_lo *= h_xi / h_zeta;
  
     // On x-edges, one-sided
     Real pz_md_hi_lo  = myhalf * ( sol(i,j+1,k  ) + sol(i,j,k  ));
          h_eta        = fourth * ( zp(i  ,j+2,k  ) - zp(i  ,j  ,k  )
                                       +zp(i+1,j+2,k  ) - zp(i+1,j  ,k  ) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j+1,k+1) - zp(i  ,j+1,k  )
                                       +zp(i+1,j+1,k+1) - zp(i+1,j+1,k  ) );
          pz_md_hi_lo *= h_eta/ h_zeta;
  
     // On x-edges, one-sided
     Real pz_md_lo_lo  = myhalf * ( sol(i,j,k  ) + sol(i,j-1,k  ));
          h_eta        = fourth * ( zp(i  ,j+1,k  ) - zp(i  ,j-1,k  )
                                       +zp(i+1,j+1,k  ) - zp(i+1,j-1,k  ) ) * dyinv;
          h_zeta       = fourth * ( zp(i  ,j,k+1) - zp(i  ,j  ,k  )
                                       +zp(i+1,j,k+1) - zp(i+1,j  ,k  ) );
          pz_md_lo_lo *= h_eta/ h_zeta;
  
     // On z-face
     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;
     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;
  
     pz_lo -= myhalf * h_xi_on_zlo  * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
     pz_lo -= myhalf * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
  
     return -pz_lo;
 }

Enumerations

Functions

Enumeration Type Documentation

◆ TerrainPoissonFaceBit

Function Documentation

◆ terrpoisson_adotx()

◆ terrpoisson_flux_x()

◆ terrpoisson_flux_y()

◆ terrpoisson_flux_z()

◆ terrpoisson_flux_zlo_dir()