Skip to content

Commit

Permalink
Changing computation at outflow for the EBTensor (AMReX-Codes#1187)
Browse files Browse the repository at this point in the history 
Computation at outflow boundaries for the EBTensor has been changed to not take into account cells outside the domain.
  • Loading branch information
@OscarAntepara
OscarAntepara authored Aug 2, 2020
1 parent ccaebcc commit 9cd5388
Show file tree
Hide file tree
Showing 3 changed files with 164 additions and 9 deletions.
22 changes: 18 additions & 4 deletions Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -417,6 +417,10 @@ 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 auto dlo = amrex::lbound(domain);
const auto dhi = amrex::ubound(domain);
const auto& bcondloc = *m_bcondloc[amrlev][mglev];

Array<MultiFab,AMREX_SPACEDIM> const& etamf = m_b_coeffs[amrlev][mglev];
Array<MultiFab,AMREX_SPACEDIM> const& kapmf = m_kappa[amrlev][mglev];
Expand Down Expand Up @@ -499,19 +503,29 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const
Array4<Real const> const& apy = area[1]->const_array(mfi);,
Array4<Real const> const& apz = area[2]->const_array(mfi););
Array4<EBCellFlag const> const& flag = flags->const_array(mfi);


const auto & bdcv = bcondloc.bndryConds(mfi);

Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> 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,
{
mlebtensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,apx,flag,dxinv);
mlebtensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,apx,flag,dxinv,dlo,dhi,bct);
}
, ybx, tybx,
{
mlebtensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,apy,flag,dxinv);
mlebtensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,apy,flag,dxinv,dlo,dhi,bct);
}
, zbx, tzbx,
{
mlebtensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,apz,flag,dxinv);
mlebtensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,apz,flag,dxinv,dlo,dhi,bct);
}
);
}
Expand Down
43 changes: 41 additions & 2 deletions Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,9 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
Array4<Real const> const& kapx,
Array4<Real const> const& apx,
Array4<EBCellFlag const> const& flag,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
const Dim3& dlo, const Dim3& dhi,
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
{
const Real dyi = dxinv[1];
const auto lo = amrex::lbound(box);
Expand All @@ -40,12 +42,30 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);

Real whi = mlebtensor_weight(jhip-jhim);
Real wlo = mlebtensor_weight(jlop-jlom);

Real dudy = (0.5*dyi) * ((vel(i ,jhip,0,0)-vel(i ,jhim,0,0))*whi
+(vel(i-1,jlop,0,0)-vel(i-1,jlom,0,0))*wlo);

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==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);
}

Real divu = dvdy;
Real xif = kapx(i,j,0);
Real mun = 0.75*(etax(i,j,0,0)-xif); // restore the original eta
Expand All @@ -64,7 +84,9 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
Array4<Real const> const& kapy,
Array4<Real const> const& apy,
Array4<EBCellFlag const> const& flag,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
const Dim3& dlo, const Dim3& dhi,
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
{
const Real dxi = dxinv[0];
const auto lo = amrex::lbound(box);
Expand All @@ -85,13 +107,30 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);

Real whi = mlebtensor_weight(ihip-ihim);
Real wlo = mlebtensor_weight(ilop-ilom);

Real dudx = (0.5*dxi) * ((vel(ihip,j ,0,0)-vel(ihim,j ,0,0))*whi
+(vel(ilop,j-1,0,0)-vel(ilom,j-1,0,0))*wlo);
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==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);
}


Real divu = dudx;
Real xif = kapy(i,j,0);
Real mun = 0.75*(etay(i,j,0,1)-xif); // restore the original eta
Expand Down
108 changes: 105 additions & 3 deletions Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,9 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
Array4<Real const> const& kapx,
Array4<Real const> const& apx,
Array4<EBCellFlag const> const& flag,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
const Dim3& dlo, const Dim3& dhi,
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
{
const Real dyi = dxinv[1];
const Real dzi = dxinv[2];
Expand All @@ -39,28 +41,61 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
}
else
{

int jhip = j + flag(i ,j,k).isConnected(0, 1,0);
int jhim = j - flag(i ,j,k).isConnected(0,-1,0);
int jlop = j + flag(i-1,j,k).isConnected(0, 1,0);
int jlom = j - flag(i-1,j,k).isConnected(0,-1,0);

Real whi = mlebtensor_weight(jhip-jhim);
Real wlo = mlebtensor_weight(jlop-jlom);

Real dudy = (0.5*dyi) * ((vel(i ,jhip,k,0)-vel(i ,jhim,k,0))*whi
+(vel(i-1,jlop,k,0)-vel(i-1,jlom,k,0))*wlo);
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==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);
}

int khip = k + flag(i ,j,k).isConnected(0,0, 1);
int khim = k - flag(i ,j,k).isConnected(0,0,-1);
int klop = k + flag(i-1,j,k).isConnected(0,0, 1);
int klom = k - flag(i-1,j,k).isConnected(0,0,-1);

whi = mlebtensor_weight(khip-khim);
wlo = mlebtensor_weight(klop-klom);

Real dudz = (0.5*dzi) * ((vel(i ,j,khip,0)-vel(i ,j,khim,0))*whi
+(vel(i-1,j,klop,0)-vel(i-1,j,klom,0))*wlo);
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==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);
}

Real divu = dvdy + dwdz;
Real xif = kapx(i,j,k);
Real mun = 0.75*(etax(i,j,k,0)-xif); // restore the original eta
Expand All @@ -81,7 +116,9 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
Array4<Real const> const& kapy,
Array4<Real const> const& apy,
Array4<EBCellFlag const> const& flag,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
const Dim3& dlo, const Dim3& dhi,
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
{
const Real dxi = dxinv[0];
const Real dzi = dxinv[2];
Expand All @@ -105,24 +142,56 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
int ihim = i - flag(i,j ,k).isConnected(-1,0,0);
int ilop = i + flag(i,j-1,k).isConnected( 1,0,0);
int ilom = i - flag(i,j-1,k).isConnected(-1,0,0);

Real whi = mlebtensor_weight(ihip-ihim);
Real wlo = mlebtensor_weight(ilop-ilom);

Real dudx = (0.5*dxi) * ((vel(ihip,j ,k,0)-vel(ihim,j ,k,0))*whi
+(vel(ilop,j-1,k,0)-vel(ilom,j-1,k,0))*wlo);
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==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);
}

int khip = k + flag(i,j ,k).isConnected(0,0, 1);
int khim = k - flag(i,j ,k).isConnected(0,0,-1);
int klop = k + flag(i,j-1,k).isConnected(0,0, 1);
int klom = k - flag(i,j-1,k).isConnected(0,0,-1);

whi = mlebtensor_weight(khip-khim);
wlo = mlebtensor_weight(klop-klom);

Real dvdz = (0.5*dzi) * ((vel(i,j ,khip,1)-vel(i,j ,khim,1))*whi
+(vel(i,j-1,klop,1)-vel(i,j-1,klom,1))*wlo);
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==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);
}

Real divu = dudx + dwdz;
Real xif = kapy(i,j,k);
Real mun = 0.75*(etay(i,j,k,1)-xif); // restore the original eta
Expand All @@ -143,7 +212,9 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
Array4<Real const> const& kapz,
Array4<Real const> const& apz,
Array4<EBCellFlag const> const& flag,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
const Dim3& dlo, const Dim3& dhi,
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
{
const Real dxi = dxinv[0];
const Real dyi = dxinv[1];
Expand All @@ -167,6 +238,7 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
int ihim = i - flag(i,j,k ).isConnected(-1,0,0);
int ilop = i + flag(i,j,k-1).isConnected( 1,0,0);
int ilom = i - flag(i,j,k-1).isConnected(-1,0,0);

Real whi = mlebtensor_weight(ihip-ihim);
Real wlo = mlebtensor_weight(ilop-ilom);

Expand All @@ -175,17 +247,47 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> 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==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);
}

int jhip = j + flag(i,j,k ).isConnected(0, 1,0);
int jhim = j - flag(i,j,k ).isConnected(0,-1,0);
int jlop = j + flag(i,j,k-1).isConnected(0, 1,0);
int jlom = j - flag(i,j,k-1).isConnected(0,-1,0);

whi = mlebtensor_weight(jhip-jhim);
wlo = mlebtensor_weight(jlop-jlom);

Real dvdy = (0.5*dyi) * ((vel(i,jhip,k ,1)-vel(i,jhim,k ,1))*whi
+(vel(i,jlop,k-1,1)-vel(i,jlom,k-1,1))*wlo);
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==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);
}

Real divu = dudx + dvdy;
Real xif = kapz(i,j,k);
Real mun = 0.75*(etaz(i,j,k,2)-xif); // restore the original eta
Expand Down

0 comments on commit 9cd5388

Please sign in to comment.