From ed9709cb299dd2233b29a3bb569dbf717bfd1609 Mon Sep 17 00:00:00 2001 From: Ankith A Das Date: Wed, 19 Aug 2026 08:46:38 +0530 Subject: [PATCH] Index-based mlebabeclap_adotx kernel --- .../MLMG/AMReX_MLEBABecLap_2D_K.H | 209 ++++----- .../MLMG/AMReX_MLEBABecLap_3D_K.H | 429 +++++++++--------- .../MLMG/AMReX_MLEBABecLap_F.cpp | 27 +- 3 files changed, 336 insertions(+), 329 deletions(-) diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H index 1c4e6153927..8c19a69b917 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H @@ -163,17 +163,14 @@ void mlebabeclap_adotx_centroid (Box const& box, Array4 const& y, } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlebabeclap_adotx (Box const& box, Array4 const& y, +void mlebabeclap_adotx (int i, int j, int k, int n, Array4 const& y, Array4 const& x, Array4 const& a, Array4 const& bX, Array4 const& bY, - Array4 const& ccm, Array4 const& flag, - Array4 const& vfrc, Array4 const& apx, - Array4 const& apy, Array4 const& fcx, - Array4 const& fcy, Array4 const& ba, - Array4 const& bc, Array4 const& beb, + Array4 const& ccm, EBData const& ebdata, + Array4 const& beb, bool is_dirichlet, Array4 const& phieb, bool is_inhomog, GpuArray const& dxinv, - Real alpha, Real beta, int ncomp, + Real alpha, Real beta, bool beta_on_centroid, bool phi_on_centroid) noexcept { Real dhx = beta*dxinv[0]*dxinv[0]; @@ -184,120 +181,126 @@ void mlebabeclap_adotx (Box const& box, Array4 const& y, bool beta_on_center = !(beta_on_centroid); bool phi_on_center = !( phi_on_centroid); - amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept + auto const& flag = ebdata.get(); + auto const& vfrc = ebdata.get(); + auto const& apx = ebdata.get(); + auto const& apy = ebdata.get(); + auto const& fcx = ebdata.get(); + auto const& fcy = ebdata.get(); + auto const& ba = ebdata.get(); + auto const& bc = ebdata.get(); + + if (flag(i,j,k).isCovered()) { - if (flag(i,j,k).isCovered()) - { - y(i,j,k,n) = Real(0.0); - } - else if (flag(i,j,k).isRegular()) - { - y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) - - dhx * (bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i ,j,k,n)) - - bX(i ,j,k,n)*(x(i ,j,k,n) - x(i-1,j,k,n))) - - dhy * (bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j ,k,n)) - - bY(i,j ,k,n)*(x(i,j ,k,n) - x(i,j-1,k,n))); - } - else - { - Real kappa = vfrc(i,j,k); - Real apxm = apx(i,j,k); - Real apxp = apx(i+1,j,k); - Real apym = apy(i,j,k); - Real apyp = apy(i,j+1,k); + y(i,j,k,n) = Real(0.0); + } + else if (flag(i,j,k).isRegular()) + { + y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + - dhx * (bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i ,j,k,n)) + - bX(i ,j,k,n)*(x(i ,j,k,n) - x(i-1,j,k,n))) + - dhy * (bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j ,k,n)) + - bY(i,j ,k,n)*(x(i,j ,k,n) - x(i,j-1,k,n))); + } + else + { + Real kappa = vfrc(i,j,k); + Real apxm = apx(i,j,k); + Real apxp = apx(i+1,j,k); + Real apym = apy(i,j,k); + Real apyp = apy(i,j+1,k); - Real fxm = bX(i,j,k,n) * (x(i,j,k,n)-x(i-1,j,k,n)); - if (apxm != Real(0.0) && apxm != Real(1.0)) { - int jj = j + static_cast(std::copysign(Real(1.0),fcx(i,j,k))); - Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k)) : Real(0.0); - if (beta_on_center && phi_on_center) { - fxm = (Real(1.0)-fracy)*fxm + fracy*bX(i,jj,k,n)*(x(i,jj,k,n)-x(i-1,jj,k,n)); - } else if (beta_on_centroid && phi_on_center) { - fxm = bX(i,j,k,n) * ( (Real(1.0)-fracy)*(x(i, j,k,n)-x(i-1, j,k,n)) - + fracy *(x(i,jj,k,n)-x(i-1,jj,k,n)) ); - } + Real fxm = bX(i,j,k,n) * (x(i,j,k,n)-x(i-1,j,k,n)); + if (apxm != Real(0.0) && apxm != Real(1.0)) { + int jj = j + static_cast(std::copysign(Real(1.0),fcx(i,j,k))); + Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k)) : Real(0.0); + if (beta_on_center && phi_on_center) { + fxm = (Real(1.0)-fracy)*fxm + fracy*bX(i,jj,k,n)*(x(i,jj,k,n)-x(i-1,jj,k,n)); + } else if (beta_on_centroid && phi_on_center) { + fxm = bX(i,j,k,n) * ( (Real(1.0)-fracy)*(x(i, j,k,n)-x(i-1, j,k,n)) + + fracy *(x(i,jj,k,n)-x(i-1,jj,k,n)) ); } + } - Real fxp = bX(i+1,j,k,n)*(x(i+1,j,k,n)-x(i,j,k,n)); - if (apxp != Real(0.0) && apxp != Real(1.0)) { - int jj = j + static_cast(std::copysign(Real(1.0),fcx(i+1,j,k))); - Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k)) : Real(0.0); - if (beta_on_center && phi_on_center) { - fxp = (Real(1.0)-fracy)*fxp + fracy*bX(i+1,jj,k,n)*(x(i+1,jj,k,n)-x(i,jj,k,n)); - } else if (beta_on_centroid && phi_on_center) { - fxp = bX(i+1,j,k,n) * ( (Real(1.0)-fracy)*(x(i+1, j,k,n)-x(i, j,k,n)) - + fracy *(x(i+1,jj,k,n)-x(i,jj,k,n)) ); - } + Real fxp = bX(i+1,j,k,n)*(x(i+1,j,k,n)-x(i,j,k,n)); + if (apxp != Real(0.0) && apxp != Real(1.0)) { + int jj = j + static_cast(std::copysign(Real(1.0),fcx(i+1,j,k))); + Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k)) : Real(0.0); + if (beta_on_center && phi_on_center) { + fxp = (Real(1.0)-fracy)*fxp + fracy*bX(i+1,jj,k,n)*(x(i+1,jj,k,n)-x(i,jj,k,n)); + } else if (beta_on_centroid && phi_on_center) { + fxp = bX(i+1,j,k,n) * ( (Real(1.0)-fracy)*(x(i+1, j,k,n)-x(i, j,k,n)) + + fracy *(x(i+1,jj,k,n)-x(i,jj,k,n)) ); } + } - Real fym = bY(i,j,k,n)*(x(i,j,k,n)-x(i,j-1,k,n)); - if (apym != Real(0.0) && apym != Real(1.0)) { - int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j,k))); - Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k)) : Real(0.0); - if (beta_on_center && phi_on_center) { - fym = (Real(1.0)-fracx)*fym + fracx*bY(ii,j,k,n)*(x(ii,j,k,n)-x(ii,j-1,k,n)); - } else if (beta_on_centroid && phi_on_center) { - fym = bY(i,j,k,n) * ( (Real(1.0)-fracx)*(x( i,j,k,n)-x( i,j-1,k,n)) - + fracx *(x(ii,j,k,n)-x(ii,j-1,k,n)) ); - } + Real fym = bY(i,j,k,n)*(x(i,j,k,n)-x(i,j-1,k,n)); + if (apym != Real(0.0) && apym != Real(1.0)) { + int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j,k))); + Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k)) : Real(0.0); + if (beta_on_center && phi_on_center) { + fym = (Real(1.0)-fracx)*fym + fracx*bY(ii,j,k,n)*(x(ii,j,k,n)-x(ii,j-1,k,n)); + } else if (beta_on_centroid && phi_on_center) { + fym = bY(i,j,k,n) * ( (Real(1.0)-fracx)*(x( i,j,k,n)-x( i,j-1,k,n)) + + fracx *(x(ii,j,k,n)-x(ii,j-1,k,n)) ); } + } - Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n)-x(i,j,k,n)); - if (apyp != Real(0.0) && apyp != Real(1.0)) { - int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j+1,k))); - Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k)) : Real(0.0); - if (beta_on_center && phi_on_center) { - fyp = (Real(1.0)-fracx)*fyp + fracx*bY(ii,j+1,k,n)*(x(ii,j+1,k,n)-x(ii,j,k,n)); - } else if (beta_on_centroid && phi_on_center) { - fyp = bY(i,j+1,k,n) * ( (Real(1.0)-fracx)*(x( i,j+1,k,n)-x( i,j,k,n)) - + fracx *(x(ii,j+1,k,n)-x(ii,j,k,n)) ); - } + Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n)-x(i,j,k,n)); + if (apyp != Real(0.0) && apyp != Real(1.0)) { + int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j+1,k))); + Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k)) : Real(0.0); + if (beta_on_center && phi_on_center) { + fyp = (Real(1.0)-fracx)*fyp + fracx*bY(ii,j+1,k,n)*(x(ii,j+1,k,n)-x(ii,j,k,n)); + } else if (beta_on_centroid && phi_on_center) { + fyp = bY(i,j+1,k,n) * ( (Real(1.0)-fracx)*(x( i,j+1,k,n)-x( i,j,k,n)) + + fracx *(x(ii,j+1,k,n)-x(ii,j,k,n)) ); } + } - Real feb = Real(0.0); - if (is_dirichlet) { - Real dapx = (apxm-apxp)/dxinv[1]; - Real dapy = (apym-apyp)/dxinv[0]; - Real anorm = std::hypot(dapx,dapy); - Real anorminv = Real(1.0)/anorm; - Real anrmx = dapx * anorminv; - Real anrmy = dapy * anorminv; - - Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0); + Real feb = Real(0.0); + if (is_dirichlet) { + Real dapx = (apxm-apxp)/dxinv[1]; + Real dapy = (apym-apyp)/dxinv[0]; + Real anorm = std::hypot(dapx,dapy); + Real anorminv = Real(1.0)/anorm; + Real anrmx = dapx * anorminv; + Real anrmy = dapy * anorminv; - Real bctx = bc(i,j,k,0); - Real bcty = bc(i,j,k,1); - Real dx_eb = get_dx_eb(kappa); + Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0); - Real dg, gx, gy, sx, sy; - if (std::abs(anrmx) > std::abs(anrmy)) { - dg = dx_eb / std::abs(anrmx); - } else { - dg = dx_eb / std::abs(anrmy); - } - gx = (bctx - dg*anrmx); - gy = (bcty - dg*anrmy); - sx = std::copysign(Real(1.0),anrmx); - sy = std::copysign(Real(1.0),anrmy); + Real bctx = bc(i,j,k,0); + Real bcty = bc(i,j,k,1); + Real dx_eb = get_dx_eb(kappa); - int ii = i - static_cast(sx); - int jj = j - static_cast(sy); + Real dg, gx, gy, sx, sy; + if (std::abs(anrmx) > std::abs(anrmy)) { + dg = dx_eb / std::abs(anrmx); + } else { + dg = dx_eb / std::abs(anrmy); + } + gx = (bctx - dg*anrmx); + gy = (bcty - dg*anrmy); + sx = std::copysign(Real(1.0),anrmx); + sy = std::copysign(Real(1.0),anrmy); - Real phig = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy) * x(i ,j ,k,n) - + ( - gx*sx - gx*gy*sx*sy) * x(ii,j ,k,n) - + ( - gy*sy - gx*gy*sx*sy) * x(i ,jj,k,n) - + ( + gx*gy*sx*sy) * x(ii,jj,k,n) ; + int ii = i - static_cast(sx); + int jj = j - static_cast(sy); - Real dphidn = (phib-phig) / dg; + Real phig = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy) * x(i ,j ,k,n) + + ( - gx*sx - gx*gy*sx*sy) * x(ii,j ,k,n) + + ( - gy*sy - gx*gy*sx*sy) * x(i ,jj,k,n) + + ( + gx*gy*sx*sy) * x(ii,jj,k,n) ; - feb = dphidn * ba(i,j,k) * beb(i,j,k,n); - } + Real dphidn = (phib-phig) / dg; - y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) * - (dhx*(apxm*fxm-apxp*fxp) + - dhy*(apym*fym-apyp*fyp) - dh*feb); + feb = dphidn * ba(i,j,k) * beb(i,j,k,n); } - }); + + y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) * + (dhx*(apxm*fxm-apxp*fxp) + + dhy*(apym*fym-apyp*fyp) - dh*feb); + } } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H index 538437cd580..abfd8bfe8a6 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H @@ -213,19 +213,14 @@ void mlebabeclap_adotx_centroid (Box const& box, Array4 const& y, } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlebabeclap_adotx (Box const& box, Array4 const& y, +void mlebabeclap_adotx (int i, int j, int k, int n, Array4 const& y, Array4 const& x, Array4 const& a, Array4 const& bX, Array4 const& bY, Array4 const& bZ, Array4 const& ccm, - Array4 const& flag, - Array4 const& vfrc, Array4 const& apx, - Array4 const& apy, Array4 const& apz, - Array4 const& fcx, Array4 const& fcy, - Array4 const& fcz, Array4 const& ba, - Array4 const& bc, Array4 const& beb, + EBData const& ebdata, Array4 const& beb, bool is_dirichlet, Array4 const& phieb, bool is_inhomog, GpuArray const& dxinv, - Real alpha, Real beta, int ncomp, + Real alpha, Real beta, bool beta_on_centroid, bool phi_on_centroid) noexcept { Real dhx = beta*dxinv[0]*dxinv[0]; @@ -235,232 +230,240 @@ void mlebabeclap_adotx (Box const& box, Array4 const& y, bool beta_on_center = !(beta_on_centroid); bool phi_on_center = !( phi_on_centroid); - amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept + auto const& flag = ebdata.get(); + auto const& vfrc = ebdata.get(); + auto const& apx = ebdata.get(); + auto const& apy = ebdata.get(); + auto const& apz = ebdata.get(); + auto const& fcx = ebdata.get(); + auto const& fcy = ebdata.get(); + auto const& fcz = ebdata.get(); + auto const& ba = ebdata.get(); + auto const& bc = ebdata.get(); + + if (flag(i,j,k).isCovered()) { - if (flag(i,j,k).isCovered()) - { - y(i,j,k,n) = Real(0.0); - } - else if (flag(i,j,k).isRegular()) - { - y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) - - dhx * (bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i ,j,k,n)) - -bX(i ,j,k,n)*(x(i ,j,k,n) - x(i-1,j,k,n))) - - dhy * (bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j ,k,n)) - -bY(i,j ,k,n)*(x(i,j ,k,n) - x(i,j-1,k,n))) - - dhz * (bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k ,n)) - -bZ(i,j,k ,n)*(x(i,j,k ,n) - x(i,j,k-1,n))); - } - else - { - Real kappa = vfrc(i,j,k); - Real apxm = apx(i,j,k); - Real apxp = apx(i+1,j,k); - Real apym = apy(i,j,k); - Real apyp = apy(i,j+1,k); - Real apzm = apz(i,j,k); - Real apzp = apz(i,j,k+1); + y(i,j,k,n) = Real(0.0); + } + else if (flag(i,j,k).isRegular()) + { + y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + - dhx * (bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i ,j,k,n)) + -bX(i ,j,k,n)*(x(i ,j,k,n) - x(i-1,j,k,n))) + - dhy * (bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j ,k,n)) + -bY(i,j ,k,n)*(x(i,j ,k,n) - x(i,j-1,k,n))) + - dhz * (bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k ,n)) + -bZ(i,j,k ,n)*(x(i,j,k ,n) - x(i,j,k-1,n))); + } + else + { + Real kappa = vfrc(i,j,k); + Real apxm = apx(i,j,k); + Real apxp = apx(i+1,j,k); + Real apym = apy(i,j,k); + Real apyp = apy(i,j+1,k); + Real apzm = apz(i,j,k); + Real apzp = apz(i,j,k+1); - Real fxm = bX(i,j,k,n)*(x(i,j,k,n) - x(i-1,j,k,n)); - if (apxm != Real(0.0) && apxm != Real(1.0)) { - int jj = j + static_cast(std::copysign(Real(1.0), fcx(i,j,k,0))); - int kk = k + static_cast(std::copysign(Real(1.0), fcx(i,j,k,1))); - Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0); - Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0); - if (beta_on_center && phi_on_center) - { - fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm + - fracy*(Real(1.0)-fracz)*bX(i,jj,k ,n)*(x(i,jj,k ,n)-x(i-1,jj,k ,n)) + - fracz*(Real(1.0)-fracy)*bX(i,j ,kk,n)*(x(i,j ,kk,n)-x(i-1,j ,kk,n)) + - fracy* fracz *bX(i,jj,kk,n)*(x(i,jj,kk,n)-x(i-1,jj,kk,n)); - } - else if (beta_on_centroid && phi_on_center) - { - fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(x(i, j, k,n)-x(i-1, j, k,n)) + - fracy *(Real(1.0)-fracz)*(x(i,jj, k,n)-x(i-1,jj, k,n)) + - fracz *(Real(1.0)-fracy)*(x(i, j,kk,n)-x(i-1, j,kk,n)) + - fracy * fracz *(x(i,jj,kk,n)-x(i-1,jj,kk,n)); - fxm *= bX(i,j,k,n); - } + Real fxm = bX(i,j,k,n)*(x(i,j,k,n) - x(i-1,j,k,n)); + if (apxm != Real(0.0) && apxm != Real(1.0)) { + int jj = j + static_cast(std::copysign(Real(1.0), fcx(i,j,k,0))); + int kk = k + static_cast(std::copysign(Real(1.0), fcx(i,j,k,1))); + Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0); + Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0); + if (beta_on_center && phi_on_center) + { + fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm + + fracy*(Real(1.0)-fracz)*bX(i,jj,k ,n)*(x(i,jj,k ,n)-x(i-1,jj,k ,n)) + + fracz*(Real(1.0)-fracy)*bX(i,j ,kk,n)*(x(i,j ,kk,n)-x(i-1,j ,kk,n)) + + fracy* fracz *bX(i,jj,kk,n)*(x(i,jj,kk,n)-x(i-1,jj,kk,n)); } - - Real fxp = bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i,j,k,n)); - if (apxp != Real(0.0) && apxp != Real(1.0)) { - int jj = j + static_cast(std::copysign(Real(1.0),fcx(i+1,j,k,0))); - int kk = k + static_cast(std::copysign(Real(1.0),fcx(i+1,j,k,1))); - Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k,0)) : Real(0.0); - Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk)) ? std::abs(fcx(i+1,j,k,1)) : Real(0.0); - if (beta_on_center && phi_on_center) - { - fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp + - fracy*(Real(1.0)-fracz)*bX(i+1,jj,k ,n)*(x(i+1,jj,k ,n)-x(i,jj,k ,n)) + - fracz*(Real(1.0)-fracy)*bX(i+1,j ,kk,n)*(x(i+1,j ,kk,n)-x(i,j ,kk,n)) + - fracy* fracz *bX(i+1,jj,kk,n)*(x(i+1,jj,kk,n)-x(i,jj,kk,n)); - } - else if (beta_on_centroid && phi_on_center) - { - fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(x(i+1, j, k,n)-x(i, j, k,n)) + - fracy *(Real(1.0)-fracz)*(x(i+1,jj, k,n)-x(i,jj, k,n)) + - fracz *(Real(1.0)-fracy)*(x(i+1, j,kk,n)-x(i, j,kk,n)) + - fracy * fracz *(x(i+1,jj,kk,n)-x(i,jj,kk,n)); - fxp *= bX(i+1,j,k,n); - - } + else if (beta_on_centroid && phi_on_center) + { + fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(x(i, j, k,n)-x(i-1, j, k,n)) + + fracy *(Real(1.0)-fracz)*(x(i,jj, k,n)-x(i-1,jj, k,n)) + + fracz *(Real(1.0)-fracy)*(x(i, j,kk,n)-x(i-1, j,kk,n)) + + fracy * fracz *(x(i,jj,kk,n)-x(i-1,jj,kk,n)); + fxm *= bX(i,j,k,n); } + } - Real fym = bY(i,j,k,n)*(x(i,j,k,n) - x(i,j-1,k,n)); - if (apym != Real(0.0) && apym != Real(1.0)) { - int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j,k,0))); - int kk = k + static_cast(std::copysign(Real(1.0),fcy(i,j,k,1))); - Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0); - Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0); - if (beta_on_center && phi_on_center) - { - fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym + - fracx*(Real(1.0)-fracz)*bY(ii,j,k ,n)*(x(ii,j,k ,n)-x(ii,j-1,k ,n)) + - fracz*(Real(1.0)-fracx)*bY(i ,j,kk,n)*(x(i ,j,kk,n)-x(i ,j-1,kk,n)) + - fracx* fracz *bY(ii,j,kk,n)*(x(ii,j,kk,n)-x(ii,j-1,kk,n)); - } - else if (beta_on_centroid && phi_on_center) - { - fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(x( i,j, k,n)-x( i,j-1, k,n)) + - fracx *(Real(1.0)-fracz)*(x(ii,j, k,n)-x(ii,j-1, k,n)) + - fracz *(Real(1.0)-fracx)*(x(i ,j,kk,n)-x( i,j-1,kk,n)) + - fracx * fracz *(x(ii,j,kk,n)-x(ii,j-1,kk,n)); - fym *= bY(i,j,k,n); - - } + Real fxp = bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i,j,k,n)); + if (apxp != Real(0.0) && apxp != Real(1.0)) { + int jj = j + static_cast(std::copysign(Real(1.0),fcx(i+1,j,k,0))); + int kk = k + static_cast(std::copysign(Real(1.0),fcx(i+1,j,k,1))); + Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k,0)) : Real(0.0); + Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk)) ? std::abs(fcx(i+1,j,k,1)) : Real(0.0); + if (beta_on_center && phi_on_center) + { + fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp + + fracy*(Real(1.0)-fracz)*bX(i+1,jj,k ,n)*(x(i+1,jj,k ,n)-x(i,jj,k ,n)) + + fracz*(Real(1.0)-fracy)*bX(i+1,j ,kk,n)*(x(i+1,j ,kk,n)-x(i,j ,kk,n)) + + fracy* fracz *bX(i+1,jj,kk,n)*(x(i+1,jj,kk,n)-x(i,jj,kk,n)); } + else if (beta_on_centroid && phi_on_center) + { + fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(x(i+1, j, k,n)-x(i, j, k,n)) + + fracy *(Real(1.0)-fracz)*(x(i+1,jj, k,n)-x(i,jj, k,n)) + + fracz *(Real(1.0)-fracy)*(x(i+1, j,kk,n)-x(i, j,kk,n)) + + fracy * fracz *(x(i+1,jj,kk,n)-x(i,jj,kk,n)); + fxp *= bX(i+1,j,k,n); - Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j,k,n)); - if (apyp != Real(0.0) && apyp != Real(1.0)) { - int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j+1,k,0))); - int kk = k + static_cast(std::copysign(Real(1.0),fcy(i,j+1,k,1))); - Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k,0)) : Real(0.0); - Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk)) ? std::abs(fcy(i,j+1,k,1)) : Real(0.0); - if (beta_on_center && phi_on_center) - { - fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp + - fracx*(Real(1.0)-fracz)*bY(ii,j+1,k ,n)*(x(ii,j+1,k ,n)-x(ii,j,k ,n)) + - fracz*(Real(1.0)-fracx)*bY(i ,j+1,kk,n)*(x(i ,j+1,kk,n)-x(i ,j,kk,n)) + - fracx* fracz *bY(ii,j+1,kk,n)*(x(ii,j+1,kk,n)-x(ii,j,kk,n)); - } - else if (beta_on_centroid && phi_on_center) - { - fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(x( i,j+1, k,n)-x( i,j, k,n)) + - fracx *(Real(1.0)-fracz)*(x(ii,j+1, k,n)-x(ii,j, k,n)) + - fracz *(Real(1.0)-fracx)*(x( i,j+1,kk,n)-x( i,j,kk,n)) + - fracx * fracz *(x(ii,j+1,kk,n)-x(ii,j,kk,n)); - fyp *= bY(i,j+1,k,n); - - } } + } - Real fzm = bZ(i,j,k,n)*(x(i,j,k,n) - x(i,j,k-1,n)); - if (apzm != Real(0.0) && apzm != Real(1.0)) { - int ii = i + static_cast(std::copysign(Real(1.0),fcz(i,j,k,0))); - int jj = j + static_cast(std::copysign(Real(1.0),fcz(i,j,k,1))); - Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0); - Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0); - if (beta_on_center && phi_on_center) - { - fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm + - fracx*(Real(1.0)-fracy)*bZ(ii,j ,k,n)*(x(ii,j ,k,n)-x(ii,j ,k-1,n)) + - fracy*(Real(1.0)-fracx)*bZ(i ,jj,k,n)*(x(i ,jj,k,n)-x(i ,jj,k-1,n)) + - fracx* fracy *bZ(ii,jj,k,n)*(x(ii,jj,k,n)-x(ii,jj,k-1,n)); - } - else if (beta_on_centroid && phi_on_center) - { - fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(x( i, j,k,n)-x( i, j,k-1,n)) + - fracx *(Real(1.0)-fracy)*(x(ii, j,k,n)-x(ii, j,k-1,n)) + - fracy *(Real(1.0)-fracx)*(x( i,jj,k,n)-x( i,jj,k-1,n)) + - fracx * fracy *(x(ii,jj,k,n)-x(ii,jj,k-1,n)); - fzm *= bZ(i,j,k,n); - - } + Real fym = bY(i,j,k,n)*(x(i,j,k,n) - x(i,j-1,k,n)); + if (apym != Real(0.0) && apym != Real(1.0)) { + int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j,k,0))); + int kk = k + static_cast(std::copysign(Real(1.0),fcy(i,j,k,1))); + Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0); + Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0); + if (beta_on_center && phi_on_center) + { + fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym + + fracx*(Real(1.0)-fracz)*bY(ii,j,k ,n)*(x(ii,j,k ,n)-x(ii,j-1,k ,n)) + + fracz*(Real(1.0)-fracx)*bY(i ,j,kk,n)*(x(i ,j,kk,n)-x(i ,j-1,kk,n)) + + fracx* fracz *bY(ii,j,kk,n)*(x(ii,j,kk,n)-x(ii,j-1,kk,n)); } + else if (beta_on_centroid && phi_on_center) + { + fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(x( i,j, k,n)-x( i,j-1, k,n)) + + fracx *(Real(1.0)-fracz)*(x(ii,j, k,n)-x(ii,j-1, k,n)) + + fracz *(Real(1.0)-fracx)*(x(i ,j,kk,n)-x( i,j-1,kk,n)) + + fracx * fracz *(x(ii,j,kk,n)-x(ii,j-1,kk,n)); + fym *= bY(i,j,k,n); - Real fzp = bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k,n)); - if (apzp != Real(0.0) && apzp != Real(1.0)) { - int ii = i + static_cast(std::copysign(Real(1.0),fcz(i,j,k+1,0))); - int jj = j + static_cast(std::copysign(Real(1.0),fcz(i,j,k+1,1))); - Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1)) ? std::abs(fcz(i,j,k+1,0)) : Real(0.0); - Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1)) ? std::abs(fcz(i,j,k+1,1)) : Real(0.0); - if (beta_on_center && phi_on_center) - { - fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp + - fracx*(Real(1.0)-fracy)*bZ(ii,j ,k+1,n)*(x(ii,j ,k+1,n)-x(ii,j ,k,n)) + - fracy*(Real(1.0)-fracx)*bZ(i ,jj,k+1,n)*(x(i ,jj,k+1,n)-x(i ,jj,k,n)) + - fracx* fracy *bZ(ii,jj,k+1,n)*(x(ii,jj,k+1,n)-x(ii,jj,k,n)); - } - else if (beta_on_centroid && phi_on_center) - { - fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(x( i, j,k+1,n)-x( i, j,k,n)) + - fracx *(Real(1.0)-fracy)*(x(ii, j,k+1,n)-x(ii, j,k,n)) + - fracy *(Real(1.0)-fracx)*(x( i,jj,k+1,n)-x( i,jj,k,n)) + - fracx * fracy *(x(ii,jj,k+1,n)-x(ii,jj,k,n)); - fzp *= bZ(i,j,k+1,n); + } + } - } + Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j,k,n)); + if (apyp != Real(0.0) && apyp != Real(1.0)) { + int ii = i + static_cast(std::copysign(Real(1.0),fcy(i,j+1,k,0))); + int kk = k + static_cast(std::copysign(Real(1.0),fcy(i,j+1,k,1))); + Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k,0)) : Real(0.0); + Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk)) ? std::abs(fcy(i,j+1,k,1)) : Real(0.0); + if (beta_on_center && phi_on_center) + { + fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp + + fracx*(Real(1.0)-fracz)*bY(ii,j+1,k ,n)*(x(ii,j+1,k ,n)-x(ii,j,k ,n)) + + fracz*(Real(1.0)-fracx)*bY(i ,j+1,kk,n)*(x(i ,j+1,kk,n)-x(i ,j,kk,n)) + + fracx* fracz *bY(ii,j+1,kk,n)*(x(ii,j+1,kk,n)-x(ii,j,kk,n)); } + else if (beta_on_centroid && phi_on_center) + { + fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(x( i,j+1, k,n)-x( i,j, k,n)) + + fracx *(Real(1.0)-fracz)*(x(ii,j+1, k,n)-x(ii,j, k,n)) + + fracz *(Real(1.0)-fracx)*(x( i,j+1,kk,n)-x( i,j,kk,n)) + + fracx * fracz *(x(ii,j+1,kk,n)-x(ii,j,kk,n)); + fyp *= bY(i,j+1,k,n); - Real feb = Real(0.0); - if (is_dirichlet) { - Real dapx = apxm-apxp; - Real dapy = apym-apyp; - Real dapz = apzm-apzp; - Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz); - Real anorminv = Real(1.0)/anorm; - Real anrmx = dapx * anorminv; - Real anrmy = dapy * anorminv; - Real anrmz = dapz * anorminv; + } + } - Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0); + Real fzm = bZ(i,j,k,n)*(x(i,j,k,n) - x(i,j,k-1,n)); + if (apzm != Real(0.0) && apzm != Real(1.0)) { + int ii = i + static_cast(std::copysign(Real(1.0),fcz(i,j,k,0))); + int jj = j + static_cast(std::copysign(Real(1.0),fcz(i,j,k,1))); + Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0); + Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0); + if (beta_on_center && phi_on_center) + { + fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm + + fracx*(Real(1.0)-fracy)*bZ(ii,j ,k,n)*(x(ii,j ,k,n)-x(ii,j ,k-1,n)) + + fracy*(Real(1.0)-fracx)*bZ(i ,jj,k,n)*(x(i ,jj,k,n)-x(i ,jj,k-1,n)) + + fracx* fracy *bZ(ii,jj,k,n)*(x(ii,jj,k,n)-x(ii,jj,k-1,n)); + } + else if (beta_on_centroid && phi_on_center) + { + fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(x( i, j,k,n)-x( i, j,k-1,n)) + + fracx *(Real(1.0)-fracy)*(x(ii, j,k,n)-x(ii, j,k-1,n)) + + fracy *(Real(1.0)-fracx)*(x( i,jj,k,n)-x( i,jj,k-1,n)) + + fracx * fracy *(x(ii,jj,k,n)-x(ii,jj,k-1,n)); + fzm *= bZ(i,j,k,n); - Real bctx = bc(i,j,k,0); - Real bcty = bc(i,j,k,1); - Real bctz = bc(i,j,k,2); - Real dx_eb = get_dx_eb(kappa); + } + } - Real dg = dx_eb / amrex::max(std::abs(anrmx), std::abs(anrmy), - std::abs(anrmz)); - Real gx = bctx - dg*anrmx; - Real gy = bcty - dg*anrmy; - Real gz = bctz - dg*anrmz; - Real sx = std::copysign(Real(1.0),anrmx); - Real sy = std::copysign(Real(1.0),anrmy); - Real sz = std::copysign(Real(1.0),anrmz); - int ii = i - static_cast(sx); - int jj = j - static_cast(sy); - int kk = k - static_cast(sz); + Real fzp = bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k,n)); + if (apzp != Real(0.0) && apzp != Real(1.0)) { + int ii = i + static_cast(std::copysign(Real(1.0),fcz(i,j,k+1,0))); + int jj = j + static_cast(std::copysign(Real(1.0),fcz(i,j,k+1,1))); + Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1)) ? std::abs(fcz(i,j,k+1,0)) : Real(0.0); + Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1)) ? std::abs(fcz(i,j,k+1,1)) : Real(0.0); + if (beta_on_center && phi_on_center) + { + fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp + + fracx*(Real(1.0)-fracy)*bZ(ii,j ,k+1,n)*(x(ii,j ,k+1,n)-x(ii,j ,k,n)) + + fracy*(Real(1.0)-fracx)*bZ(i ,jj,k+1,n)*(x(i ,jj,k+1,n)-x(i ,jj,k,n)) + + fracx* fracy *bZ(ii,jj,k+1,n)*(x(ii,jj,k+1,n)-x(ii,jj,k,n)); + } + else if (beta_on_centroid && phi_on_center) + { + fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(x( i, j,k+1,n)-x( i, j,k,n)) + + fracx *(Real(1.0)-fracy)*(x(ii, j,k+1,n)-x(ii, j,k,n)) + + fracy *(Real(1.0)-fracx)*(x( i,jj,k+1,n)-x( i,jj,k,n)) + + fracx * fracy *(x(ii,jj,k+1,n)-x(ii,jj,k,n)); + fzp *= bZ(i,j,k+1,n); - gx = sx*gx; - gy = sy*gy; - gz = sz*gz; - Real gxy = gx*gy; - Real gxz = gx*gz; - Real gyz = gy*gz; - Real gxyz = gx*gy*gz; - Real phig = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz) * x(i ,j ,k ,n) - + (-gz - gxz - gyz - gxyz) * x(i ,j ,kk,n) - + (-gy - gxy - gyz - gxyz) * x(i ,jj,k ,n) - + (gyz + gxyz) * x(i ,jj,kk,n) - + (-gx - gxy - gxz - gxyz) * x(ii,j ,k ,n) - + (gxz + gxyz) * x(ii,j ,kk,n) - + (gxy + gxyz) * x(ii,jj,k ,n) - + (-gxyz) * x(ii,jj,kk,n); - - Real dphidn = (phib-phig)/dg; - - feb = dphidn * ba(i,j,k) * beb(i,j,k,n); } + } - y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) * - (dhx*(apxm*fxm - apxp*fxp) + - dhy*(apym*fym - apyp*fyp) + - dhz*(apzm*fzm - apzp*fzp) - dhx*feb); + Real feb = Real(0.0); + if (is_dirichlet) { + Real dapx = apxm-apxp; + Real dapy = apym-apyp; + Real dapz = apzm-apzp; + Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz); + Real anorminv = Real(1.0)/anorm; + Real anrmx = dapx * anorminv; + Real anrmy = dapy * anorminv; + Real anrmz = dapz * anorminv; + + Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0); + + Real bctx = bc(i,j,k,0); + Real bcty = bc(i,j,k,1); + Real bctz = bc(i,j,k,2); + Real dx_eb = get_dx_eb(kappa); + + Real dg = dx_eb / amrex::max(std::abs(anrmx), std::abs(anrmy), + std::abs(anrmz)); + Real gx = bctx - dg*anrmx; + Real gy = bcty - dg*anrmy; + Real gz = bctz - dg*anrmz; + Real sx = std::copysign(Real(1.0),anrmx); + Real sy = std::copysign(Real(1.0),anrmy); + Real sz = std::copysign(Real(1.0),anrmz); + int ii = i - static_cast(sx); + int jj = j - static_cast(sy); + int kk = k - static_cast(sz); + + gx = sx*gx; + gy = sy*gy; + gz = sz*gz; + Real gxy = gx*gy; + Real gxz = gx*gz; + Real gyz = gy*gz; + Real gxyz = gx*gy*gz; + Real phig = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz) * x(i ,j ,k ,n) + + (-gz - gxz - gyz - gxyz) * x(i ,j ,kk,n) + + (-gy - gxy - gyz - gxyz) * x(i ,jj,k ,n) + + (gyz + gxyz) * x(i ,jj,kk,n) + + (-gx - gxy - gxz - gxyz) * x(ii,j ,k ,n) + + (gxz + gxyz) * x(ii,j ,kk,n) + + (gxy + gxyz) * x(ii,jj,k ,n) + + (-gxyz) * x(ii,jj,kk,n); + + Real dphidn = (phib-phig)/dg; + + feb = dphidn * ba(i,j,k) * beb(i,j,k,n); } - }); + + y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) * + (dhx*(apxm*fxm - apxp*fxp) + + dhy*(apym*fym - apyp*fyp) + + dhz*(apzm*fzm - apzp*fzp) - dhx*feb); + } } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp index a8e80c954b7..85253075874 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp @@ -130,7 +130,9 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c amrex::ignore_unused(AMREX_D_DECL(domlo_x, domlo_y, domlo_z), AMREX_D_DECL(domhi_x, domhi_y, domhi_z), AMREX_D_DECL(extdir_x, extdir_y, extdir_z)); - amrex::ignore_unused(ccfab); + amrex::ignore_unused(ccfab, flagfab, vfracfab, bafab, bcfab, + AMREX_D_DECL(apxfab,apyfab,apzfab), + AMREX_D_DECL(fcxfab,fcyfab,fczfab)); #else AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bx, tbx, { @@ -147,18 +149,17 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c }); #endif } else { - AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bx, tbx, - { - mlebabeclap_adotx(tbx, yfab, xfab, afab, AMREX_D_DECL(bxfab,byfab,bzfab), - ccmfab, flagfab, vfracfab, - AMREX_D_DECL(apxfab,apyfab,apzfab), - AMREX_D_DECL(fcxfab,fcyfab,fczfab), - bafab, bcfab, bebfab, - is_eb_dirichlet, - phiebfab, - is_eb_inhomog, dxinvarr, - ascalar, bscalar, ncomp, beta_on_centroid, phi_on_centroid); - }); + auto const& ebdata = factory->getEBData(mfi); + AMREX_HOST_DEVICE_PARALLEL_FOR_4D( bx, ncomp, i, j, k, n, + { + mlebabeclap_adotx(i, j, k, n, yfab, xfab, afab, + AMREX_D_DECL(bxfab,byfab,bzfab), + ccmfab, ebdata, bebfab, + is_eb_dirichlet, + phiebfab, + is_eb_inhomog, dxinvarr, + ascalar, bscalar, beta_on_centroid, phi_on_centroid); + }); } if (has_overset) { Array4 const& osm = m_overset_mask[amrlev][mglev]->const_array(mfi);