diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp index 0d5a914c195..b871ee5305d 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp @@ -417,7 +417,7 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const const Geometry& geom = m_geom[amrlev][mglev]; const auto dxinv = geom.InvCellSizeArray(); - const Box& domain = amrex::enclosedCells(geom.Domain()); + const Box& domain = geom.Domain(); const auto dlo = amrex::lbound(domain); const auto dhi = amrex::ubound(domain); const auto& bcondloc = *m_bcondloc[amrlev][mglev]; @@ -479,21 +479,30 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const AMREX_D_TERM(Array4 const kapxfab = kapmf[0].const_array(mfi);, Array4 const kapyfab = kapmf[1].const_array(mfi);, Array4 const kapzfab = kapmf[2].const_array(mfi);); - + + const auto & bdcv = bcondloc.bndryConds(mfi); + Array2D bct; + for (int icomp = 0; icomp < AMREX_SPACEDIM; ++icomp) { + for (OrientationIter face; face; ++face) { + Orientation ori = face(); + bct(ori,icomp) = bdcv[icomp][ori]; + } + } + if (fabtyp == FabType::regular) { AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM ( xbx, txbx, { - mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv); + mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct); } , ybx, tybx, { - mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv); + mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct); } , zbx, tzbx, { - mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv); + mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct); } ); } @@ -504,16 +513,6 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const Array4 const& apz = area[2]->const_array(mfi);); Array4 const& flag = flags->const_array(mfi); - const auto & bdcv = bcondloc.bndryConds(mfi); - - Array2D bct; - for (OrientationIter face; face; ++face) { - Orientation ori = face(); - for (int icomp = 0; icomp < AMREX_SPACEDIM; ++icomp) { - bct(ori,icomp) = bdcv[icomp][ori]; - } - } - AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM ( xbx, txbx, { diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H index f2e7cf73843..419330db8bb 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H @@ -52,18 +52,33 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4 const& fx, Real dvdy = (0.5*dyi) * ((vel(i ,jhip,0,1)-vel(i ,jhim,0,1))*whi +(vel(i-1,jlop,0,1)-vel(i-1,jlom,0,1))*wlo); - if (i==dhi.x+1 and bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ - dudy = (dyi) * ((vel(i-1,jlop,0,0)-vel(i-1,jlom,0,0))*wlo); - } - if (i==dhi.x+1 and bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_NEUMANN){ - dvdy = (dyi) * ((vel(i-1,jlop,0,1)-vel(i-1,jlom,0,1))*wlo); - } - - if (i==dlo.x and bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ - dudy = (dyi) * ((vel(i ,jhip,0,0)-vel(i ,jhim,0,0))*whi); + if (i==dhi.x+1){ + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudy = (dyi) * ((vel(i-1,jlop,0,0)-vel(i-1,jlom,0,0))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdy = (dyi) * ((vel(i-1,jlop,0,1)-vel(i-1,jlom,0,1))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudy = (dyi) * ((vel(i ,jhip,0,0)-vel(i ,jhim,0,0))*whi); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdy = (dyi) * ((vel(i ,jhip,0,1)-vel(i ,jhim,0,1))*whi); + } } - if (i==dlo.x and bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_NEUMANN){ - dvdy = (dyi) * ((vel(i ,jhip,0,1)-vel(i ,jhim,0,1))*whi); + if (i==dlo.x){ + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudy = (dyi) * ((vel(i ,jhip,0,0)-vel(i ,jhim,0,0))*whi); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdy = (dyi) * ((vel(i ,jhip,0,1)-vel(i ,jhim,0,1))*whi); + } + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudy = (dyi) * ((vel(i-1,jlop,0,0)-vel(i-1,jlom,0,0))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdy = (dyi) * ((vel(i-1,jlop,0,1)-vel(i-1,jlom,0,1))*wlo); + } } Real divu = dvdy; @@ -116,21 +131,35 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4 const& fy, Real dvdx = (0.5*dxi) * ((vel(ihip,j ,0,1)-vel(ihim,j ,0,1))*whi +(vel(ilop,j-1,0,1)-vel(ilom,j-1,0,1))*wlo); - if (j==dhi.y+1 and bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){ - dudx = (dxi) * ((vel(ilop,j-1,0,0)-vel(ilom,j-1,0,0))*wlo); - } - if (j==dhi.y+1 and bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ - dvdx = (dxi) * ((vel(ilop,j-1,0,1)-vel(ilom,j-1,0,1))*wlo); - } - - if (j==dlo.y and bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){ - dudx = (dxi) * ((vel(ihip,j ,0,0)-vel(ihim,j ,0,0))*whi); + if (j==dhi.y+1){ + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudx = (dxi) * ((vel(ilop,j-1,0,0)-vel(ilom,j-1,0,0))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdx = (dxi) * ((vel(ilop,j-1,0,1)-vel(ilom,j-1,0,1))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudx = (dxi) * ((vel(ihip,j ,0,0)-vel(ihim,j ,0,0))*whi); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdx = (dxi) * ((vel(ihip,j ,0,1)-vel(ihim,j ,0,1))*whi); + } } - if (j==dlo.y and bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ - dvdx = (dxi) * ((vel(ihip,j ,0,1)-vel(ihim,j ,0,1))*whi); + if (j==dlo.y){ + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudx = (dxi) * ((vel(ihip,j ,0,0)-vel(ihim,j ,0,0))*whi); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdx = (dxi) * ((vel(ihip,j ,0,1)-vel(ihim,j ,0,1))*whi); + } + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudx = (dxi) * ((vel(ilop,j-1,0,0)-vel(ilom,j-1,0,0))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdx = (dxi) * ((vel(ilop,j-1,0,1)-vel(ilom,j-1,0,1))*wlo); + } } - Real divu = dudx; Real xif = kapy(i,j,0); Real mun = 0.75*(etay(i,j,0,1)-xif); // restore the original eta diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H index 0379da80ecc..6afba949892 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H @@ -55,18 +55,33 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4 const& fx, Real dvdy = (0.5*dyi) * ((vel(i ,jhip,k,1)-vel(i ,jhim,k,1))*whi +(vel(i-1,jlop,k,1)-vel(i-1,jlom,k,1))*wlo); - if (i==dhi.x+1 and bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ - dudy = (dyi) * ((vel(i-1,jlop,k,0)-vel(i-1,jlom,k,0))*wlo); - } - if (i==dhi.x+1 and bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_NEUMANN){ - dvdy = (dyi) * ((vel(i-1,jlop,k,1)-vel(i-1,jlom,k,1))*wlo); - } - - if (i==dlo.x and bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ - dudy = (dyi) * ((vel(i ,jhip,k,0)-vel(i ,jhim,k,0))*whi); + if (i==dhi.x+1){ + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudy = (dyi) * ((vel(i-1,jlop,k,0)-vel(i-1,jlom,k,0))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdy = (dyi) * ((vel(i-1,jlop,k,1)-vel(i-1,jlom,k,1))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudy = (dyi) * ((vel(i ,jhip,k,0)-vel(i ,jhim,k,0))*whi); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdy = (dyi) * ((vel(i ,jhip,k,1)-vel(i ,jhim,k,1))*whi); + } } - if (i==dlo.x and bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_NEUMANN){ - dvdy = (dyi) * ((vel(i ,jhip,k,1)-vel(i ,jhim,k,1))*whi); + if (i==dlo.x){ + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudy = (dyi) * ((vel(i ,jhip,k,0)-vel(i ,jhim,k,0))*whi); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdy = (dyi) * ((vel(i ,jhip,k,1)-vel(i ,jhim,k,1))*whi); + } + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudy = (dyi) * ((vel(i-1,jlop,k,0)-vel(i-1,jlom,k,0))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdy = (dyi) * ((vel(i-1,jlop,k,1)-vel(i-1,jlom,k,1))*wlo); + } } int khip = k + flag(i ,j,k).isConnected(0,0, 1); @@ -82,18 +97,33 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4 const& fx, Real dwdz = (0.5*dzi) * ((vel(i ,j,khip,2)-vel(i ,j,khim,2))*whi +(vel(i-1,j,klop,2)-vel(i-1,j,klom,2))*wlo); - if (i==dhi.x+1 and bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ - dudz = (dzi) * ((vel(i-1,j,klop,0)-vel(i-1,j,klom,0))*wlo); - } - if (i==dhi.x+1 and bct(Orientation(Direction::x,Orientation::high),2)==AMREX_LO_NEUMANN){ - dwdz = (dzi) * ((vel(i-1,j,klop,2)-vel(i-1,j,klom,2))*wlo); - } - - if (i==dlo.x and bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ - dudz = (dzi) * ((vel(i ,j,khip,0)-vel(i ,j,khim,0))*whi); + if (i==dhi.x+1){ + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudz = (dzi) * ((vel(i-1,j,klop,0)-vel(i-1,j,klom,0))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdz = (dzi) * ((vel(i-1,j,klop,2)-vel(i-1,j,klom,2))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudz = (dzi) * ((vel(i ,j,khip,0)-vel(i ,j,khim,0))*whi); + } + if (bct(Orientation(Direction::x,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdz = (dzi) * ((vel(i ,j,khip,2)-vel(i ,j,khim,2))*whi); + } } - if (i==dlo.x and bct(Orientation(Direction::x,Orientation::low),2)==AMREX_LO_NEUMANN){ - dwdz = (dzi) * ((vel(i ,j,khip,2)-vel(i ,j,khim,2))*whi); + if (i==dlo.x){ + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudz = (dzi) * ((vel(i ,j,khip,0)-vel(i ,j,khim,0))*whi); + } + if (bct(Orientation(Direction::x,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdz = (dzi) * ((vel(i ,j,khip,2)-vel(i ,j,khim,2))*whi); + } + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudz = (dzi) * ((vel(i-1,j,klop,0)-vel(i-1,j,klom,0))*wlo); + } + if (bct(Orientation(Direction::x,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdz = (dzi) * ((vel(i-1,j,klop,2)-vel(i-1,j,klom,2))*wlo); + } } Real divu = dvdy + dwdz; @@ -151,18 +181,33 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4 const& fy, Real dvdx = (0.5*dxi) * ((vel(ihip,j ,k,1)-vel(ihim,j ,k,1))*whi +(vel(ilop,j-1,k,1)-vel(ilom,j-1,k,1))*wlo); - if (j==dhi.y+1 and bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){ - dudx = (dxi) * ((vel(ilop,j-1,k,0)-vel(ilom,j-1,k,0))*wlo); - } - if (j==dhi.y+1 and bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ - dvdx = (dxi) * ((vel(ilop,j-1,k,1)-vel(ilom,j-1,k,1))*wlo); - } - - if (j==dlo.y and bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){ - dudx = (dxi) * ((vel(ihip,j ,k,0)-vel(ihim,j ,k,0))*whi); + if (j==dhi.y+1){ + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudx = (dxi) * ((vel(ilop,j-1,k,0)-vel(ilom,j-1,k,0))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdx = (dxi) * ((vel(ilop,j-1,k,1)-vel(ilom,j-1,k,1))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudx = (dxi) * ((vel(ihip,j ,k,0)-vel(ihim,j ,k,0))*whi); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdx = (dxi) * ((vel(ihip,j ,k,1)-vel(ihim,j ,k,1))*whi); + } } - if (j==dlo.y and bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ - dvdx = (dxi) * ((vel(ihip,j ,k,1)-vel(ihim,j ,k,1))*whi); + if (j==dlo.y){ + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudx = (dxi) * ((vel(ihip,j ,k,0)-vel(ihim,j ,k,0))*whi); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdx = (dxi) * ((vel(ihip,j ,k,1)-vel(ihim,j ,k,1))*whi); + } + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudx = (dxi) * ((vel(ilop,j-1,k,0)-vel(ilom,j-1,k,0))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdx = (dxi) * ((vel(ilop,j-1,k,1)-vel(ilom,j-1,k,1))*wlo); + } } int khip = k + flag(i,j ,k).isConnected(0,0, 1); @@ -178,18 +223,33 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4 const& fy, Real dwdz = (0.5*dzi) * ((vel(i,j ,khip,2)-vel(i,j ,khim,2))*whi +(vel(i,j-1,klop,2)-vel(i,j-1,klom,2))*wlo); - if (j==dhi.y+1 and bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ - dvdz = (dzi) * ((vel(i,j-1,klop,1)-vel(i,j-1,klom,1))*wlo); - } - if (j==dhi.y+1 and bct(Orientation(Direction::y,Orientation::high),2)==AMREX_LO_NEUMANN){ - dwdz = (dzi) * ((vel(i,j-1,klop,2)-vel(i,j-1,klom,2))*wlo); - } - - if (j==dlo.y and bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ - dvdz = (dzi) * ((vel(i,j ,khip,1)-vel(i,j ,khim,1))*whi); + if (j==dhi.y+1){ + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdz = (dzi) * ((vel(i,j-1,klop,1)-vel(i,j-1,klom,1))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdz = (dzi) * ((vel(i,j-1,klop,2)-vel(i,j-1,klom,2))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdz = (dzi) * ((vel(i,j ,khip,1)-vel(i,j ,khim,1))*whi); + } + if (bct(Orientation(Direction::y,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdz = (dzi) * ((vel(i,j ,khip,2)-vel(i,j ,khim,2))*whi); + } } - if (j==dlo.y and bct(Orientation(Direction::y,Orientation::low),2)==AMREX_LO_NEUMANN){ - dwdz = (dzi) * ((vel(i,j ,khip,2)-vel(i,j ,khim,2))*whi); + if (j==dlo.y){ + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdz = (dzi) * ((vel(i,j ,khip,1)-vel(i,j ,khim,1))*whi); + } + if (bct(Orientation(Direction::y,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdz = (dzi) * ((vel(i,j ,khip,2)-vel(i,j ,khim,2))*whi); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdz = (dzi) * ((vel(i,j-1,klop,1)-vel(i,j-1,klom,1))*wlo); + } + if (bct(Orientation(Direction::y,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdz = (dzi) * ((vel(i,j-1,klop,2)-vel(i,j-1,klom,2))*wlo); + } } Real divu = dudx + dwdz; @@ -247,18 +307,33 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4 const& fz, Real dwdx = (0.5*dxi) * ((vel(ihip,j,k ,2)-vel(ihim,j,k ,2))*whi +(vel(ilop,j,k-1,2)-vel(ilom,j,k-1,2))*wlo); - if (k==dhi.z+1 and bct(Orientation(Direction::z,Orientation::high),0)==AMREX_LO_NEUMANN){ - dudx = (dxi) * ((vel(ilop,j,k-1,0)-vel(ilom,j,k-1,0))*wlo); - } - if (k==dhi.z+1 and bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_NEUMANN){ - dwdx = (dxi) * ((vel(ilop,j,k-1,2)-vel(ilom,j,k-1,2))*wlo); - } - - if (k==dlo.z and bct(Orientation(Direction::z,Orientation::low),0)==AMREX_LO_NEUMANN){ - dudx = (dxi) * ((vel(ihip,j,k ,0)-vel(ihim,j,k ,0))*whi); + if (k==dhi.z+1){ + if (bct(Orientation(Direction::z,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudx = (dxi) * ((vel(ilop,j,k-1,0)-vel(ilom,j,k-1,0))*wlo); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdx = (dxi) * ((vel(ilop,j,k-1,2)-vel(ilom,j,k-1,2))*wlo); + } + if (bct(Orientation(Direction::z,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudx = (dxi) * ((vel(ihip,j,k ,0)-vel(ihim,j,k ,0))*whi); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdx = (dxi) * ((vel(ihip,j,k ,2)-vel(ihim,j,k ,2))*whi); + } } - if (k==dlo.z and bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_NEUMANN){ - dwdx = (dxi) * ((vel(ihip,j,k ,2)-vel(ihim,j,k ,2))*whi); + if (k==dlo.z){ + if (bct(Orientation(Direction::z,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudx = (dxi) * ((vel(ihip,j,k ,0)-vel(ihim,j,k ,0))*whi); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdx = (dxi) * ((vel(ihip,j,k ,2)-vel(ihim,j,k ,2))*whi); + } + if (bct(Orientation(Direction::z,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudx = (dxi) * ((vel(ilop,j,k-1,0)-vel(ilom,j,k-1,0))*wlo); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdx = (dxi) * ((vel(ilop,j,k-1,2)-vel(ilom,j,k-1,2))*wlo); + } } int jhip = j + flag(i,j,k ).isConnected(0, 1,0); @@ -274,18 +349,33 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4 const& fz, Real dwdy = (0.5*dyi) * ((vel(i,jhip,k ,2)-vel(i,jhim,k ,2))*whi +(vel(i,jlop,k-1,2)-vel(i,jlom,k-1,2))*wlo); - if (k==dhi.z+1 and bct(Orientation(Direction::z,Orientation::high),1)==AMREX_LO_NEUMANN){ - dvdy = (dyi) * ((vel(i,jlop,k-1,1)-vel(i,jlom,k-1,1))*wlo); - } - if (k==dhi.z+1 and bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_NEUMANN){ - dwdy = (dyi) * ((vel(i,jlop,k-1,2)-vel(i,jlom,k-1,2))*wlo); - } - - if (k==dlo.z and bct(Orientation(Direction::z,Orientation::low),1)==AMREX_LO_NEUMANN){ - dvdy = (dyi) * ((vel(i,jhip,k ,1)-vel(i,jhim,k ,1))*whi); + if (k==dhi.z+1){ + if (bct(Orientation(Direction::z,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdy = (dyi) * ((vel(i,jlop,k-1,1)-vel(i,jlom,k-1,1))*wlo); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdy = (dyi) * ((vel(i,jlop,k-1,2)-vel(i,jlom,k-1,2))*wlo); + } + if (bct(Orientation(Direction::z,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdy = (dyi) * ((vel(i,jhip,k ,1)-vel(i,jhim,k ,1))*whi); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdy = (dyi) * ((vel(i,jhip,k ,2)-vel(i,jhim,k ,2))*whi); + } } - if (k==dlo.z and bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_NEUMANN){ - dwdy = (dyi) * ((vel(i,jhip,k ,2)-vel(i,jhim,k ,2))*whi); + if (k==dlo.z){ + if (bct(Orientation(Direction::z,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdy = (dyi) * ((vel(i,jhip,k ,1)-vel(i,jhim,k ,1))*whi); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdy = (dyi) * ((vel(i,jhip,k ,2)-vel(i,jhim,k ,2))*whi); + } + if (bct(Orientation(Direction::z,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdy = (dyi) * ((vel(i,jlop,k-1,1)-vel(i,jlom,k-1,1))*wlo); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdy = (dyi) * ((vel(i,jlop,k-1,2)-vel(i,jlom,k-1,2))*wlo); + } } Real divu = dudx + dvdy; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp index 3030a0f4bba..396069dae33 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp @@ -211,6 +211,11 @@ MLTensorOp::apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc applyBCTensor(amrlev, mglev, in, bc_mode, s_mode, bndry ); const auto dxinv = m_geom[amrlev][mglev].InvCellSizeArray(); + const Geometry& geom = m_geom[amrlev][mglev]; + const Box& domain = geom.Domain(); + const auto dlo = amrex::lbound(domain); + const auto dhi = amrex::ubound(domain); + const auto& bcondloc = *m_bcondloc[amrlev][mglev]; Array const& etamf = m_b_coeffs[amrlev][mglev]; Array const& kapmf = m_kappa[amrlev][mglev]; @@ -245,18 +250,27 @@ MLTensorOp::apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc Array4 const fyfab = fluxfab_tmp[1].array();, Array4 const fzfab = fluxfab_tmp[2].array();); + const auto & bdcv = bcondloc.bndryConds(mfi); + Array2D bct; + for (int icomp = 0; icomp < AMREX_SPACEDIM; ++icomp) { + for (OrientationIter face; face; ++face) { + Orientation ori = face(); + bct(ori,icomp) = bdcv[icomp][ori]; + } + } + AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM ( xbx, txbx, { - mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv); + mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct); } , ybx, tybx, { - mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv); + mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct); } , zbx, tzbx, { - mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv); + mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct); } ); @@ -403,6 +417,11 @@ MLTensorOp::compFlux (int amrlev, const Array& fluxes, applyBCTensor(amrlev, mglev, sol, BCMode::Inhomogeneous, StateMode::Solution, m_bndry_sol[amrlev].get()); const auto dxinv = m_geom[amrlev][mglev].InvCellSizeArray(); + const Geometry& geom = m_geom[amrlev][mglev]; + const Box& domain = geom.Domain(); + const auto dlo = amrex::lbound(domain); + const auto dhi = amrex::ubound(domain); + const auto& bcondloc = *m_bcondloc[amrlev][mglev]; Array const& etamf = m_b_coeffs[amrlev][mglev]; Array const& kapmf = m_kappa[amrlev][mglev]; @@ -436,18 +455,27 @@ MLTensorOp::compFlux (int amrlev, const Array& fluxes, Array4 const fyfab = fluxfab_tmp[1].array();, Array4 const fzfab = fluxfab_tmp[2].array();); + const auto & bdcv = bcondloc.bndryConds(mfi); + Array2D bct; + for (int icomp = 0; icomp < AMREX_SPACEDIM; ++icomp) { + for (OrientationIter face; face; ++face) { + Orientation ori = face(); + bct(ori,icomp) = bdcv[icomp][ori]; + } + } + AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM ( xbx, txbx, { - mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv); + mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct); } , ybx, tybx, { - mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv); + mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct); } , zbx, tzbx, { - mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv); + mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct); } ); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H index 3cc7b4caead..716af64a34e 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H @@ -129,7 +129,9 @@ void mltensor_cross_terms_fx (Box const& box, Array4 const& fx, Array4 const& vel, Array4 const& etax, Array4 const& kapx, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + const Dim3& dlo, const Dim3& dhi, + Array2D const& bct) noexcept { const Real dyi = dxinv[1]; const auto lo = amrex::lbound(box); @@ -141,6 +143,34 @@ void mltensor_cross_terms_fx (Box const& box, Array4 const& fx, for (int i = lo.x; i <= hi.x; ++i) { Real dudy = (vel(i,j+1,0,0)+vel(i-1,j+1,0,0)-vel(i,j-1,0,0)-vel(i-1,j-1,0,0))*(0.25*dyi); Real dvdy = (vel(i,j+1,0,1)+vel(i-1,j+1,0,1)-vel(i,j-1,0,1)-vel(i-1,j-1,0,1))*(0.25*dyi); + if (i==dhi.x+1){ + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudy = (0.5*dyi) * (vel(i-1,j+1,0,0)-vel(i-1,j-1,0,0)); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdy = (0.5*dyi) * (vel(i-1,j+1,0,1)-vel(i-1,j-1,0,1)); + } + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudy = (0.5*dyi) * (vel(i ,j+1,0,0)-vel(i ,j-1,0,0)); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdy = (0.5*dyi) * (vel(i ,j+1,0,1)-vel(i ,j-1,0,1)); + } + } + if (i==dlo.x){ + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudy = (0.5*dyi) * (vel(i ,j+1,0,0)-vel(i ,j-1,0,0)); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdy = (0.5*dyi) * (vel(i ,j+1,0,1)-vel(i ,j-1,0,1)); + } + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudy = (0.5*dyi) * (vel(i-1,j+1,0,0)-vel(i-1,j-1,0,0)); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdy = (0.5*dyi) * (vel(i-1,j+1,0,1)-vel(i-1,j-1,0,1)); + } + } Real divu = dvdy; Real xif = kapx(i,j,0); Real mun = 0.75*(etax(i,j,0,0)-xif); // restore the original eta @@ -156,7 +186,9 @@ void mltensor_cross_terms_fy (Box const& box, Array4 const& fy, Array4 const& vel, Array4 const& etay, Array4 const& kapy, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + const Dim3& dlo, const Dim3& dhi, + Array2D const& bct) noexcept { const Real dxi = dxinv[0]; const auto lo = amrex::lbound(box); @@ -168,6 +200,35 @@ void mltensor_cross_terms_fy (Box const& box, Array4 const& fy, for (int i = lo.x; i <= hi.x; ++i) { Real dudx = (vel(i+1,j,0,0)+vel(i+1,j-1,0,0)-vel(i-1,j,0,0)-vel(i-1,j-1,0,0))*(0.25*dxi); Real dvdx = (vel(i+1,j,0,1)+vel(i+1,j-1,0,1)-vel(i-1,j,0,1)-vel(i-1,j-1,0,1))*(0.25*dxi); + if (j==dhi.y+1){ + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudx = (0.5*dxi) * (vel(i+1,j-1,0,0)-vel(i-1,j-1,0,0)); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdx = (0.5*dxi) * (vel(i+1,j-1,0,1)-vel(i-1,j-1,0,1)); + } + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudx = (0.5*dxi) * (vel(i+1,j ,0,0)-vel(i-1,j ,0,0)); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdx = (0.5*dxi) * (vel(i+1,j ,0,1)-vel(i-1,j ,0,1)); + } + } + if (j==dlo.y){ + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudx = (0.5*dxi) * (vel(i+1,j ,0,0)-vel(i-1,j ,0,0)); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdx = (0.5*dxi) * (vel(i+1,j ,0,1)-vel(i-1,j ,0,1)); + } + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudx = (0.5*dxi) * (vel(i+1,j-1,0,0)-vel(i-1,j-1,0,0)); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdx = (0.5*dxi) * (vel(i+1,j-1,0,1)-vel(i-1,j-1,0,1)); + } + } + Real divu = dudx; Real xif = kapy(i,j,0); Real mun = 0.75*(etay(i,j,0,1)-xif); // restore the original eta diff --git a/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H index 52e46191dcf..c10c1c3e7da 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H @@ -1056,7 +1056,9 @@ void mltensor_cross_terms_fx (Box const& box, Array4 const& fx, Array4 const& vel, Array4 const& etax, Array4 const& kapx, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + const Dim3& dlo, const Dim3& dhi, + Array2D const& bct) noexcept { const Real dyi = dxinv[1]; const Real dzi = dxinv[2]; @@ -1070,8 +1072,65 @@ void mltensor_cross_terms_fx (Box const& box, Array4 const& fx, for (int i = lo.x; i <= hi.x; ++i) { Real dudy = (vel(i,j+1,k,0)+vel(i-1,j+1,k,0)-vel(i,j-1,k,0)-vel(i-1,j-1,k,0))*(0.25*dyi); Real dvdy = (vel(i,j+1,k,1)+vel(i-1,j+1,k,1)-vel(i,j-1,k,1)-vel(i-1,j-1,k,1))*(0.25*dyi); + if (i==dhi.x+1){ + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudy = (0.5*dyi) * (vel(i-1,j+1,k,0)-vel(i-1,j-1,k,0)); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdy = (0.5*dyi) * (vel(i-1,j+1,k,1)-vel(i-1,j-1,k,1)); + } + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudy = (0.5*dyi) * (vel(i ,j+1,k,0)-vel(i ,j-1,k,0)); + } + if (bct(Orientation(Direction::x,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdy = (0.5*dyi) * (vel(i ,j+1,k,1)-vel(i ,j-1,k,1)); + } + } + if (i==dlo.x){ + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudy = (0.5*dyi) * (vel(i ,j+1,k,0)-vel(i ,j-1,k,0)); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdy = (0.5*dyi) * (vel(i ,j+1,k,1)-vel(i ,j-1,k,1)); + } + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudy = (0.5*dyi) * (vel(i-1,j+1,k,0)-vel(i-1,j-1,k,0)); + } + if (bct(Orientation(Direction::x,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdy = (0.5*dyi) * (vel(i-1,j+1,k,1)-vel(i-1,j-1,k,1)); + } + } Real dudz = (vel(i,j,k+1,0)+vel(i-1,j,k+1,0)-vel(i,j,k-1,0)-vel(i-1,j,k-1,0))*(0.25*dzi); Real dwdz = (vel(i,j,k+1,2)+vel(i-1,j,k+1,2)-vel(i,j,k-1,2)-vel(i-1,j,k-1,2))*(0.25*dzi); + if (i==dhi.x+1){ + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudz = (0.5*dzi) * (vel(i-1,j,k+1,0)-vel(i-1,j,k-1,0)); + } + if (bct(Orientation(Direction::x,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdz = (0.5*dzi) * (vel(i-1,j,k+1,2)-vel(i-1,j,k-1,2)); + } + if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudz = (0.5*dzi) * (vel(i ,j,k+1,0)-vel(i ,j,k-1,0)); + } + if (bct(Orientation(Direction::x,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdz = (0.5*dzi) * (vel(i ,j,k+1,2)-vel(i ,j,k-1,2)); + } + } + if (i==dlo.x){ + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudz = (0.5*dzi) * (vel(i ,j,k+1,0)-vel(i ,j,k-1,0)); + } + if (bct(Orientation(Direction::x,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdz = (0.5*dzi) * (vel(i ,j,k+1,2)-vel(i ,j,k-1,2)); + } + if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudz = (0.5*dzi) * (vel(i-1,j,k+1,0)-vel(i-1,j,k-1,0)); + } + if (bct(Orientation(Direction::x,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdz = (0.5*dzi) * (vel(i-1,j,k+1,2)-vel(i-1,j,k-1,2)); + } + } + Real divu = dvdy + dwdz; Real xif = kapx(i,j,k); Real mun = 0.75*(etax(i,j,k,0)-xif); // restore the original eta @@ -1089,7 +1148,9 @@ void mltensor_cross_terms_fy (Box const& box, Array4 const& fy, Array4 const& vel, Array4 const& etay, Array4 const& kapy, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + const Dim3& dlo, const Dim3& dhi, + Array2D const& bct) noexcept { const Real dxi = dxinv[0]; const Real dzi = dxinv[2]; @@ -1103,8 +1164,64 @@ void mltensor_cross_terms_fy (Box const& box, Array4 const& fy, for (int i = lo.x; i <= hi.x; ++i) { Real dudx = (vel(i+1,j,k,0)+vel(i+1,j-1,k,0)-vel(i-1,j,k,0)-vel(i-1,j-1,k,0))*(0.25*dxi); Real dvdx = (vel(i+1,j,k,1)+vel(i+1,j-1,k,1)-vel(i-1,j,k,1)-vel(i-1,j-1,k,1))*(0.25*dxi); + if (j==dhi.y+1){ + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudx = (0.5*dxi) * (vel(i+1,j-1,k,0)-vel(i-1,j-1,k,0)); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdx = (0.5*dxi) * (vel(i+1,j-1,k,1)-vel(i-1,j-1,k,1)); + } + if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudx = (0.5*dxi) * (vel(i+1,j ,k,0)-vel(i-1,j ,k,0)); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdx = (0.5*dxi) * (vel(i+1,j ,k,1)-vel(i-1,j ,k,1)); + } + } + if (j==dlo.y){ + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudx = (0.5*dxi) * (vel(i+1,j ,k,0)-vel(i-1,j ,k,0)); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdx = (0.5*dxi) * (vel(i+1,j ,k,1)-vel(i-1,j ,k,1)); + } + if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudx = (0.5*dxi) * (vel(i+1,j-1,k,0)-vel(i-1,j-1,k,0)); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdx = (0.5*dxi) * (vel(i+1,j-1,k,1)-vel(i-1,j-1,k,1)); + } + } Real dvdz = (vel(i,j,k+1,1)+vel(i,j-1,k+1,1)-vel(i,j,k-1,1)-vel(i,j-1,k-1,1))*(0.25*dzi); Real dwdz = (vel(i,j,k+1,2)+vel(i,j-1,k+1,2)-vel(i,j,k-1,2)-vel(i,j-1,k-1,2))*(0.25*dzi); + if (j==dhi.y+1){ + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdz = (0.5*dzi) * (vel(i,j-1,k+1,1)-vel(i,j-1,k-1,1)); + } + if (bct(Orientation(Direction::y,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdz = (0.5*dzi) * (vel(i,j-1,k+1,2)-vel(i,j-1,k-1,2)); + } + if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdz = (0.5*dzi) * (vel(i,j ,k+1,1)-vel(i,j ,k-1,1)); + } + if (bct(Orientation(Direction::y,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdz = (0.5*dzi) * (vel(i,j ,k+1,2)-vel(i,j ,k-1,2)); + } + } + if (j==dlo.y){ + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdz = (0.5*dzi) * (vel(i,j ,k+1,1)-vel(i,j ,k-1,1)); + } + if (bct(Orientation(Direction::y,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdz = (0.5*dzi) * (vel(i,j ,k+1,2)-vel(i,j ,k-1,2)); + } + if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdz = (0.5*dzi) * (vel(i,j-1,k+1,1)-vel(i,j-1,k-1,1)); + } + if (bct(Orientation(Direction::y,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdz = (0.5*dzi) * (vel(i,j-1,k+1,2)-vel(i,j-1,k-1,2)); + } + } Real divu = dudx + dwdz; Real xif = kapy(i,j,k); Real mun = 0.75*(etay(i,j,k,1)-xif); // restore the original eta @@ -1122,7 +1239,9 @@ void mltensor_cross_terms_fz (Box const& box, Array4 const& fz, Array4 const& vel, Array4 const& etaz, Array4 const& kapz, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + const Dim3& dlo, const Dim3& dhi, + Array2D const& bct) noexcept { const Real dxi = dxinv[0]; const Real dyi = dxinv[1]; @@ -1136,8 +1255,65 @@ void mltensor_cross_terms_fz (Box const& box, Array4 const& fz, for (int i = lo.x; i <= hi.x; ++i) { Real dudx = (vel(i+1,j,k,0)+vel(i+1,j,k-1,0)-vel(i-1,j,k,0)-vel(i-1,j,k-1,0))*(0.25*dxi); Real dwdx = (vel(i+1,j,k,2)+vel(i+1,j,k-1,2)-vel(i-1,j,k,2)-vel(i-1,j,k-1,2))*(0.25*dxi); + if (k==dhi.z+1){ + if (bct(Orientation(Direction::z,Orientation::high),0)==AMREX_LO_NEUMANN){ + dudx = (0.5*dxi) * (vel(i+1,j,k-1,0)-vel(i-1,j,k-1,0)); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdx = (0.5*dxi) * (vel(i+1,j,k-1,2)-vel(i-1,j,k-1,2)); + } + if (bct(Orientation(Direction::z,Orientation::high),0)==AMREX_LO_DIRICHLET){ + dudx = (0.5*dxi) * (vel(i+1,j,k ,0)-vel(i-1,j,k ,0)); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdx = (0.5*dxi) * (vel(i+1,j,k ,2)-vel(i-1,j,k ,2)); + } + } + if (k==dlo.z){ + if (bct(Orientation(Direction::z,Orientation::low),0)==AMREX_LO_NEUMANN){ + dudx = (0.5*dxi) * (vel(i+1,j,k ,0)-vel(i-1,j,k ,0)); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdx = (0.5*dxi) * (vel(i+1,j,k ,2)-vel(i-1,j,k ,2)); + } + if (bct(Orientation(Direction::z,Orientation::low),0)==AMREX_LO_DIRICHLET){ + dudx = (0.5*dxi) * (vel(i+1,j,k-1,0)-vel(i-1,j,k-1,0)); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdx = (0.5*dxi) * (vel(i+1,j,k-1,2)-vel(i-1,j,k-1,2)); + } + } Real dvdy = (vel(i,j+1,k,1)+vel(i,j+1,k-1,1)-vel(i,j-1,k,1)-vel(i,j-1,k-1,1))*(0.25*dyi); Real dwdy = (vel(i,j+1,k,2)+vel(i,j+1,k-1,2)-vel(i,j-1,k,2)-vel(i,j-1,k-1,2))*(0.25*dyi); + if (k==dhi.z+1){ + if (bct(Orientation(Direction::z,Orientation::high),1)==AMREX_LO_NEUMANN){ + dvdy = (0.5*dyi) * (vel(i,j+1,k-1,1)-vel(i,j-1,k-1,1)); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_NEUMANN){ + dwdy = (0.5*dyi) * (vel(i,j+1,k-1,2)-vel(i,j-1,k-1,2)); + } + if (bct(Orientation(Direction::z,Orientation::high),1)==AMREX_LO_DIRICHLET){ + dvdy = (0.5*dyi) * (vel(i,j+1,k ,1)-vel(i,j-1,k ,1)); + } + if (bct(Orientation(Direction::z,Orientation::high),2)==AMREX_LO_DIRICHLET){ + dwdy = (0.5*dyi) * (vel(i,j+1,k ,2)-vel(i,j-1,k ,2)); + } + } + if (k==dlo.z){ + if (bct(Orientation(Direction::z,Orientation::low),1)==AMREX_LO_NEUMANN){ + dvdy = (0.5*dyi) * (vel(i,j+1,k ,1)-vel(i,j-1,k ,1)); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_NEUMANN){ + dwdy = (0.5*dyi) * (vel(i,j+1,k ,2)-vel(i,j-1,k ,2)); + } + if (bct(Orientation(Direction::z,Orientation::low),1)==AMREX_LO_DIRICHLET){ + dvdy = (0.5*dyi) * (vel(i,j+1,k-1,1)-vel(i,j-1,k-1,1)); + } + if (bct(Orientation(Direction::z,Orientation::low),2)==AMREX_LO_DIRICHLET){ + dwdy = (0.5*dyi) * (vel(i,j+1,k-1,2)-vel(i,j-1,k-1,2)); + } + } + Real divu = dudx + dvdy; Real xif = kapz(i,j,k); Real mun = 0.75*(etaz(i,j,k,2)-xif); // restore the original eta