ERFPhysBCFunct_u Class Reference

#include <ERF_PhysBCFunct.H>

Collaboration diagram for ERFPhysBCFunct_u:

Public Member Functions
	ERFPhysBCFunct_u (const int lev, const amrex::Geometry &geom, const amrex::Vector< amrex::BCRec > &domain_bcs_type, const amrex::Gpu::DeviceVector< amrex::BCRec > &domain_bcs_type_d, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > bc_extdir_vals, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > bc_neumann_vals, const TerrainType &terrain_type, std::unique_ptr< amrex::MultiFab > &z_phys_nd, const bool use_real_bcs, amrex::Real *u_bc_data)
	~ERFPhysBCFunct_u ()
void	operator() (amrex::MultiFab &mf, amrex::MultiFab &xvel, amrex::MultiFab &yvel, amrex::IntVect const &nghost, const double time, int bccomp, bool do_fb)
void	impose_lateral_xvel_bcs (const amrex::Array4< amrex::Real > &dest_arr, const amrex::Array4< amrex::Real const > &xvel_arr, const amrex::Array4< amrex::Real const > &yvel_arr, const amrex::Box &bx, const amrex::Box &domain, int bccomp, const double time)
void	impose_vertical_xvel_bcs (const amrex::Array4< amrex::Real > &dest_arr, const amrex::Box &bx, const amrex::Box &domain, const amrex::Array4< amrex::Real const > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > dxInv, int bccomp, const double time)

Private Attributes
int	m_lev
amrex::Geometry	m_geom
amrex::Vector< amrex::BCRec >	m_domain_bcs_type
amrex::Gpu::DeviceVector< amrex::BCRec >	m_domain_bcs_type_d
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max >	m_bc_extdir_vals
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max >	m_bc_neumann_vals
TerrainType	m_terrain_type
amrex::MultiFab *	m_z_phys_nd
bool	m_use_real_bcs
amrex::Real *	m_u_bc_data
BelowGroundGhostSync	m_below_ground_sync

Constructor & Destructor Documentation

ERFPhysBCFunct_u()

ERFPhysBCFunct_u::ERFPhysBCFunct_u ( const int  lev, const amrex::Geometry &  geom, const amrex::Vector< amrex::BCRec > &  domain_bcs_type, const amrex::Gpu::DeviceVector< amrex::BCRec > &  domain_bcs_type_d, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max >  bc_extdir_vals, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max >  bc_neumann_vals, const TerrainType &  terrain_type, std::unique_ptr< amrex::MultiFab > &  z_phys_nd, const bool  use_real_bcs, amrex::Real *  u_bc_data  )

inline

Construct the x-velocity physical boundary-condition functor.

Parameters

[in]	lev	level index
[in]	geom	geometry for the level
[in]	domain_bcs_type	host boundary-condition records
[in]	domain_bcs_type_d	device boundary-condition records
[in]	bc_extdir_vals	external Dirichlet values for boundary faces
[in]	bc_neumann_vals	Neumann values for boundary faces
[in]	terrain_type	terrain type for the level
[in]	z_phys_nd	nodal physical height MultiFab
[in]	use_real_bcs	whether real-data boundary conditions are active
[in]	u_bc_data	optional x-velocity boundary data

        : m_lev(lev), m_geom(geom),
          m_domain_bcs_type(domain_bcs_type),
          m_domain_bcs_type_d(domain_bcs_type_d),
          m_bc_extdir_vals(bc_extdir_vals),
          m_bc_neumann_vals(bc_neumann_vals),
          m_terrain_type(terrain_type),
          m_z_phys_nd(z_phys_nd.get()),
          m_use_real_bcs(use_real_bcs),
          m_u_bc_data(u_bc_data)
    { }

~ERFPhysBCFunct_u()

ERFPhysBCFunct_u::~ERFPhysBCFunct_u ( )

inline

Destroy the x-velocity physical boundary-condition functor.

{}

Member Function Documentation

impose_lateral_xvel_bcs()

void ERFPhysBCFunct_u::impose_lateral_xvel_bcs	(	const amrex::Array4< amrex::Real > &	dest_arr,
		const amrex::Array4< amrex::Real const > &	xvel_arr,
		const amrex::Array4< amrex::Real const > &	yvel_arr,
		const amrex::Box &	bx,
		const amrex::Box &	domain,
		int	bccomp,
		const double	time
	)

Impose lateral boundary conditions on x velocity.

Parameters

[in,out]	dest_arr	x-velocity data to be filled
[in]	xvel_arr	x-velocity used to determine upwind inflow
[in]	yvel_arr	y-velocity used to determine upwind inflow
[in]	bx	box holding data to be filled
[in]	domain	simulation domain
[in]	bccomp	boundary-condition component
[in]	time	time at which data should be filled

Impose lateral boundary conditions on x-component of velocity

Parameters

[in,out]	dest_arr	Array4 of the quantity to be filled
[in]	xvel_arr	x-velocity used to determine upwind Dirichlet inflow
[in]	yvel_arr	y-velocity used to determine upwind Dirichlet inflow
[in]	bx	box associated with this data
[in]	domain	computational domain
[in]	bccomp	index into m_domain_bcs_type
[in]	time	time at which the data should be filled

{
    BL_PROFILE_VAR("impose_lateral_xvel_bcs()",impose_lateral_xvel_bcs);
    const auto& dom_lo = lbound(domain);
    const auto& dom_hi = ubound(domain);
 
    //
    // An ext_dir_upwind face decides inflow against the face velocity, but the box we
    // loop over does not always have a face velocity at every one of its cells: the
    // state carries one more lateral ghost than the velocities do (ngrow_state =
    // ComputeGhostCells()+1 against ngrow_vels = ComputeGhostCells()), and a velocity
    // box is nodal in its own direction where the array being read is cell centered.
    // Either way the outermost ghost cells have no velocity of their own, so sample
    // the nearest one that exists. Cells already inside the velocity array are
    // unaffected; only that outermost layer changes, and it previously read out of
    // bounds.
    //
    const auto& xvel_lo = lbound(xvel_arr);
    const auto& xvel_hi = ubound(xvel_arr);
    const auto& yvel_lo = lbound(yvel_arr);
    const auto& yvel_hi = ubound(yvel_arr);
 
    // xlo: ori = 0
    // ylo: ori = 1
    // zlo: ori = 2
    // xhi: ori = 3
    // yhi: ori = 4
    // zhi: ori = 5
 
    // Based on BCRec for the domain, we need to make BCRec for this Box
    // bccomp is used as starting index for m_domain_bcs_type
    //      0 is used as starting index for bcrs
    Vector<BCRec> bcrs(1);
    setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs);
 
    Gpu::DeviceVector<BCRec> bcrs_d(1);
    Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
    const BCRec* bc_ptr = bcrs_d.data();
 
    GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,1> l_bc_extdir_vals_d;
 
    for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
        l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori];
    }
 
    GeometryData const& geomdata = m_geom.data();
    bool is_periodic_in_x = geomdata.isPeriodic(0);
    bool is_periodic_in_y = geomdata.isPeriodic(1);
 
    // First do all ext_dir bcs
    if (!is_periodic_in_x)
    {
        Real* xvel_bc_ptr = m_u_bc_data;
        Box bx_xlo(bx);  bx_xlo.setBig  (0,dom_lo.x-1);
        Box bx_xhi(bx);  bx_xhi.setSmall(0,dom_hi.x+2);
        Box bx_xlo_face(bx); bx_xlo_face.setSmall(0,dom_lo.x  ); bx_xlo_face.setBig(0,dom_lo.x  );
        Box bx_xhi_face(bx); bx_xhi_face.setSmall(0,dom_hi.x+1); bx_xhi_face.setBig(0,dom_hi.x+1);
        ParallelFor(bx_xlo, bx_xlo_face,
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                int iflip = dom_lo.x - i;
                if (bc_ptr[0].lo(0) == ERFBCType::ext_dir) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0];
                } else if (bc_ptr[0].lo(0) == ERFBCType::ext_dir_upwind && xvel_arr(dom_lo.x,Clamp(j,xvel_lo.y,xvel_hi.y),Clamp(k,xvel_lo.z,xvel_hi.z)) >= zero) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0];
                } else if (bc_ptr[0].lo(0) == ERFBCType::foextrap) {
                    dest_arr(i,j,k) =  dest_arr(dom_lo.x,j,k);
                } else if (bc_ptr[0].lo(0) == ERFBCType::open) {
                    dest_arr(i,j,k) =  dest_arr(dom_lo.x,j,k);
                } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_even) {
                    dest_arr(i,j,k) =  dest_arr(iflip,j,k);
                } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_odd) {
                    dest_arr(i,j,k) = -dest_arr(iflip,j,k);
                } else if (bc_ptr[0].lo(0) == ERFBCType::neumann_int) {
                    dest_arr(i,j,k) = (Real(4.0)*dest_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+2,j,k))/three;
                } else if (bc_ptr[0].lo(0) == ERFBCType::hoextrap) {
                    Real delta_i = (dom_lo.x - i);
                    dest_arr(i,j,k) = (one + delta_i)*dest_arr(dom_lo.x,j,k) - delta_i*dest_arr(dom_lo.x+1,j,k);
                }
            },
            // We only set the values on the domain faces themselves if EXT_DIR or neumann_int
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                if (bc_ptr[0].lo(0) == ERFBCType::ext_dir) {
                      dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0];
                } else if (bc_ptr[0].lo(0) == ERFBCType::ext_dir_upwind && xvel_arr(dom_lo.x,Clamp(j,xvel_lo.y,xvel_hi.y),Clamp(k,xvel_lo.z,xvel_hi.z)) >= zero) {
                      dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0];
                } else if (bc_ptr[0].lo(0) == ERFBCType::neumann_int) {
                      dest_arr(i,j,k) = (Real(4.0)*dest_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+2,j,k))/three;
                }
            }
        );
        ParallelFor(bx_xhi, bx_xhi_face,
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                int iflip =  2*(dom_hi.x + 1) - i;
                if (bc_ptr[0].hi(0) == ERFBCType::ext_dir) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3];
                } else if (bc_ptr[0].hi(0) == ERFBCType::ext_dir_upwind && xvel_arr(dom_hi.x+1,Clamp(j,xvel_lo.y,xvel_hi.y),Clamp(k,xvel_lo.z,xvel_hi.z)) <= zero) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3];
                } else if (bc_ptr[0].hi(0) == ERFBCType::foextrap) {
                    dest_arr(i,j,k) =  dest_arr(dom_hi.x+1,j,k);
                } else if (bc_ptr[0].hi(0) == ERFBCType::open) {
                    dest_arr(i,j,k) =  dest_arr(dom_hi.x+1,j,k);
                } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_even) {
                    dest_arr(i,j,k) =  dest_arr(iflip,j,k);
                } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_odd) {
                    dest_arr(i,j,k) = -dest_arr(iflip,j,k);
                } else if (bc_ptr[0].hi(0) == ERFBCType::neumann_int) {
                    dest_arr(i,j,k) = (Real(4.0)*dest_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x-1,j,k))/three;
                } else if (bc_ptr[0].hi(0) == ERFBCType::hoextrap) {
                    Real delta_i = (i - dom_hi.x - 1);
                    dest_arr(i,j,k) = (one + delta_i)*dest_arr(dom_hi.x+1,j,k) - delta_i*dest_arr(dom_hi.x,j,k);
                }
            },
            // We only set the values on the domain faces themselves if EXT_DIR or neumann_int
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                if (bc_ptr[0].hi(0) == ERFBCType::ext_dir) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3];
                } else if (bc_ptr[0].hi(0) == ERFBCType::ext_dir_upwind && xvel_arr(dom_hi.x+1,Clamp(j,xvel_lo.y,xvel_hi.y),Clamp(k,xvel_lo.z,xvel_hi.z)) <= zero) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3];
                } else if (bc_ptr[0].hi(0) == ERFBCType::neumann_int) {
                    dest_arr(i,j,k) = (Real(4.0)*dest_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x-1,j,k))/three;
                }
            }
        );
    } // not periodic in x
 
    if (!is_periodic_in_y)
    {
        // Populate ghost cells on lo-y and hi-y domain boundaries
        Real* xvel_bc_ptr = m_u_bc_data;
        Box bx_ylo(bx);  bx_ylo.setBig  (1,dom_lo.y-1);
        Box bx_yhi(bx);  bx_yhi.setSmall(1,dom_hi.y+1);
        ParallelFor(bx_ylo, bx_yhi,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                int jflip = dom_lo.y - 1 - j;
                if (bc_ptr[0].lo(1) == ERFBCType::ext_dir) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1];
                } else if (bc_ptr[0].lo(1) == ERFBCType::ext_dir_upwind && yvel_arr(Clamp(i,yvel_lo.x,yvel_hi.x),dom_lo.y,Clamp(k,yvel_lo.z,yvel_hi.z)) >= zero) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1];
                } else if (bc_ptr[0].lo(1) == ERFBCType::foextrap) {
                    dest_arr(i,j,k) =  dest_arr(i,dom_lo.y,k);
                } else if (bc_ptr[0].lo(1) == ERFBCType::open) {
                    dest_arr(i,j,k) =  dest_arr(i,dom_lo.y,k);
                } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_even) {
                    dest_arr(i,j,k) =  dest_arr(i,jflip,k);
                } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_odd) {
                    dest_arr(i,j,k) = -dest_arr(i,jflip,k);
                } else if (bc_ptr[0].lo(1) == ERFBCType::hoextrap) {
                    Real delta_j = (dom_lo.y - j);
                    dest_arr(i,j,k) = (one + delta_j)*dest_arr(i,dom_lo.y,k) - delta_j*dest_arr(i,dom_lo.y+1,k);
                }
            },
            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                int jflip =  2*dom_hi.y + 1 - j;
                if (bc_ptr[0].hi(1) == ERFBCType::ext_dir) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4];
                } else if (bc_ptr[0].hi(1) == ERFBCType::ext_dir_upwind && yvel_arr(Clamp(i,yvel_lo.x,yvel_hi.x),dom_hi.y+1,Clamp(k,yvel_lo.z,yvel_hi.z)) <= zero) {
                    dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4];
                } else if (bc_ptr[0].hi(1) == ERFBCType::foextrap) {
                    dest_arr(i,j,k) =  dest_arr(i,dom_hi.y,k);
                } else if (bc_ptr[0].hi(1) == ERFBCType::open) {
                    dest_arr(i,j,k) =  dest_arr(i,dom_hi.y,k);
                } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_even) {
                    dest_arr(i,j,k) =  dest_arr(i,jflip,k);
                } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_odd) {
                    dest_arr(i,j,k) = -dest_arr(i,jflip,k);
                } else if (bc_ptr[0].hi(1) == ERFBCType::hoextrap) {
                    Real delta_j = (j - dom_hi.y);
                    dest_arr(i,j,k) = (one + delta_j)*dest_arr(i,dom_hi.y,k) - delta_j*dest_arr(i,dom_hi.y-1,k);
                }
            }
        );
    } // not periodic in y
 
    Gpu::streamSynchronize();
}

Here is the call graph for this function:

impose_vertical_xvel_bcs()

void ERFPhysBCFunct_u::impose_vertical_xvel_bcs	(	const amrex::Array4< amrex::Real > &	dest_arr,
		const amrex::Box &	bx,
		const amrex::Box &	domain,
		const amrex::Array4< amrex::Real const > &	z_nd,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM >	dxInv,
		int	bccomp,
		const double	time
	)

Impose vertical boundary conditions on x velocity.

Parameters

[in,out]	dest_arr	x-velocity data to be filled
[in]	bx	box holding data to be filled
[in]	domain	simulation domain
[in]	z_nd	nodal physical height
[in]	dxInv	inverse cell size
[in]	bccomp	boundary-condition component
[in]	time	time at which data should be filled

Impose vertical boundary conditions on x-component of velocity

Parameters

[in,out]	dest_arr	the Array4 of the quantity to be filled
[in]	bx	the box associated with this data
[in]	domain	the computational domain
[in]	z_phys_nd	height coordinate at nodes
[in]	dxInv	inverse cell size array
[in]	bccomp	index into m_domain_bcs_type
[in]	time	time at which the data should be filled

{
    BL_PROFILE_VAR("impose_vertical_xvel_bcs()",impose_vertical_xvel_bcs);
    const auto& dom_lo = lbound(domain);
    const auto& dom_hi = ubound(domain);
 
    Box per_grown_domain(domain);
    int growx = (m_geom.isPeriodic(0)) ? 1 : 0;
    int growy = (m_geom.isPeriodic(1)) ? 1 : 0;
    per_grown_domain.grow(IntVect(growx,growy,0));
    const auto& perdom_lo = lbound(per_grown_domain);
    const auto& perdom_hi = ubound(per_grown_domain);
 
    GeometryData const& geomdata = m_geom.data();
 
    // Based on BCRec for the domain, we need to make BCRec for this Box
    // bccomp is used as starting index for m_domain_bcs_type
    //      0 is used as starting index for bcrs
    Vector<BCRec> bcrs(1);
    setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs);
 
    Gpu::DeviceVector<BCRec> bcrs_d(1);
    Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin());
    const BCRec* bc_ptr = bcrs_d.data();
 
    GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,1> l_bc_extdir_vals_d;
 
    for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
        l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori];
    }
 
    {
        // Populate ghost cells on lo-z and hi-z domain boundaries
        Box bx_zlo(bx);  bx_zlo.setBig  (2,dom_lo.z-1);
        Box bx_zhi(bx);  bx_zhi.setSmall(2,dom_hi.z+1);
        ParallelFor(bx_zlo, bx_zhi,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                int kflip = dom_lo.z - 1 - k;
                if (bc_ptr[0].lo(2) == ERFBCType::ext_dir) {
                    #ifdef ERF_USE_TERRAIN_VELOCITY
                        dest_arr(i,j,k) = prob->compute_terrain_velocity(time);
                    #else
                        dest_arr(i,j,k) = l_bc_extdir_vals_d[0][2];
                    #endif
                } else if (bc_ptr[0].lo(2) == ERFBCType::foextrap) {
                    dest_arr(i,j,k) =  dest_arr(i,j,dom_lo.z);
                } else if (bc_ptr[0].lo(2) == ERFBCType::open) {
                    dest_arr(i,j,k) =  dest_arr(i,j,dom_lo.z);
                } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_even) {
                    dest_arr(i,j,k) =  dest_arr(i,j,kflip);
                } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_odd) {
                    dest_arr(i,j,k) = -dest_arr(i,j,kflip);
                } else if (bc_ptr[0].lo(2) == ERFBCType::hoextrap) {
                    Real delta_k = (dom_lo.z - k);
                    dest_arr(i,j,k) = (one + delta_k)*dest_arr(i,j,dom_lo.z) - delta_k*dest_arr(i,j,dom_lo.z+1);
                }
            },
            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                int kflip =  2*dom_hi.z + 1 - k;
                if (bc_ptr[0].hi(2) == ERFBCType::ext_dir) {
                    dest_arr(i,j,k) = l_bc_extdir_vals_d[0][5];
                } else if (bc_ptr[0].hi(2) == ERFBCType::foextrap) {
                    dest_arr(i,j,k) =  dest_arr(i,j,dom_hi.z);
                } else if (bc_ptr[0].hi(2) == ERFBCType::open) {
                    dest_arr(i,j,k) =  dest_arr(i,j,dom_hi.z);
                } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_even) {
                    dest_arr(i,j,k) =  dest_arr(i,j,kflip);
                } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_odd) {
                    dest_arr(i,j,k) = -dest_arr(i,j,kflip);
                } else if (bc_ptr[0].hi(2) == ERFBCType::hoextrap){
                    Real delta_k = (k - dom_hi.z);
                    dest_arr(i,j,k) = (one + delta_k)*dest_arr(i,j,dom_hi.z) - delta_k*dest_arr(i,j,dom_hi.z-1);
                }
            }
        );
    } // z
 
    if (m_z_phys_nd) {
        const auto&  bx_lo = lbound(bx);
        const auto&  bx_hi = ubound(bx);
 
        const auto& zphys_lo = lbound(Box(z_phys_nd));
        const auto& zphys_hi = ubound(Box(z_phys_nd));
 
        // Neumann conditions (d<var>/dn = 0) must be aware of the surface normal with terrain.
        // An additional source term arises from d<var>/dx & d<var>/dy & met_h_xi/eta/zeta.
        //=====================================================================================
        // Only modify scalars, U, or V
        // Loop over each component
        // Hit for Neumann condition at kmin
        if (bcrs[0].lo(2) == ERFBCType::foextrap) {
            // Loop over ghost cells in bottom XY-plane (valid box)
            Box xybx = bx;
            xybx.setBig(2,-1);
            xybx.setSmall(2,bx.smallEnd()[2]);
            int k0 = 0;
 
            // Get the dz cell size
            Real dz = geomdata.CellSize(2);
 
            // Fill all the Neumann srcs with terrain
            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                // Clip indices for ghost-cells
                int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x);
                    ii = amrex::min(amrex::max(ii,zphys_lo.x),zphys_hi.x);
                int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y);
                    jj = amrex::min(amrex::max(jj,zphys_lo.y),zphys_hi.y);
 
                // Get metrics
                Real met_h_xi   = Compute_h_xi_AtIface  (ii, jj, k0, dxInv, z_phys_nd);
                Real met_h_eta  = Compute_h_eta_AtIface (ii, jj, k0, dxInv, z_phys_nd);
                Real met_h_zeta = Compute_h_zeta_AtIface(ii, jj, k0, dxInv, z_phys_nd);
 
                // GradX at IJK location inside domain -- this relies on the assumption that we have
                // used foextrap for quantities outside the domain to define the gradient as zero
                Real GradVarx, GradVary;
                if ( i < dom_lo.x-1 || i > dom_hi.x+1 || (i+1 > bx_hi.x && i-1 < bx_lo.x) ) {
                    GradVarx = zero;
                } else if (i+1 > bx_hi.x) {
                    GradVarx =       dxInv[0] * (dest_arr(i  ,j,k0) - dest_arr(i-1,j,k0));
                } else if (i-1 < bx_lo.x) {
                    GradVarx =       dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i  ,j,k0));
                } else {
                    GradVarx = myhalf * dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i-1,j,k0));
                }
 
                // GradY at IJK location inside domain -- this relies on the assumption that we have
                // used foextrap for quantities outside the domain to define the gradient as zero
                if ( j < dom_lo.y-1 || j > dom_hi.y+1 || (j+1 > bx_hi.y && j-1 < bx_lo.y) ) {
                    GradVary = zero;
                } else if (j+1 > bx_hi.y) {
                    GradVary =       dxInv[1] * (dest_arr(i,j  ,k0) - dest_arr(i,j-1,k0));
                } else if (j-1 < bx_lo.y) {
                    GradVary =       dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j  ,k0));
                } else {
                    GradVary = myhalf * dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j-1,k0));
                }
 
                // Prefactor
                Real met_fac =  met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + one );
 
                // Accumulate in bottom ghost cell (EXTRAP already populated)
                dest_arr(i,j,k) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary );
            });
        } // foextrap
    } //m_z_phys_nd
    Gpu::streamSynchronize();
}

Here is the call graph for this function:

operator()()

void ERFPhysBCFunct_u::operator()	(	amrex::MultiFab &	mf,
		amrex::MultiFab &	xvel,
		amrex::MultiFab &	yvel,
		amrex::IntVect const &	nghost,
		const double	time,
		int	bccomp,
		bool	do_fb
	)

Impose physical boundary conditions at domain boundaries

Parameters

[in,out]	mf	x-velocity MultiFab to be filled
[in]	xvel	x-velocity MultiFab used by upwind boundary conditions
[in]	yvel	y-velocity MultiFab used by upwind boundary conditions
[in]	nghost	number of ghost cells to be filled for velocity components
[in]	time	time at which the data should be filled
[in]	bccomp	index into m_domain_bcs_type
[in]	do_fb	whether to fill periodic/interior ghost cells first

Impose physical boundary conditions on x velocity.

Parameters

[in,out]	mf	x-velocity MultiFab to be filled
[in]	xvel	x-velocity MultiFab used by upwind boundary conditions
[in]	yvel	y-velocity MultiFab used by upwind boundary conditions
[in]	nghost	number of ghost cells to be filled
[in]	time	time at which data should be filled
[in]	bccomp	index into m_domain_bcs_type
[in]	do_fb	whether to fill periodic/interior ghost cells first

{
    BL_PROFILE("ERFPhysBCFunct_u::()");
 
    //
    // We fill all of the interior and periodic ghost cells first, so we can fill
    //    those directly inside the lateral and vertical calls.
    // If triply periodic this is all we do
    //
    if (do_fb) {
        mf.FillBoundary(m_geom.periodicity());
    }
 
    if (m_geom.isAllPeriodic()) return;
 
    const auto& domain = m_geom.Domain();
    const auto dxInv   = m_geom.InvCellSizeArray();
 
    Box gdomainx = surroundingNodes(domain,0);
    for (int i = 0; i < AMREX_SPACEDIM; ++i) {
        if (m_geom.isPeriodic(i)) {
            gdomainx.grow(i, nghost[i]);
        }
    }
 
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    {
        for (MFIter mfi(mf,false); mfi.isValid(); ++mfi)
        {
            //
            // This is the box we pass to the different routines
            // NOTE -- this is the full grid NOT the tile box
            //
            Box bx  = mfi.validbox();
 
            //
            // These are the boxes we use to test on relative to the domain
            //
            Box xbx1 = surroundingNodes(bx,0); xbx1.grow(nghost);
            if(xbx1.smallEnd(2) < domain.smallEnd(2)) xbx1.setSmall(2,domain.smallEnd(2));
            if(xbx1.bigEnd(2)   > domain.bigEnd(2))   xbx1.setBig(2,domain.bigEnd(2));
 
            Box xbx2 = surroundingNodes(bx,0); xbx2.grow(nghost);
 
            Array4<const Real> z_nd_arr;
 
            if (m_z_phys_nd) {
                z_nd_arr = m_z_phys_nd->const_array(mfi);
            }
 
            if (!gdomainx.contains(xbx2)) {
                Array4<      Real> const& dest_arr = mf.array(mfi);
                Array4<const Real> const& velx_arr = xvel.const_array(mfi);
                Array4<const Real> const& vely_arr = yvel.const_array(mfi);
 
                if (!gdomainx.contains(xbx1) && !m_use_real_bcs) {
                    impose_lateral_xvel_bcs(dest_arr,velx_arr,vely_arr,xbx1,domain,bccomp,time);
                }
 
                impose_vertical_xvel_bcs(dest_arr,xbx2,domain,z_nd_arr,dxInv,bccomp,time);
            }
        } // MFIter
    } // OpenMP
 
    // The terrain correction below the ground is computed box by box; make the
    // copies of each ghost cell agree.  Only a terrain-fitted mesh has a
    // non-zero correction -- z_phys_nd is allocated (and zeroed) for every
    // terrain type, so testing it alone would pay for the sync on flat runs.
    if (m_z_phys_nd && terrain_is_fitted(m_terrain_type)) {
        m_below_ground_sync(mf, 0, 1, nghost, m_geom);
    }
} // operator()

Member Data Documentation

m_bc_extdir_vals

amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NBCVAR_max> ERFPhysBCFunct_u::m_bc_extdir_vals

private

m_bc_neumann_vals

amrex::Array<amrex::Array<amrex::Real, AMREX_SPACEDIM*2>,AMREX_SPACEDIM+NBCVAR_max> ERFPhysBCFunct_u::m_bc_neumann_vals

private

m_below_ground_sync

BelowGroundGhostSync ERFPhysBCFunct_u::m_below_ground_sync

private

m_domain_bcs_type

amrex::Vector<amrex::BCRec> ERFPhysBCFunct_u::m_domain_bcs_type

private

m_domain_bcs_type_d

amrex::Gpu::DeviceVector<amrex::BCRec> ERFPhysBCFunct_u::m_domain_bcs_type_d

private

m_geom

amrex::Geometry ERFPhysBCFunct_u::m_geom

private

m_lev

int ERFPhysBCFunct_u::m_lev

private

m_terrain_type

TerrainType ERFPhysBCFunct_u::m_terrain_type

private

m_u_bc_data

amrex::Real* ERFPhysBCFunct_u::m_u_bc_data

private

m_use_real_bcs

bool ERFPhysBCFunct_u::m_use_real_bcs

private

m_z_phys_nd

amrex::MultiFab* ERFPhysBCFunct_u::m_z_phys_nd

private

---

The documentation for this class was generated from the following files:

• Source/BoundaryConditions/ERF_PhysBCFunct.H
• Source/BoundaryConditions/ERF_BoundaryConditionsXvel.cpp
• Source/BoundaryConditions/ERF_PhysBCFunct.cpp