Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Changing computation at inflow/outflow for tensor solve #1235

Merged
merged 9 commits into from
Aug 10, 2020
29 changes: 14 additions & 15 deletions Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down Expand Up @@ -479,21 +479,30 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const
AMREX_D_TERM(Array4<Real const> const kapxfab = kapmf[0].const_array(mfi);,
Array4<Real const> const kapyfab = kapmf[1].const_array(mfi);,
Array4<Real const> const kapzfab = kapmf[2].const_array(mfi););


const auto & bdcv = bcondloc.bndryConds(mfi);
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> 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);
}
);
}
Expand All @@ -504,16 +513,6 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const
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,
{
Expand Down
75 changes: 52 additions & 23 deletions Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H
Original file line number Diff line number Diff line change
Expand Up @@ -52,18 +52,33 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> 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;
Expand Down Expand Up @@ -116,21 +131,35 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> 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
Expand Down
Loading