diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H index 90a145ac44e..769e725636a 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H @@ -122,7 +122,10 @@ void mlndlap_interpadd_ha (int /*i*/, int /*j*/, int /*k*/, Array4 const& AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu (int /*i*/, int /*j*/, int /*k*/, Array4 const& rhs, Array4 const& vel, Array4 const& msk, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + Box const& nddom, + GpuArray const& bclo, + GpuArray const& bchi) noexcept {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -138,7 +141,10 @@ void mlndlap_mknewu (int /*i*/, int /*j*/, int /*k*/, Array4 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 const& frh, Array4 const& vel, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + Box const& nddom, + GpuArray const& bclo, + GpuArray const& bchi) noexcept {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_2D_K.H index a7f052e2790..a5c80da4b10 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_2D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_2D_K.H @@ -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 const& rhs, Array4 const& vel, Array4 const& msk, GpuArray const& dxinv, + Box const& nodal_domain, + GpuArray const& bclo, + GpuArray 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(6*i-3); Real fp = facy / static_cast(6*i+3); rhs(i,j,k) += fm*(vel(i-1,j,k,1)-vel(i-1,j-1,k,1)) @@ -884,17 +908,42 @@ void mlndlap_mknewu (int i, int j, int k, Array4 const& u, Array4 const& frh, Array4 const& vel, - bool is_rz, GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + Box const& nodal_domain, + GpuArray const& bclo, + GpuArray 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(6*i-3); Real fp = facy / static_cast(6*i+3); frh(i,j,0) += fm*(vel(i-1,j,0,1)-vel(i-1,j-1,0,1)) @@ -1656,19 +1705,43 @@ void mlndlap_set_stencil_eb (int i, int j, int, Array4 const& sten, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu_eb (int i, int j, int, Array4 const& rhs, Array4 const& vel, Array4 const& vfrac, Array4 const& intg, - Array4 const& msk, GpuArray const& dxinv) noexcept + Array4 const& msk, GpuArray const& dxinv, + Box const& nodal_domain, + GpuArray const& bclo, + GpuArray 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; } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_3D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_3D_K.H index 860f0921a06..9ae27495c93 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_3D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_3D_K.H @@ -1779,27 +1779,58 @@ void mlndlap_interpadd_ha (int i, int j, int k, Array4 const& fine, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu (int i, int j, int k, Array4 const& rhs, Array4 const& vel, Array4 const& msk, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + Box const& nodal_domain, + GpuArray const& bclo, + GpuArray const& bchi) noexcept { Real facx = 0.25*dxinv[0]; Real facy = 0.25*dxinv[1]; Real facz = 0.25*dxinv[2]; + 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-1,0)+vel(i,j-1,k-1,0) - -vel(i-1,j ,k-1,0)+vel(i,j ,k-1,0) - -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,1)-vel(i,j-1,k-1,1) - +vel(i-1,j ,k-1,1)+vel(i,j ,k-1,1) - -vel(i-1,j-1,k ,1)-vel(i,j-1,k ,1) - +vel(i-1,j ,k ,1)+vel(i,j ,k ,1)) - + facz*(-vel(i-1,j-1,k-1,2)-vel(i,j-1,k-1,2) - -vel(i-1,j ,k-1,2)-vel(i,j ,k-1,2) - +vel(i-1,j-1,k ,2)+vel(i,j-1,k ,2) - +vel(i-1,j ,k ,2)+vel(i,j ,k ,2)); + + Real zero_ilo = 1.0; + Real zero_ihi = 1.0; + Real zero_jlo = 1.0; + Real zero_jhi = 1.0; + Real zero_klo = 1.0; + Real zero_khi = 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; + if ((bclo[2] == LinOpBCType::Neumann or bclo[2] == LinOpBCType::inflow) + and k == domlo.z) zero_klo = 0.0; + if ((bchi[2] == LinOpBCType::Neumann or bchi[2] == LinOpBCType::inflow) + and k == domhi.z) zero_khi = 0.0; + + rhs(i,j,k) = facx*(-vel(i-1,j-1,k-1,0)*zero_jlo*zero_klo+vel(i,j-1,k-1,0)*zero_jlo*zero_klo + -vel(i-1,j ,k-1,0)*zero_jhi*zero_klo+vel(i,j ,k-1,0)*zero_jhi*zero_klo + -vel(i-1,j-1,k ,0)*zero_jlo*zero_khi+vel(i,j-1,k ,0)*zero_jlo*zero_khi + -vel(i-1,j ,k ,0)*zero_jhi*zero_khi+vel(i,j ,k ,0)*zero_jhi*zero_khi) + + + facy*(-vel(i-1,j-1,k-1,1)*zero_ilo*zero_klo-vel(i,j-1,k-1,1)*zero_ihi*zero_klo + +vel(i-1,j ,k-1,1)*zero_ilo*zero_klo+vel(i,j ,k-1,1)*zero_ihi*zero_klo + -vel(i-1,j-1,k ,1)*zero_ilo*zero_khi-vel(i,j-1,k ,1)*zero_ihi*zero_khi + +vel(i-1,j ,k ,1)*zero_ilo*zero_khi+vel(i,j ,k ,1)*zero_ihi*zero_khi) + + + facz*(-vel(i-1,j-1,k-1,2)*zero_ilo*zero_jlo-vel(i,j-1,k-1,2)*zero_ihi*zero_jlo + -vel(i-1,j ,k-1,2)*zero_ilo*zero_jhi-vel(i,j ,k-1,2)*zero_ihi*zero_jhi + +vel(i-1,j-1,k ,2)*zero_ilo*zero_jlo+vel(i,j-1,k ,2)*zero_ihi*zero_jlo + +vel(i-1,j ,k ,2)*zero_ilo*zero_jhi+vel(i,j ,k ,2)*zero_ihi*zero_jhi); } } @@ -1838,23 +1869,51 @@ void mlndlap_mknewu (int i, int j, int k, Array4 const& u, Array4 const& frh, Array4 const& vel, - GpuArray const& dxinv) noexcept + GpuArray const& dxinv, + Box const& nodal_domain, + GpuArray const& bclo, + GpuArray const& bchi) noexcept { + const auto domlo = amrex::lbound(nodal_domain); + const auto domhi = amrex::ubound(nodal_domain); + IntVect iv(i,j,k); if (fvbx.contains(iv) and !fvbx.strictly_contains(iv)) { - frh(i,j,k) = 0.25*dxinv[0]*(-vel(i-1,j-1,k-1,0)+vel(i,j-1,k-1,0) - -vel(i-1,j ,k-1,0)+vel(i,j ,k-1,0) - -vel(i-1,j-1,k ,0)+vel(i,j-1,k ,0) - -vel(i-1,j ,k ,0)+vel(i,j ,k ,0)) - + 0.25*dxinv[1]*(-vel(i-1,j-1,k-1,1)-vel(i,j-1,k-1,1) - +vel(i-1,j ,k-1,1)+vel(i,j ,k-1,1) - -vel(i-1,j-1,k ,1)-vel(i,j-1,k ,1) - +vel(i-1,j ,k ,1)+vel(i,j ,k ,1)) - + 0.25*dxinv[2]*(-vel(i-1,j-1,k-1,2)-vel(i,j-1,k-1,2) - -vel(i-1,j ,k-1,2)-vel(i,j ,k-1,2) - +vel(i-1,j-1,k ,2)+vel(i,j-1,k ,2) - +vel(i-1,j ,k ,2)+vel(i,j ,k ,2)); + Real zero_ilo = 1.0; + Real zero_ihi = 1.0; + Real zero_jlo = 1.0; + Real zero_jhi = 1.0; + Real zero_klo = 1.0; + Real zero_khi = 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; + if ((bclo[2] == LinOpBCType::Neumann or bclo[2] == LinOpBCType::inflow) + and k == domlo.z) zero_klo = 0.0; + if ((bchi[2] == LinOpBCType::Neumann or bchi[2] == LinOpBCType::inflow) + and k == domhi.z) zero_khi = 0.0; + + frh(i,j,k) = 0.25*dxinv[0]*(-vel(i-1,j-1,k-1,0)*zero_jlo*zero_klo+vel(i,j-1,k-1,0)*zero_jlo*zero_klo + -vel(i-1,j ,k-1,0)*zero_jhi*zero_klo+vel(i,j ,k-1,0)*zero_jhi*zero_klo + -vel(i-1,j-1,k ,0)*zero_jlo*zero_khi+vel(i,j-1,k ,0)*zero_jlo*zero_khi + -vel(i-1,j ,k ,0)*zero_jhi*zero_khi+vel(i,j ,k ,0)*zero_jhi*zero_khi) + + 0.25*dxinv[1]*(-vel(i-1,j-1,k-1,1)*zero_ilo*zero_klo-vel(i,j-1,k-1,1)*zero_ihi*zero_klo + +vel(i-1,j ,k-1,1)*zero_ilo*zero_klo+vel(i,j ,k-1,1)*zero_ihi*zero_klo + -vel(i-1,j-1,k ,1)*zero_ilo*zero_khi-vel(i,j-1,k ,1)*zero_ihi*zero_khi + +vel(i-1,j ,k ,1)*zero_ilo*zero_khi+vel(i,j ,k ,1)*zero_ihi*zero_khi) + + 0.25*dxinv[2]*(-vel(i-1,j-1,k-1,2)*zero_ilo*zero_jlo-vel(i,j-1,k-1,2)*zero_ihi*zero_jlo + -vel(i-1,j ,k-1,2)*zero_ilo*zero_jhi-vel(i,j ,k-1,2)*zero_ihi*zero_jhi + +vel(i-1,j-1,k ,2)*zero_ilo*zero_jlo+vel(i,j-1,k ,2)*zero_ihi*zero_jlo + +vel(i-1,j ,k ,2)*zero_ilo*zero_jhi+vel(i,j ,k ,2)*zero_ihi*zero_jhi); } } @@ -6238,112 +6297,141 @@ void mlndlap_set_stencil_eb (int i, int j, int k, Array4 const& sten, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu_eb (int i, int j, int k, Array4 const& rhs, Array4 const& vel, Array4 const& vfrac, Array4 const& intg, - Array4 const& msk, GpuArray const& dxinv) noexcept + Array4 const& msk, GpuArray const& dxinv, + Box const& nodal_domain, + GpuArray const& bclo, + GpuArray const& bchi) noexcept { Real facx = 0.25_rt*dxinv[0]; Real facy = 0.25_rt*dxinv[1]; Real facz = 0.25_rt*dxinv[2]; + const auto domlo = amrex::lbound(nodal_domain); + const auto domhi = amrex::ubound(nodal_domain); + if (!msk(i,j,k)) { + + Real zero_ilo = 1.0; + Real zero_ihi = 1.0; + Real zero_jlo = 1.0; + Real zero_jhi = 1.0; + Real zero_klo = 1.0; + Real zero_khi = 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; + if ((bclo[2] == LinOpBCType::Neumann or bclo[2] == LinOpBCType::inflow) + and k == domlo.z) zero_klo = 0.0; + if ((bchi[2] == LinOpBCType::Neumann or bchi[2] == LinOpBCType::inflow) + and k == domhi.z) zero_khi = 0.0; + rhs(i,j,k) = facx*( vel(i-1,j-1,k ,0)*( -vfrac(i-1,j-1,k ) -2._rt*intg(i-1,j-1,k ,i_S_y) +2._rt*intg(i-1,j-1,k ,i_S_z) - +4._rt*intg(i-1,j-1,k ,i_S_y_z)) + +4._rt*intg(i-1,j-1,k ,i_S_y_z))*zero_jlo*zero_khi +vel(i ,j-1,k ,0)*( vfrac(i ,j-1,k ) +2._rt*intg(i ,j-1,k ,i_S_y) -2._rt*intg(i ,j-1,k ,i_S_z) - -4._rt*intg(i ,j-1,k ,i_S_y_z)) + -4._rt*intg(i ,j-1,k ,i_S_y_z))*zero_jlo*zero_khi +vel(i-1,j ,k ,0)*( -vfrac(i-1,j ,k ) +2._rt*intg(i-1,j ,k ,i_S_y) +2._rt*intg(i-1,j ,k ,i_S_z) - -4._rt*intg(i-1,j ,k ,i_S_y_z)) + -4._rt*intg(i-1,j ,k ,i_S_y_z))*zero_jhi*zero_khi +vel(i ,j ,k ,0)*( vfrac(i ,j ,k ) -2._rt*intg(i ,j ,k ,i_S_y) -2._rt*intg(i ,j ,k ,i_S_z) - +4._rt*intg(i ,j ,k ,i_S_y_z)) + +4._rt*intg(i ,j ,k ,i_S_y_z))*zero_jhi*zero_khi +vel(i-1,j-1,k-1,0)*( -vfrac(i-1,j-1,k-1) -2._rt*intg(i-1,j-1,k-1,i_S_y) -2._rt*intg(i-1,j-1,k-1,i_S_z) - -4._rt*intg(i-1,j-1,k-1,i_S_y_z)) + -4._rt*intg(i-1,j-1,k-1,i_S_y_z))*zero_jlo*zero_klo +vel(i ,j-1,k-1,0)*( vfrac(i ,j-1,k-1) +2._rt*intg(i ,j-1,k-1,i_S_y) +2._rt*intg(i ,j-1,k-1,i_S_z) - +4._rt*intg(i ,j-1,k-1,i_S_y_z)) + +4._rt*intg(i ,j-1,k-1,i_S_y_z))*zero_jlo*zero_klo +vel(i-1,j ,k-1,0)*( -vfrac(i-1,j ,k-1) +2._rt*intg(i-1,j ,k-1,i_S_y) -2._rt*intg(i-1,j ,k-1,i_S_z) - +4._rt*intg(i-1,j ,k-1,i_S_y_z)) + +4._rt*intg(i-1,j ,k-1,i_S_y_z))*zero_jhi*zero_klo +vel(i ,j ,k-1,0)*( vfrac(i ,j ,k-1) -2._rt*intg(i ,j ,k-1,i_S_y) +2._rt*intg(i ,j ,k-1,i_S_z) - -4._rt*intg(i ,j ,k-1,i_S_y_z)) ) + -4._rt*intg(i ,j ,k-1,i_S_y_z))*zero_jhi*zero_klo ) + facy*( vel(i-1,j-1,k ,1)*( -vfrac(i-1,j-1,k ) -2._rt*intg(i-1,j-1,k ,i_S_x) +2._rt*intg(i-1,j-1,k ,i_S_z) - +4._rt*intg(i-1,j-1,k ,i_S_x_z)) + +4._rt*intg(i-1,j-1,k ,i_S_x_z))*zero_ilo*zero_khi +vel(i ,j-1,k ,1)*( -vfrac(i ,j-1,k ) +2._rt*intg(i ,j-1,k ,i_S_x) +2._rt*intg(i ,j-1,k ,i_S_z) - -4._rt*intg(i ,j-1,k ,i_S_x_z)) + -4._rt*intg(i ,j-1,k ,i_S_x_z))*zero_ihi*zero_khi +vel(i-1,j ,k ,1)*( vfrac(i-1,j ,k ) +2._rt*intg(i-1,j ,k ,i_S_x) -2._rt*intg(i-1,j ,k ,i_S_z) - -4._rt*intg(i-1,j ,k ,i_S_x_z)) + -4._rt*intg(i-1,j ,k ,i_S_x_z))*zero_ilo*zero_khi +vel(i ,j ,k ,1)*( vfrac(i ,j ,k ) -2._rt*intg(i ,j ,k ,i_S_x) -2._rt*intg(i ,j ,k ,i_S_z) - +4._rt*intg(i ,j ,k ,i_S_x_z)) + +4._rt*intg(i ,j ,k ,i_S_x_z))*zero_ihi*zero_khi +vel(i-1,j-1,k-1,1)*( -vfrac(i-1,j-1,k-1) -2._rt*intg(i-1,j-1,k-1,i_S_x) -2._rt*intg(i-1,j-1,k-1,i_S_z) - -4._rt*intg(i-1,j-1,k-1,i_S_x_z)) + -4._rt*intg(i-1,j-1,k-1,i_S_x_z))*zero_ilo*zero_klo +vel(i ,j-1,k-1,1)*( -vfrac(i ,j-1,k-1) +2._rt*intg(i ,j-1,k-1,i_S_x) -2._rt*intg(i ,j-1,k-1,i_S_z) - +4._rt*intg(i ,j-1,k-1,i_S_x_z)) + +4._rt*intg(i ,j-1,k-1,i_S_x_z))*zero_ihi*zero_klo +vel(i-1,j ,k-1,1)*( vfrac(i-1,j ,k-1) +2._rt*intg(i-1,j ,k-1,i_S_x) +2._rt*intg(i-1,j ,k-1,i_S_z) - +4._rt*intg(i-1,j ,k-1,i_S_x_z)) + +4._rt*intg(i-1,j ,k-1,i_S_x_z))*zero_ilo*zero_klo +vel(i ,j ,k-1,1)*( vfrac(i ,j ,k-1) -2._rt*intg(i ,j ,k-1,i_S_x) +2._rt*intg(i ,j ,k-1,i_S_z) - -4._rt*intg(i ,j ,k-1,i_S_x_z)) ) + -4._rt*intg(i ,j ,k-1,i_S_x_z))*zero_ihi*zero_klo ) + facz*( vel(i-1,j-1,k ,2)*( vfrac(i-1,j-1,k ) +2._rt*intg(i-1,j-1,k ,i_S_x) +2._rt*intg(i-1,j-1,k ,i_S_y) - +4._rt*intg(i-1,j-1,k ,i_S_x_y)) + +4._rt*intg(i-1,j-1,k ,i_S_x_y))*zero_ilo*zero_jlo +vel(i ,j-1,k ,2)*( vfrac(i ,j-1,k ) -2._rt*intg(i ,j-1,k ,i_S_x) +2._rt*intg(i ,j-1,k ,i_S_y) - -4._rt*intg(i ,j-1,k ,i_S_x_y)) + -4._rt*intg(i ,j-1,k ,i_S_x_y))*zero_ihi*zero_jlo +vel(i-1,j ,k ,2)*( vfrac(i-1,j ,k ) +2._rt*intg(i-1,j ,k ,i_S_x) -2._rt*intg(i-1,j ,k ,i_S_y) - -4._rt*intg(i-1,j ,k ,i_S_x_y)) + -4._rt*intg(i-1,j ,k ,i_S_x_y))*zero_ilo*zero_jhi +vel(i ,j ,k ,2)*( vfrac(i ,j ,k ) -2._rt*intg(i ,j ,k ,i_S_x) -2._rt*intg(i ,j ,k ,i_S_y) - +4._rt*intg(i ,j ,k ,i_S_x_y)) + +4._rt*intg(i ,j ,k ,i_S_x_y))*zero_ihi*zero_jhi +vel(i-1,j-1,k-1,2)*( -vfrac(i-1,j-1,k-1) -2._rt*intg(i-1,j-1,k-1,i_S_x) -2._rt*intg(i-1,j-1,k-1,i_S_y) - -4._rt*intg(i-1,j-1,k-1,i_S_x_y)) + -4._rt*intg(i-1,j-1,k-1,i_S_x_y))*zero_ilo*zero_jlo +vel(i ,j-1,k-1,2)*( -vfrac(i ,j-1,k-1) +2._rt*intg(i ,j-1,k-1,i_S_x) -2._rt*intg(i ,j-1,k-1,i_S_y) - +4._rt*intg(i ,j-1,k-1,i_S_x_y)) + +4._rt*intg(i ,j-1,k-1,i_S_x_y))*zero_ihi*zero_jlo +vel(i-1,j ,k-1,2)*( -vfrac(i-1,j ,k-1) -2._rt*intg(i-1,j ,k-1,i_S_x) +2._rt*intg(i-1,j ,k-1,i_S_y) - +4._rt*intg(i-1,j ,k-1,i_S_x_y)) + +4._rt*intg(i-1,j ,k-1,i_S_x_y))*zero_ilo*zero_jhi +vel(i ,j ,k-1,2)*( -vfrac(i ,j ,k-1) +2._rt*intg(i ,j ,k-1,i_S_x) +2._rt*intg(i ,j ,k-1,i_S_y) - -4._rt*intg(i ,j ,k-1,i_S_x_y)) ); + -4._rt*intg(i ,j ,k-1,i_S_x_y))*zero_ihi*zero_jhi ); } else { rhs(i,j,k) = 0._rt; } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp index 70f2986f066..96d4d33c91f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp @@ -230,7 +230,7 @@ MLNodeLaplacian::compRHS (const Vector& rhs, const Vector& Array4 const& intgarr = intg->const_array(mfi); AMREX_HOST_DEVICE_FOR_3D(bx, i, j, k, { - mlndlap_divu_eb(i,j,k,rhsarr,velarr,vfracarr,intgarr,dmskarr,dxinvarr); + mlndlap_divu_eb(i,j,k,rhsarr,velarr,vfracarr,intgarr,dmskarr,dxinvarr,nddom,lobc,hibc); }); } else @@ -244,12 +244,14 @@ MLNodeLaplacian::compRHS (const Vector& rhs, const Vector& #if (AMREX_SPACEDIM == 2) AMREX_HOST_DEVICE_PARALLEL_FOR_3D (bx, i, j, k, { - mlndlap_divu(i,j,k,rhsarr,velarr,dmskarr,dxinvarr,is_rz); + mlndlap_divu(i,j,k,rhsarr,velarr,dmskarr,dxinvarr, + nddom,lobc,hibc,is_rz); }); #else AMREX_HOST_DEVICE_PARALLEL_FOR_3D (bx, i, j, k, { - mlndlap_divu(i,j,k,rhsarr,velarr,dmskarr,dxinvarr); + mlndlap_divu(i,j,k,rhsarr,velarr,dmskarr,dxinvarr, + nddom,lobc,hibc); }); #endif } @@ -303,6 +305,7 @@ MLNodeLaplacian::compRHS (const Vector& rhs, const Vector& frhs[ilev]->setVal(0.0); const Box& cccdom = cgeom.Domain(); + const Box& nddom = amrex::surroundingNodes(fgeom.Domain()); const auto fdxinv = fgeom.InvCellSizeArray(); const iMultiFab& fdmsk = *m_dirichlet_mask[ilev+1][0]; @@ -370,7 +373,8 @@ MLNodeLaplacian::compRHS (const Vector& rhs, const Vector& } else { rarr(i,j,k) = 0.0; } - mlndlap_divu_compute_fine_contrib(i,j,k,fvbx,rarr,varr,is_rz,fdxinv); + mlndlap_divu_compute_fine_contrib(i,j,k,fvbx,rarr,varr,fdxinv, + nddom,lobc,hibc,is_rz); }); #else AMREX_HOST_DEVICE_FOR_3D(bx_rhs, i, j, k, @@ -380,7 +384,8 @@ MLNodeLaplacian::compRHS (const Vector& rhs, const Vector& } else { rarr(i,j,k) = 0.0; } - mlndlap_divu_compute_fine_contrib(i,j,k,fvbx,rarr,varr,fdxinv); + mlndlap_divu_compute_fine_contrib(i,j,k,fvbx,rarr,varr,fdxinv, + nddom,lobc,hibc); }); #endif @@ -1745,6 +1750,7 @@ MLNodeLaplacian::compSyncResidualCoarse (MultiFab& sync_resid, const MultiFab& a const BoxArray& ccba = m_grids[0][0]; const BoxArray& ndba = amrex::convert(ccba, IntVect::TheNodeVector()); const BoxArray& ccfba = amrex::convert(fine_grids, IntVect::TheZeroVector()); + const auto lobc = LoBC(); const auto hibc = HiBC(); @@ -1895,7 +1901,7 @@ MLNodeLaplacian::compSyncResidualCoarse (MultiFab& sync_resid, const MultiFab& a Array4 const& intgarr = intg->const_array(mfi); AMREX_HOST_DEVICE_FOR_3D(bx, i, j, k, { - mlndlap_divu_eb(i,j,k,rhsarr,uarr,vfracarr,intgarr,dmskarr,dxinv); + mlndlap_divu_eb(i,j,k,rhsarr,uarr,vfracarr,intgarr,dmskarr,dxinv,nddom,lobc,hibc); }); } else @@ -1904,12 +1910,12 @@ MLNodeLaplacian::compSyncResidualCoarse (MultiFab& sync_resid, const MultiFab& a #if (AMREX_SPACEDIM == 2) AMREX_HOST_DEVICE_PARALLEL_FOR_3D (bx, i, j, k, { - mlndlap_divu(i,j,k,rhsarr,uarr,dmskarr,dxinv,is_rz); + mlndlap_divu(i,j,k,rhsarr,uarr,dmskarr,dxinv,nddom,lobc,hibc,is_rz); }); #else AMREX_HOST_DEVICE_PARALLEL_FOR_3D (bx, i, j, k, { - mlndlap_divu(i,j,k,rhsarr,uarr,dmskarr,dxinv); + mlndlap_divu(i,j,k,rhsarr,uarr,dmskarr,dxinv,nddom,lobc,hibc); }); #endif } @@ -2044,6 +2050,9 @@ MLNodeLaplacian::compSyncResidualFine (MultiFab& sync_resid, const MultiFab& phi const MultiFab& sigma_orig = *m_sigma[0][0][0]; const iMultiFab& dmsk = *m_dirichlet_mask[0][0]; + const auto lobc = LoBC(); + const auto hibc = HiBC(); + #ifdef AMREX_USE_EB auto factory = dynamic_cast(m_factory[0][0].get()); const FabArray* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr; @@ -2053,6 +2062,7 @@ MLNodeLaplacian::compSyncResidualFine (MultiFab& sync_resid, const MultiFab& phi const Geometry& geom = m_geom[0][0]; const Box& ccdom = geom.Domain(); + const Box& nddom = amrex::surroundingNodes(ccdom); const auto dxinv = geom.InvCellSizeArray(); #if (AMREX_SPACEDIM == 2) bool is_rz = m_is_rz; @@ -2142,7 +2152,7 @@ MLNodeLaplacian::compSyncResidualFine (MultiFab& sync_resid, const MultiFab& phi Array4 const& intgarr = intg->const_array(mfi); AMREX_HOST_DEVICE_FOR_3D(bx, i, j, k, { - mlndlap_divu_eb(i,j,k,rhsarr,uarr,vfracarr,intgarr,tmpmaskarr,dxinv); + mlndlap_divu_eb(i,j,k,rhsarr,uarr,vfracarr,intgarr,tmpmaskarr,dxinv,nddom,lobc,hibc); }); } else @@ -2151,12 +2161,12 @@ MLNodeLaplacian::compSyncResidualFine (MultiFab& sync_resid, const MultiFab& phi #if (AMREX_SPACEDIM == 2) AMREX_HOST_DEVICE_PARALLEL_FOR_3D (bx, i, j, k, { - mlndlap_divu(i,j,k,rhsarr,uarr,tmpmaskarr,dxinv,is_rz); + mlndlap_divu(i,j,k,rhsarr,uarr,tmpmaskarr,dxinv,nddom,lobc,hibc,is_rz); }); #else AMREX_HOST_DEVICE_PARALLEL_FOR_3D (bx, i, j, k, { - mlndlap_divu(i,j,k,rhsarr,uarr,tmpmaskarr,dxinv); + mlndlap_divu(i,j,k,rhsarr,uarr,tmpmaskarr,dxinv,nddom,lobc,hibc); }); #endif }