ERF_AdvectionSrcForMom.cpp File Reference

#include "AMReX_BCRec.H"
#include "ERF_Advection.H"
#include "ERF_AdvectionSrcForMom_N.H"
#include "ERF_AdvectionSrcForMom_T.H"
#include "ERF_EB.H"

Include dependency graph for ERF_AdvectionSrcForMom.cpp:

Functions
void	AdvectionSrcForMom (const MFIter &mfi, const Box &bxx, const Box &bxy, const Box &bxz, const Vector< Box > &bxx_grown, const Vector< Box > &bxy_grown, const Vector< Box > &bxz_grown, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &cell_data, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &omega, const Array4< const Real > &z_nd, const Array4< const Real > &ax, const Array4< const Real > &ay, const Array4< const Real > &az, const Array4< const Real > &detJ, Gpu::DeviceVector< Real > &stretched_dz_d, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_ux, const Array4< const Real > &mf_vx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vy, const AdvType horiz_adv_type, const AdvType vert_adv_type, const Real horiz_upw_frac, const Real vert_upw_frac, MeshType &mesh_type, TerrainType &terrain_type, const eb_ &ebfact, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_u_arr, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_v_arr, GpuArray< Array4< Real >, AMREX_SPACEDIM > &flx_w_arr, const Vector< iMultiFab > &physbnd_mask, const bool already_on_centroids, const int lo_z_face, const int hi_z_face, const Box &domain, const BCRec *bc_ptr_h)

Function Documentation

AdvectionSrcForMom()

void AdvectionSrcForMom	(	const MFIter &	mfi,
		const Box &	bxx,
		const Box &	bxy,
		const Box &	bxz,
		const Vector< Box > &	bxx_grown,
		const Vector< Box > &	bxy_grown,
		const Vector< Box > &	bxz_grown,
		const Array4< Real > &	rho_u_rhs,
		const Array4< Real > &	rho_v_rhs,
		const Array4< Real > &	rho_w_rhs,
		const Array4< const Real > &	cell_data,
		const Array4< const Real > &	u,
		const Array4< const Real > &	v,
		const Array4< const Real > &	w,
		const Array4< const Real > &	rho_u,
		const Array4< const Real > &	rho_v,
		const Array4< const Real > &	omega,
		const Array4< const Real > &	z_nd,
		const Array4< const Real > &	ax,
		const Array4< const Real > &	ay,
		const Array4< const Real > &	az,
		const Array4< const Real > &	detJ,
		Gpu::DeviceVector< Real > &	stretched_dz_d,
		const GpuArray< Real, AMREX_SPACEDIM > &	cellSizeInv,
		const Array4< const Real > &	mf_mx,
		const Array4< const Real > &	mf_ux,
		const Array4< const Real > &	mf_vx,
		const Array4< const Real > &	mf_my,
		const Array4< const Real > &	mf_uy,
		const Array4< const Real > &	mf_vy,
		const AdvType	horiz_adv_type,
		const AdvType	vert_adv_type,
		const Real	horiz_upw_frac,
		const Real	vert_upw_frac,
		MeshType &	mesh_type,
		TerrainType &	terrain_type,
		const eb_ &	ebfact,
		GpuArray< Array4< Real >, AMREX_SPACEDIM > &	flx_u_arr,
		GpuArray< Array4< Real >, AMREX_SPACEDIM > &	flx_v_arr,
		GpuArray< Array4< Real >, AMREX_SPACEDIM > &	flx_w_arr,
		const Vector< iMultiFab > &	physbnd_mask,
		const bool	already_on_centroids,
		const int	lo_z_face,
		const int	hi_z_face,
		const Box &	domain,
		const BCRec *	bc_ptr_h
	)

Function for computing the advective tendency for the momentum equations This routine has explicit expressions for all cases (terrain or not) when the horizontal and vertical spatial orders are <= 2, and calls more specialized functions when either (or both) spatial order(s) is greater than two

Parameters

[in]	mfi	MultiFab Iterator
[in]	bxx	box over which the x-momentum is updated
[in]	bxy	box over which the y-momentum is updated
[in]	bxz	box over which the z-momentum is updated
[in]	bxx_grown	grown boxes of bxx to loop over the nodal grids of bxx
[in]	bxy_grown	grown boxes of bxy to loop over the nodal grids of bxy
[in]	bxz_grown	grown boxes of bxz to loop over the nodal grids of bxz
[out]	rho_u_rhs	tendency for the x-momentum equation
[out]	rho_v_rhs	tendency for the y-momentum equation
[out]	rho_w_rhs	tendency for the z-momentum equation
[in]	cell_data	cell-centered conservative state used by open-boundary terms
[in]	u	x-component of the velocity
[in]	v	y-component of the velocity
[in]	w	z-component of the velocity
[in]	rho_u	x-component of the momentum
[in]	rho_v	y-component of the momentum
[in]	omega	component of the momentum normal to the z-coordinate surface
[in]	z_nd	height coordinate at nodes
[in]	ax	Area fraction of x-faces
[in]	ay	Area fraction of y-faces
[in]	az	Area fraction of z-faces
[in]	detJ	Jacobian of the metric transformation (= 1 if use_terrain_fitted_coords is false)
[in]	stretched_dz_d	device vector of vertical cell spacing
[in]	cellSizeInv	inverse of the grid spacing
[in]	mf_mx	x map factor at cell centers
[in]	mf_ux	x map factor at x-faces
[in]	mf_vx	x map factor at y-faces
[in]	mf_my	y map factor at cell centers
[in]	mf_uy	y map factor at x-faces
[in]	mf_vy	y map factor at y-faces
[in]	horiz_adv_type	sets the spatial order to be used for lateral derivatives
[in]	vert_adv_type	sets the spatial order to be used for vertical derivatives
[in]	horiz_upw_frac	horizontal upwind blending fraction
[in]	vert_upw_frac	vertical upwind blending fraction
[in]	mesh_type	grid spacing type
[in]	terrain_type	terrain representation type
[in]	ebfact	EB factories for cell- and face-centered variables
[in,out]	flx_u_arr	container of fluxes for x-momentum
[in,out]	flx_v_arr	container of fluxes for y-momentum
[in,out]	flx_w_arr	container of fluxes for z-momentum
[in]	physbnd_mask	Vector of masks for flux interpolation (=1 otherwise, =0 if physbnd)
[in]	already_on_centroids	flag whether flux interpolation is unnecessary
[in]	lo_z_face	minimum z-face k-index at this level
[in]	hi_z_face	maximum z-face k-index at this level
[in]	domain	computational domain
[in]	bc_ptr_h	boundary condition records on host

{
    BL_PROFILE_VAR("AdvectionSrcForMom", AdvectionSrcForMom);
 
    AMREX_ALWAYS_ASSERT(bxz.smallEnd(2) > 0);
 
    // compute mapfactor inverses
    Box box2d_u(bxx);   box2d_u.setRange(2,0);   box2d_u.grow({3,3,0});
    Box box2d_v(bxy);   box2d_v.setRange(2,0);   box2d_v.grow({3,3,0});
    FArrayBox mf_u_invFAB(box2d_u,1,The_Async_Arena());
    FArrayBox mf_v_invFAB(box2d_v,1,The_Async_Arena());
    const Array4<Real>& mf_u_inv = mf_u_invFAB.array();
    const Array4<Real>& mf_v_inv = mf_v_invFAB.array();
 
    const bool use_terrain_fitted_coords = ( terrain_type == TerrainType::StaticFittedMesh ||
                                             terrain_type == TerrainType::MovingFittedMesh);
 
    ParallelFor(box2d_u, box2d_v,
    [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
    {
        mf_u_inv(i,j,0) = one / mf_ux(i,j,0);
    },
    [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
    {
        mf_v_inv(i,j,0) = one / mf_vy(i,j,0);
    });
 
    if (mesh_type == MeshType::ConstantDz && terrain_type != TerrainType::EB)
    {
        // amrex::Print() << "ADV:CONSTANT DZ " << std::endl;
        AdvectionSrcForMom_ConstantDz(bxx, bxy, bxz,
                                      rho_u_rhs, rho_v_rhs, rho_w_rhs, u, v, w,
                                      rho_u, rho_v, omega,
                                      cellSizeInv, stretched_dz_d,
                                      mf_mx, mf_ux, mf_vx, mf_my, mf_uy, mf_vy,
                                      horiz_adv_type, vert_adv_type,
                                      horiz_upw_frac, vert_upw_frac,
                                      terrain_type, lo_z_face, hi_z_face);
    }
    else if (mesh_type == MeshType::StretchedDz && terrain_type != TerrainType::EB)
    {
        // amrex::Print() << "ADV:STRETCHED DZ " << std::endl;
        AdvectionSrcForMom_StretchedDz(bxx, bxy, bxz,
                                       rho_u_rhs, rho_v_rhs, rho_w_rhs,
                                       u, v, w, rho_u, rho_v, omega,
                                       cellSizeInv, stretched_dz_d,
                                       mf_mx, mf_ux, mf_vx, mf_my, mf_uy, mf_vy,
                                       horiz_adv_type, vert_adv_type,
                                       horiz_upw_frac, vert_upw_frac,
                                       lo_z_face, hi_z_face);
    }
    else if ( terrain_type == TerrainType::EB)
    {
        // amrex::Print() << "ADV:EB " << std::endl;
        AdvectionSrcForMom_EB(mfi, bxx, bxy, bxz, bxx_grown, bxy_grown, bxz_grown,
                              rho_u_rhs, rho_v_rhs, rho_w_rhs,
                              u, v, w,
                              rho_u, rho_v, omega,
                              cellSizeInv,
                              mf_mx, mf_ux, mf_vx, mf_my, mf_uy, mf_vy,
                              horiz_adv_type, vert_adv_type,
                              horiz_upw_frac, vert_upw_frac,
                              ebfact, flx_u_arr, flx_v_arr, flx_w_arr,
                              physbnd_mask, already_on_centroids,
                              lo_z_face, hi_z_face, domain);
    }
    else
    {
        AMREX_ALWAYS_ASSERT(use_terrain_fitted_coords);
        // amrex::Print() << "ADV:TF " << std::endl;
        AdvectionSrcForMom_TF(bxx, bxy, bxz,
                              rho_u_rhs, rho_v_rhs, rho_w_rhs,
                              u, v, w,
                              rho_u, rho_v, omega,
                              z_nd, ax, ay, az, detJ,
                              cellSizeInv,
                              mf_mx, mf_ux, mf_vx, mf_my, mf_uy, mf_vy,
                              horiz_adv_type, vert_adv_type,
                              horiz_upw_frac, vert_upw_frac,
                              lo_z_face, hi_z_face);
 
    }
 
    // Open bc will be imposed upon all vars (we only access cons here for simplicity)
    const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open);
    const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open);
    const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open);
    const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open);
 
    // The open-boundary kernels below overwrite the tendencies the operators above just
    //    wrote, so they have to run on exactly the boxes those operators covered.  Rebuilding
    //    the boxes here instead is what went wrong: the z shrink was applied unconditionally,
    //    while the caller applies it only where the box actually reaches the bottom or top of
    //    the domain.  On a grid decomposed in z that dropped one w-face at every internal
    //    split from the open-boundary treatment, leaving in place the interior stencil --
    //    which differences across the open lateral boundary, into an unfilled ghost cell.
    const Box& tbx = bxx;
    const Box& tby = bxy;
    const Box& tbz = bxz;
 
    const int domhi_z = domain.bigEnd(2);
 
    // Special advection operator for open BC (bndry normal/tangent operations)
    //
    // Where two perpendicular open boundaries meet, every face and cell on the
    //    boundary must be updated by exactly one of these kernels: they assign rather
    //    than accumulate, so an overlap would silently keep whichever ran last and
    //    throw the other away.  A boundary-normal momentum owns the corner in its own
    //    direction -- that face lies on the open boundary, so the radiation condition
    //    is what belongs there -- and the tangential boxes are trimmed to match.
    //    w is tangent to both lateral boundaries and so has no normal owner at a
    //    corner; it is covered by OpenBCTangentPatches below, which tags the corner
    //    with both sides so that the kernel differences across neither open boundary.
    if (xlo_open)
    {
        Box tbx_xlo, tby_xlo;
        if (tbx.smallEnd(0) == domain.smallEnd(0)) { tbx_xlo = makeSlab(tbx,0,domain.smallEnd(0));}
        if (tby.smallEnd(0) == domain.smallEnd(0)) { tby_xlo = makeSlab(tby,0,domain.smallEnd(0));}
 
        tby_xlo = TrimOpenBCCorner(tby_xlo, 1, domain, ylo_open, yhi_open);
 
        bool do_lo = true;
 
        AdvectionSrcForOpenBC_Normal(tbx_xlo, 0, rho_u_rhs, u, cell_data, cellSizeInv, do_lo);
        AdvectionSrcForOpenBC_Tangent_Ymom(tby_xlo, 0, rho_v_rhs, v,
                                           rho_u, rho_v, omega,
                                           ay, az, detJ, cellSizeInv,
                                           do_lo);
    }
    if (xhi_open)
    {
        Box tbx_xhi, tby_xhi;
        if (tbx.bigEnd(0) == domain.bigEnd(0)+1)   { tbx_xhi = makeSlab(tbx,0,domain.bigEnd(0)+1);}
        if (tby.bigEnd(0) == domain.bigEnd(0))     { tby_xhi = makeSlab(tby,0,domain.bigEnd(0)  );}
 
        tby_xhi = TrimOpenBCCorner(tby_xhi, 1, domain, ylo_open, yhi_open);
 
        AdvectionSrcForOpenBC_Normal(tbx_xhi, 0, rho_u_rhs, u, cell_data, cellSizeInv);
        AdvectionSrcForOpenBC_Tangent_Ymom(tby_xhi, 0, rho_v_rhs, v,
                                           rho_u, rho_v, omega,
                                           ay, az, detJ, cellSizeInv);
    }
    if (ylo_open)
    {
        Box tbx_ylo, tby_ylo;
        if (tbx.smallEnd(1) == domain.smallEnd(1)) { tbx_ylo = makeSlab(tbx,1,domain.smallEnd(1));}
        if (tby.smallEnd(1) == domain.smallEnd(1)) { tby_ylo = makeSlab(tby,1,domain.smallEnd(1));}
 
        tbx_ylo = TrimOpenBCCorner(tbx_ylo, 0, domain, xlo_open, xhi_open);
 
        bool do_lo = true;
        AdvectionSrcForOpenBC_Tangent_Xmom(tbx_ylo, 1, rho_u_rhs, u,
                                           rho_u, rho_v, omega,
                                           ax, az, detJ, cellSizeInv,
                                           do_lo);
        AdvectionSrcForOpenBC_Normal(tby_ylo, 1, rho_v_rhs, v, cell_data, cellSizeInv, do_lo);
    }
    if (yhi_open)
    {
        Box tbx_yhi, tby_yhi;
        if (tbx.bigEnd(1) == domain.bigEnd(1))     { tbx_yhi = makeSlab(tbx,1,domain.bigEnd(1)  );}
        if (tby.bigEnd(1) == domain.bigEnd(1)+1)   { tby_yhi = makeSlab(tby,1,domain.bigEnd(1)+1);}
 
        tbx_yhi = TrimOpenBCCorner(tbx_yhi, 0, domain, xlo_open, xhi_open);
 
        AdvectionSrcForOpenBC_Tangent_Xmom(tbx_yhi, 1, rho_u_rhs, u,
                                           rho_u, rho_v, omega,
                                           ax, az, detJ, cellSizeInv);
        AdvectionSrcForOpenBC_Normal(tby_yhi, 1, rho_v_rhs, v, cell_data, cellSizeInv);
    }
 
    for (const OpenBCPatch& patch : OpenBCTangentPatches(tbz, domain,
                                                         xlo_open, xhi_open, ylo_open, yhi_open))
    {
        AdvectionSrcForOpenBC_Tangent_Zmom(patch.box, patch.x_side, patch.y_side, rho_w_rhs, w,
                                           rho_u, rho_v, omega,
                                           ax, ay, az, detJ, cellSizeInv,
                                           domhi_z);
    }
}

Here is the call graph for this function: