REMORA_BoundaryConditions_zvel.cpp

Go to the documentation of this file.

#include "AMReX_PhysBCFunct.H"
#include <REMORA_PhysBCFunct.H>
#include <REMORA_prob_common.H>
 
using namespace amrex;
 
/**
 * @param[inout] dest_arr      data on which to apply BCs
 * @param[in   ] bx            box to update on
 * @param[in   ] domain        domain box
 * @param[in   ] dxInv         pm or pn
 * @param[in   ] mskr          land-sea mask on rho-points
 * @param[in   ] calc_arr      data to use in the RHS of calculations
 * @param[in   ] time          current time
 * @param[in   ] bccomp        index into both domain_bcs_type_bcr and bc_extdir_vals for icomp=0
 */
void REMORAPhysBCFunct::impose_zvel_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
                                        const GpuArray<Real,AMREX_SPACEDIM> /*dxInv*/,const Array4<const Real>& mskr,
                                        Real /*time*/, int bccomp)
{
    const auto& dom_lo = amrex::lbound(domain);
    const auto& dom_hi = amrex::ubound(domain);
 
    // 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
    int ncomp = 1;
    Vector<BCRec> bcrs(ncomp);
    amrex::setBC(bx, domain, bccomp, 0, ncomp, m_domain_bcs_type, bcrs);
 
    // 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
    amrex::setBC(bx, domain, bccomp, 0, ncomp, m_domain_bcs_type, bcrs);
 
    amrex::Gpu::DeviceVector<BCRec> bcrs_d(ncomp);
#ifdef AMREX_USE_GPU
    Gpu::htod_memcpy_async(bcrs_d.data(), bcrs.data(), sizeof(BCRec)*ncomp);
#else
    std::memcpy(bcrs_d.data(), bcrs.data(), sizeof(BCRec)*ncomp);
#endif
    const amrex::BCRec* bc_ptr = bcrs_d.data();
    const auto* bc_extdir_vals_ptr = m_bc_extdir_vals_d.data();
 
    GeometryData const& geomdata = m_geom.data();
    bool is_periodic_in_x = geomdata.isPeriodic(0);
    bool is_periodic_in_y = geomdata.isPeriodic(1);
 
    // ROMS masks the boundary values of tke and gls (tkebc) but not those of w (bc_w3d), which
    // share this fill. Selected rather than multiplied, so a cell not yet filled never enters
    // arithmetic.
    const bool masked_var = (bccomp == BCVars::tke_bc(m_ncons));
 
    // First do all ext_dir bcs
    if (!is_periodic_in_x or bccomp == BCVars::foextrap_bc(m_ncons))
    {
        Box bx_xlo(bx);  bx_xlo.setBig  (0,dom_lo.x-1);
        Box bx_xhi(bx);  bx_xhi.setSmall(0,dom_hi.x+1);
        ParallelFor(
            bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
                const bool land = masked_var && mskr.contains(i,j,0) && mskr(i,j,0) == Real(0.0);
                int iflip = dom_lo.x - 1 - i;
                if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) {
                    dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][0];
                } else if (bc_ptr[n].lo(0) == REMORABCType::foextrap || bc_ptr[n].lo(0) == REMORABCType::clamped) {
                    dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(dom_lo.x,j,k);
                } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_even) {
                    dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(iflip,j,k);
                } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_odd) {
                    dest_arr(i,j,k) = land ? Real(0.0) : -dest_arr(iflip,j,k);
                }
            },
            bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
                const bool land = masked_var && mskr.contains(i,j,0) && mskr(i,j,0) == Real(0.0);
                int iflip = 2*dom_hi.x + 1 - i;
                if (bc_ptr[n].hi(0) == REMORABCType::ext_dir) {
                    dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][3];
                } else if (bc_ptr[n].hi(0) == REMORABCType::foextrap || bc_ptr[n].hi(0) == REMORABCType::clamped) {
                    dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(dom_hi.x,j,k);
                } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_even) {
                    dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(iflip,j,k);
                } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_odd) {
                    dest_arr(i,j,k) = land ? Real(0.0) : -dest_arr(iflip,j,k);
                }
            }
        );
    }
 
    // First do all ext_dir bcs
    if (!is_periodic_in_y or bccomp == BCVars::foextrap_bc(m_ncons))
    {
        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, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
            const bool land = masked_var && mskr.contains(i,j,0) && mskr(i,j,0) == Real(0.0);
            int jflip = dom_lo.y - 1 - j;
            if (bc_ptr[n].lo(1) == REMORABCType::ext_dir) {
                dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][1];
            } else if (bc_ptr[n].lo(1) == REMORABCType::foextrap || bc_ptr[n].lo(1) == REMORABCType::clamped) {
                dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(i,dom_lo.y,k);
            } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_even) {
                dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(i,jflip,k);
            } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_odd) {
                dest_arr(i,j,k) = land ? Real(0.0) : -dest_arr(i,jflip,k);
            }
        },
        bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
            const bool land = masked_var && mskr.contains(i,j,0) && mskr(i,j,0) == Real(0.0);
            int jflip =  2*dom_hi.y + 1 - j;
            if (bc_ptr[n].hi(1) == REMORABCType::ext_dir) {
                dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][4];
            } else if (bc_ptr[n].hi(1) == REMORABCType::foextrap || bc_ptr[n].hi(1) == REMORABCType::clamped) {
                dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(i,dom_hi.y,k);
            } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_even) {
                dest_arr(i,j,k) = land ? Real(0.0) : dest_arr(i,jflip,k);
            } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_odd) {
                dest_arr(i,j,k) = land ? Real(0.0) : -dest_arr(i,jflip,k);
            }
        });
    }
 
    // Never need to apply a boundary condition to the top
 
    Gpu::streamSynchronize();
}