Skip to content

Commit

Permalink
Fix the calculation of divu at inflow face in nodal projection so tha…
Browse files Browse the repository at this point in the history 
…t it does not use tangential velocities on an inflow face (AMReX-Codes#1219)

## Summary
Mathematically the nodal projection should only see the normal velocity at inflow faces, but it was coded so that it would use non-zero tangential velocities at inflow in the calculation of divu.  Most users never set tangential velocities to anything but zero, but if they did, this would give an incorrect divergence.
## Additional background

## 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
@asalmgren @dwillcox
asalmgren authored and dwillcox committed Oct 3, 2020
1 parent 20d9831 commit 55f3f7c
Show file tree
Hide file tree
Showing 4 changed files with 257 additions and 80 deletions.
10 changes: 8 additions & 2 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,10 @@ void mlndlap_interpadd_ha (int /*i*/, int /*j*/, int /*k*/, Array4<Real> const&
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlndlap_divu (int /*i*/, int /*j*/, int /*k*/, Array4<Real> const& rhs, Array4<Real const> const& vel,
Array4<int const> const& msk,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
Box const& nddom,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
{}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand All @@ -138,7 +141,10 @@ void mlndlap_mknewu (int /*i*/, int /*j*/, int /*k*/, Array4<Real> const& u, Arr
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlndlap_divu_compute_fine_contrib (int /*i*/, int /*j*/, int /*k*/, Box const& fvbx,
Array4<Real> const& frh, Array4<Real const> const& vel,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
Box const& nddom,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
{}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand Down
105 changes: 89 additions & 16 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeLap_2D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -832,19 +832,43 @@ void mlndlap_interpadd_ha (int i, int j, int,
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlndlap_divu (int i, int j, int k, Array4<Real> const& rhs, Array4<Real const> const& vel,
Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
Box const& nodal_domain,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
bool is_rz) noexcept
{
Real facx = 0.5*dxinv[0];
Real facy = 0.5*dxinv[1];

const auto domlo = amrex::lbound(nodal_domain);
const auto domhi = amrex::ubound(nodal_domain);

if (msk(i,j,k)) {
rhs(i,j,k) = 0.0;
} else {
rhs(i,j,k) = facx*(-vel(i-1,j-1,k,0) + vel(i,j-1,k,0)
-vel(i-1,j ,k,0) + vel(i,j ,k,0))
+ facy*(-vel(i-1,j-1,k,1) - vel(i,j-1,k,1)
+vel(i-1,j ,k,1) + vel(i,j ,k,1));

Real zero_ilo = 1.0;
Real zero_ihi = 1.0;
Real zero_jlo = 1.0;
Real zero_jhi = 1.0;

// The nodal divergence operator should not see the tangential velocity
// at an inflow face
if ((bclo[0] == LinOpBCType::Neumann or bclo[0] == LinOpBCType::inflow)
and i == domlo.x) zero_ilo = 0.0;
if ((bchi[0] == LinOpBCType::Neumann or bchi[0] == LinOpBCType::inflow)
and i == domhi.x) zero_ihi = 0.0;
if ((bclo[1] == LinOpBCType::Neumann or bclo[1] == LinOpBCType::inflow)
and j == domlo.y) zero_jlo = 0.0;
if ((bchi[1] == LinOpBCType::Neumann or bchi[1] == LinOpBCType::inflow)
and j == domhi.y) zero_jhi = 0.0;

rhs(i,j,k) = facx*(-vel(i-1,j-1,k,0)*zero_jlo + vel(i,j-1,k,0)*zero_jlo
-vel(i-1,j ,k,0)*zero_jhi + vel(i,j ,k,0)*zero_jhi)
+ facy*(-vel(i-1,j-1,k,1)*zero_ilo - vel(i,j-1,k,1)*zero_ihi
+vel(i-1,j ,k,1)*zero_ilo + vel(i,j ,k,1)*zero_ihi);
if (is_rz) {
// Here we assume we can't have inflow in the radial direction
Real fm = facy / static_cast<Real>(6*i-3);
Real fp = facy / static_cast<Real>(6*i+3);
rhs(i,j,k) += fm*(vel(i-1,j,k,1)-vel(i-1,j-1,k,1))
Expand Down Expand Up @@ -884,17 +908,42 @@ void mlndlap_mknewu (int i, int j, int k, Array4<Real> const& u, Array4<Real con
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlndlap_divu_compute_fine_contrib (int i, int j, int, Box const& fvbx,
Array4<Real> const& frh, Array4<Real const> const& vel,
bool is_rz, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
Box const& nodal_domain,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
bool is_rz) noexcept
{
const auto domlo = amrex::lbound(nodal_domain);
const auto domhi = amrex::ubound(nodal_domain);

IntVect iv(i,j);
if (fvbx.contains(iv) and !fvbx.strictly_contains(iv))
{
Real zero_ilo = 1.0;
Real zero_ihi = 1.0;
Real zero_jlo = 1.0;
Real zero_jhi = 1.0;

// The nodal divergence operator should not see the tangential velocity
// at an inflow face
if ((bclo[0] == LinOpBCType::Neumann or bclo[0] == LinOpBCType::inflow)
and i == domlo.x) zero_ilo = 0.0;
if ((bchi[0] == LinOpBCType::Neumann or bchi[0] == LinOpBCType::inflow)
and i == domhi.x) zero_ihi = 0.0;
if ((bclo[1] == LinOpBCType::Neumann or bclo[1] == LinOpBCType::inflow)
and j == domlo.y) zero_jlo = 0.0;
if ((bchi[1] == LinOpBCType::Neumann or bchi[1] == LinOpBCType::inflow)
and j == domhi.y) zero_jhi = 0.0;

Real facx = 0.5_rt*dxinv[0];
Real facy = 0.5_rt*dxinv[1];
frh(i,j,0) = facx*(-vel(i-1,j-1,0,0)+vel(i,j-1,0,0)-vel(i-1,j,0,0)+vel(i,j,0,0))
+ facy*(-vel(i-1,j-1,0,1)-vel(i,j-1,0,1)+vel(i-1,j,0,1)+vel(i,j,0,1));

frh(i,j,0) = facx*(-vel(i-1,j-1,0,0)*zero_jlo+vel(i,j-1,0,0)*zero_jlo-vel(i-1,j,0,0)*zero_jhi+vel(i,j,0,0)*zero_jhi)
+ facy*(-vel(i-1,j-1,0,1)*zero_ilo-vel(i,j-1,0,1)*zero_ihi+vel(i-1,j,0,1)*zero_ilo+vel(i,j,0,1)*zero_ihi);

if (is_rz) {
// Here we assume we can't have inflow in the radial direction
Real fm = facy / static_cast<Real>(6*i-3);
Real fp = facy / static_cast<Real>(6*i+3);
frh(i,j,0) += fm*(vel(i-1,j,0,1)-vel(i-1,j-1,0,1))
Expand Down Expand Up @@ -1656,19 +1705,43 @@ void mlndlap_set_stencil_eb (int i, int j, int, Array4<Real> const& sten,
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlndlap_divu_eb (int i, int j, int, Array4<Real> const& rhs, Array4<Real const> const& vel,
Array4<Real const> const& vfrac, Array4<Real const> const& intg,
Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
Box const& nodal_domain,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
{
Real facx = 0.5_rt*dxinv[0];
Real facy = 0.5_rt*dxinv[1];

const auto domlo = amrex::lbound(nodal_domain);
const auto domhi = amrex::ubound(nodal_domain);

if (!msk(i,j,0)) {
rhs(i,j,0) = facx*(-vel(i-1,j-1,0,0)*(vfrac(i-1,j-1,0)+2._rt*intg(i-1,j-1,0,1))
+vel(i ,j-1,0,0)*(vfrac(i ,j-1,0)+2._rt*intg(i ,j-1,0,1))
-vel(i-1,j ,0,0)*(vfrac(i-1,j ,0)-2._rt*intg(i-1,j ,0,1))
+vel(i ,j ,0,0)*(vfrac(i ,j ,0)-2._rt*intg(i ,j ,0,1)))
+ facy*(-vel(i-1,j-1,0,1)*(vfrac(i-1,j-1,0)+2._rt*intg(i-1,j-1,0,0))
-vel(i ,j-1,0,1)*(vfrac(i ,j-1,0)-2._rt*intg(i ,j-1,0,0))
+vel(i-1,j ,0,1)*(vfrac(i-1,j ,0)+2._rt*intg(i-1,j ,0,0))
+vel(i ,j ,0,1)*(vfrac(i ,j ,0)-2._rt*intg(i ,j ,0,0)));

Real zero_ilo = 1.0;
Real zero_ihi = 1.0;
Real zero_jlo = 1.0;
Real zero_jhi = 1.0;

// The nodal divergence operator should not see the tangential velocity
// at an inflow face
if ((bclo[0] == LinOpBCType::Neumann or bclo[0] == LinOpBCType::inflow)
and i == domlo.x) zero_ilo = 0.0;
if ((bchi[0] == LinOpBCType::Neumann or bchi[0] == LinOpBCType::inflow)
and i == domhi.x) zero_ihi = 0.0;
if ((bclo[1] == LinOpBCType::Neumann or bclo[1] == LinOpBCType::inflow)
and j == domlo.y) zero_jlo = 0.0;
if ((bchi[1] == LinOpBCType::Neumann or bchi[1] == LinOpBCType::inflow)
and j == domhi.y) zero_jhi = 0.0;

rhs(i,j,0) = facx*(-vel(i-1,j-1,0,0)*(vfrac(i-1,j-1,0)+2._rt*intg(i-1,j-1,0,1))*zero_jlo
+vel(i ,j-1,0,0)*(vfrac(i ,j-1,0)+2._rt*intg(i ,j-1,0,1))*zero_jlo
-vel(i-1,j ,0,0)*(vfrac(i-1,j ,0)-2._rt*intg(i-1,j ,0,1))*zero_jhi
+vel(i ,j ,0,0)*(vfrac(i ,j ,0)-2._rt*intg(i ,j ,0,1))*zero_jhi)
+ facy*(-vel(i-1,j-1,0,1)*(vfrac(i-1,j-1,0)+2._rt*intg(i-1,j-1,0,0))*zero_ilo
-vel(i ,j-1,0,1)*(vfrac(i ,j-1,0)-2._rt*intg(i ,j-1,0,0))*zero_ihi
+vel(i-1,j ,0,1)*(vfrac(i-1,j ,0)+2._rt*intg(i-1,j ,0,0))*zero_ilo
+vel(i ,j ,0,1)*(vfrac(i ,j ,0)-2._rt*intg(i ,j ,0,0))*zero_ihi);
} else {
rhs(i,j,0) = 0._rt;
}
Expand Down
Loading

0 comments on commit 55f3f7c

Please sign in to comment.