Skip to content

Commit

Permalink
Tensor solver and extension of EB domain face (#1262)
Browse files Browse the repository at this point in the history 
## Summary

ParmParse parameter eb2.extend_domain_face is now true by default.  The
reason for the change is that for most applications it's hard to imagine one
would want it to be false.  It will also solve the EB tensor solver issue
when there is a cut cell just outside the domain abutting a covered valid
cell.

Recent changes to the tensor solvers are incorrect and are reverted here.
They are incorrect because the ghost cells in the modified functions
actually have values at the cell centered.  Although the users put domain
face values in the ghost cells, that is not what the solver does.  The
solver stores the boundary values in boundary registers.  Before the stencil
is being applied, bc function is called to fill the ghost cells with
properly interpolated or extrapolated values so that the stencil operations
are just like working on normal interior cells.

Revert "Correcting tensor cross terms computation for periodic bcs (#1254)"
This reverts commit 1197adc.

Revert "Changing computation at inflow/outflow for tensor solve (#1235)"
This reverts commit b391885.

Revert "Changing computation at outflow for the EBTensor (#1187)"
This reverts commit 9cd5388.

## Checklist

The proposed changes:
- [x] fix a bug or incorrect behavior in AMReX
- [ ] add new capabilities to AMReX
- [x] changes answers in the test suite to more than roundoff level
- [x] are likely to significantly affect the results of downstream AMReX users
- [ ] are described in the proposed changes to the AMReX documentation, if appropriate
  • Loading branch information
@WeiqunZhang
WeiqunZhang authored Aug 14, 2020
1 parent 5f7a294 commit 6c36b5f
Show file tree
Hide file tree
Showing 8 changed files with 29 additions and 576 deletions.
2 changes: 1 addition & 1 deletion Src/EB/AMReX_EB2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ namespace amrex { namespace EB2 {
Vector<std::unique_ptr<IndexSpace> > IndexSpace::m_instance;

int max_grid_size = 64;
bool extend_domain_face = false;
bool extend_domain_face = true;

void Initialize ()
{
Expand Down
4 changes: 4 additions & 0 deletions Src/EB/AMReX_EB2_Level.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,10 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
return;
}

if (amrex::Verbose() > 0 && extend_domain_face == false) {
amrex::Print() << "AMReX WARNING: extend_domain_face=false is not recommended!\n";
}

BL_PROFILE("EB2::GShopLevel()-fine");

Real small_volfrac = 1.e-14;
Expand Down
30 changes: 8 additions & 22 deletions Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -417,11 +417,6 @@ 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 = geom.Domain();
const auto dlo = amrex::lbound(domain);
const auto dhi = amrex::ubound(domain);
const GpuArray<int,AMREX_SPACEDIM>& is_periodic = geom.isPeriodicArray();
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 @@ -480,30 +475,21 @@ 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,dlo,dhi,bct,is_periodic[0]);
mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv);
}
, ybx, tybx,
{
mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct,is_periodic[1]);
mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv);
}
, zbx, tzbx,
{
mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct,is_periodic[2]);
mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv);
}
);
}
Expand All @@ -513,19 +499,19 @@ 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);

AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM
( xbx, txbx,
{
mlebtensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,apx,flag,dxinv,dlo,dhi,bct,is_periodic[0]);
mlebtensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,apx,flag,dxinv);
}
, ybx, tybx,
{
mlebtensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,apy,flag,dxinv,dlo,dhi,bct,is_periodic[1]);
mlebtensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,apy,flag,dxinv);
}
, zbx, tzbx,
{
mlebtensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,apz,flag,dxinv,dlo,dhi,bct,is_periodic[2]);
mlebtensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,apz,flag,dxinv);
}
);
}
Expand Down
74 changes: 2 additions & 72 deletions Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,7 @@ 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,
const Dim3& dlo, const Dim3& dhi,
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
const int& is_periodic) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
{
const Real dyi = dxinv[1];
const auto lo = amrex::lbound(box);
Expand All @@ -43,45 +40,12 @@ 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 !is_periodic){
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 !is_periodic){
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;
Real xif = kapx(i,j,0);
Real mun = 0.75*(etax(i,j,0,0)-xif); // restore the original eta
Expand All @@ -100,10 +64,7 @@ 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,
const Dim3& dlo, const Dim3& dhi,
Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
const int& is_periodic) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
{
const Real dxi = dxinv[0];
const auto lo = amrex::lbound(box);
Expand All @@ -124,44 +85,13 @@ 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 !is_periodic){
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 !is_periodic){
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

0 comments on commit 6c36b5f

Please sign in to comment.