diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H index 1c4e6153927..b158efb7efe 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H @@ -8,172 +8,172 @@ namespace amrex { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlebabeclap_adotx_centroid (Box const& box, Array4 const& y, +void mlebabeclap_adotx_centroid (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& flag, - Array4 const& vfrc, - Array4 const& apx, Array4 const& apy, - Array4 const& fcx, Array4 const& fcy, - Array4 const& ccent, Array4 const& ba, - Array4 const& bcent, Array4 const& beb, + EBData const& ebdata, + Array4 const& beb, Array4 const& phieb, const int& domlo_x, const int& domlo_y, const int& domhi_x, const int& domhi_y, const bool& on_x_face, const bool& on_y_face, bool is_eb_dirichlet, bool is_eb_inhomog, GpuArray const& dxinv, - Real alpha, Real beta, int ncomp) noexcept + Real alpha, Real beta) noexcept { Real dhx = beta*dxinv[0]*dxinv[0]; Real dhy = beta*dxinv[1]*dxinv[1]; Real dh = beta*dxinv[0]*dxinv[1]; - 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& ccent = ebdata.get(); + auto const& ba = ebdata.get(); + auto const& bcent = 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() && - ((flag(i-1,j ,k).isRegular() && flag(i+1,j ,k).isRegular() && - flag(i ,j-1,k).isRegular() && flag(i ,j+1,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))); + y(i,j,k,n) = Real(0.0); + } + else if (flag(i,j,k).isRegular() && + ((flag(i-1,j ,k).isRegular() && flag(i+1,j ,k).isRegular() && + flag(i ,j-1,k).isRegular() && flag(i ,j+1,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); + } + 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); - // First get EB-aware slope that doesn't know about extdir - bool needs_bdry_stencil = (i <= domlo_x) || (i >= domhi_x) || - (j <= domlo_y) || (j >= domhi_y); + // First get EB-aware slope that doesn't know about extdir + bool needs_bdry_stencil = (i <= domlo_x) || (i >= domhi_x) || + (j <= domlo_y) || (j >= domhi_y); - // if phi_on_centroid -- A second order least squares fit is used - // to approximate the slope on the high and low faces. Note that if - // any of the three cells --e.g., (i-1,j), (i,j), or (i-1,j)-- are - // cut, then the least squares fit is needed. This is a bit more than - // is actually needed for most cases but it will return the correct - // value in all cases. + // if phi_on_centroid -- A second order least squares fit is used + // to approximate the slope on the high and low faces. Note that if + // any of the three cells --e.g., (i-1,j), (i,j), or (i-1,j)-- are + // cut, then the least squares fit is needed. This is a bit more than + // is actually needed for most cases but it will return the correct + // value in all cases. - Real fxm = bX(i,j,k,n) * (x(i,j,k,n)-x(i-1,j,k,n)); - if ( (apxm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0)) ) - { - Real yloc_on_xface = fcx(i,j,k); + Real fxm = bX(i,j,k,n) * (x(i,j,k,n)-x(i-1,j,k,n)); + if ( (apxm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0)) ) + { + Real yloc_on_xface = fcx(i,j,k); - if(needs_bdry_stencil) { + if(needs_bdry_stencil) { - fxm = grad_x_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - yloc_on_xface,is_eb_dirichlet,is_eb_inhomog, - on_x_face, domlo_x, domhi_x, - on_y_face, domlo_y, domhi_y); + fxm = grad_x_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + yloc_on_xface,is_eb_dirichlet,is_eb_inhomog, + on_x_face, domlo_x, domhi_x, + on_y_face, domlo_y, domhi_y); - } else { - fxm = grad_x_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, - yloc_on_xface,is_eb_dirichlet,is_eb_inhomog); - } - fxm *= bX(i,j,k,n); + } else { + fxm = grad_x_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, + yloc_on_xface,is_eb_dirichlet,is_eb_inhomog); } + fxm *= bX(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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0)) ) { - Real yloc_on_xface = fcx(i+1,j,k,0); - if(needs_bdry_stencil) { - fxp = grad_x_of_phi_on_centroids_extdir(i+1,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - yloc_on_xface,is_eb_dirichlet,is_eb_inhomog, - on_x_face, domlo_x, domhi_x, - on_y_face, domlo_y, domhi_y); - - } else { - fxp = grad_x_of_phi_on_centroids(i+1,j,k,n,x,phieb,flag,ccent,bcent, - yloc_on_xface,is_eb_dirichlet,is_eb_inhomog); - } - fxp *= bX(i+1,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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0)) ) { + Real yloc_on_xface = fcx(i+1,j,k,0); + if(needs_bdry_stencil) { + fxp = grad_x_of_phi_on_centroids_extdir(i+1,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + yloc_on_xface,is_eb_dirichlet,is_eb_inhomog, + on_x_face, domlo_x, domhi_x, + on_y_face, domlo_y, domhi_y); + } else { + fxp = grad_x_of_phi_on_centroids(i+1,j,k,n,x,phieb,flag,ccent,bcent, + yloc_on_xface,is_eb_dirichlet,is_eb_inhomog); } + fxp *= bX(i+1,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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0)) ) { - Real xloc_on_yface = fcy(i,j,k,0); + } - if(needs_bdry_stencil) { + Real fym = bY(i,j,k,n)*(x(i,j,k,n)-x(i,j-1,k,n)); + if ( (apym != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0)) ) { + Real xloc_on_yface = fcy(i,j,k,0); - fym = grad_y_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - xloc_on_yface,is_eb_dirichlet,is_eb_inhomog, - on_x_face, domlo_x, domhi_x, - on_y_face, domlo_y, domhi_y); + if(needs_bdry_stencil) { - } else { - fym = grad_y_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, - xloc_on_yface,is_eb_dirichlet,is_eb_inhomog); - } - fym *= bY(i,j,k,n); + fym = grad_y_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + xloc_on_yface,is_eb_dirichlet,is_eb_inhomog, + on_x_face, domlo_x, domhi_x, + on_y_face, domlo_y, domhi_y); + + } else { + fym = grad_y_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, + xloc_on_yface,is_eb_dirichlet,is_eb_inhomog); } + fym *= bY(i,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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0)) ) { - Real xloc_on_yface = fcy(i,j+1,k,0); - if(needs_bdry_stencil) { - fyp = grad_y_of_phi_on_centroids_extdir(i,j+1,k,n,x,phieb,flag,ccent,bcent,vfrc, - xloc_on_yface,is_eb_dirichlet,is_eb_inhomog, - on_x_face, domlo_x, domhi_x, - on_y_face, domlo_y, domhi_y); + Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n)-x(i,j,k,n)); + if ( (apyp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0)) ) { + Real xloc_on_yface = fcy(i,j+1,k,0); + if(needs_bdry_stencil) { + fyp = grad_y_of_phi_on_centroids_extdir(i,j+1,k,n,x,phieb,flag,ccent,bcent,vfrc, + xloc_on_yface,is_eb_dirichlet,is_eb_inhomog, + on_x_face, domlo_x, domhi_x, + on_y_face, domlo_y, domhi_y); - } else { - fyp = grad_y_of_phi_on_centroids(i,j+1,k,n,x,phieb,flag,ccent,bcent, - xloc_on_yface,is_eb_dirichlet,is_eb_inhomog); - } - fyp *= bY(i,j+1,k,n); + } else { + fyp = grad_y_of_phi_on_centroids(i,j+1,k,n,x,phieb,flag,ccent,bcent, + xloc_on_yface,is_eb_dirichlet,is_eb_inhomog); } + fyp *= bY(i,j+1,k,n); + } - Real feb = Real(0.0); - if (is_eb_dirichlet && flag(i,j,k).isSingleValued()) - { - 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 feb = Real(0.0); + if (is_eb_dirichlet && flag(i,j,k).isSingleValued()) + { + 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; - feb = grad_eb_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - anrmx,anrmy,is_eb_inhomog, - on_x_face, domlo_x, domhi_x, - on_y_face, domlo_y, domhi_y); + feb = grad_eb_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + anrmx,anrmy,is_eb_inhomog, + on_x_face, domlo_x, domhi_x, + on_y_face, domlo_y, domhi_y); - feb *= ba(i,j,k) * beb(i,j,k,n); - } + feb *= 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); - } - }); + 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 -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 +184,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..55ae6bf50ec 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H @@ -8,224 +8,223 @@ namespace amrex { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlebabeclap_adotx_centroid (Box const& box, Array4 const& y, +void mlebabeclap_adotx_centroid (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& flag, - Array4 const& vfrc, Array4 const& apx, - Array4 const& apy, Array4 const& apz, - Array4 const& fcx, Array4 const& fcy, - Array4 const& fcz, - Array4 const& ccent, Array4 const& ba, - Array4 const& bcent, Array4 const& beb, + EBData const& ebdata, + Array4 const& beb, Array4 const& phieb, const int& domlo_x, const int& domlo_y, const int& domlo_z, const int& domhi_x, const int& domhi_y, const int& domhi_z, const bool& on_x_face, const bool& on_y_face, const bool& on_z_face, bool is_eb_dirichlet, bool is_eb_inhomog, GpuArray const& dxinv, - Real alpha, Real beta, int ncomp) noexcept + Real alpha, Real beta) noexcept { Real dhx = beta*dxinv[0]*dxinv[0]; Real dhy = beta*dxinv[1]*dxinv[1]; Real dhz = beta*dxinv[2]*dxinv[2]; - 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& ccent = ebdata.get(); + auto const& ba = ebdata.get(); + auto const& bcent = 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() && - ((flag(i-1,j ,k ).isRegular() && flag(i+1,j ,k ).isRegular() && - flag(i ,j-1,k ).isRegular() && flag(i ,j+1,k ).isRegular() && - flag(i ,j ,k-1).isRegular() && flag(i ,j ,k+1).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); - - // First get EB-aware slope that doesn't know about extdir - bool needs_bdry_stencil = (i <= domlo_x) || (i >= domhi_x) || - (j <= domlo_y) || (j >= domhi_y) || - (k <= domlo_z) || (k >= domhi_z); - - Real fxm = bX(i,j,k,n)*(x(i,j,k,n) - x(i-1,j,k,n)); - if ( (apxm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0)) ) { - Real yloc_on_xface = fcx(i,j,k,0); - Real zloc_on_xface = fcx(i,j,k,1); - - if(needs_bdry_stencil) { - fxm = grad_x_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - yloc_on_xface,zloc_on_xface, - is_eb_dirichlet,is_eb_inhomog, - on_x_face,domlo_x,domhi_x, - on_y_face,domlo_y,domhi_y, - on_z_face,domlo_z,domhi_z); - } else { - fxm = grad_x_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, - yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog); - } + y(i,j,k,n) = Real(0.0); + } + else if (flag(i,j,k).isRegular() && + ((flag(i-1,j ,k ).isRegular() && flag(i+1,j ,k ).isRegular() && + flag(i ,j-1,k ).isRegular() && flag(i ,j+1,k ).isRegular() && + flag(i ,j ,k-1).isRegular() && flag(i ,j ,k+1).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); - fxm *= bX(i,j,k,n); + // First get EB-aware slope that doesn't know about extdir + bool needs_bdry_stencil = (i <= domlo_x) || (i >= domhi_x) || + (j <= domlo_y) || (j >= domhi_y) || + (k <= domlo_z) || (k >= domhi_z); + + Real fxm = bX(i,j,k,n)*(x(i,j,k,n) - x(i-1,j,k,n)); + if ( (apxm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0)) ) { + Real yloc_on_xface = fcx(i,j,k,0); + Real zloc_on_xface = fcx(i,j,k,1); + + if(needs_bdry_stencil) { + fxm = grad_x_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + yloc_on_xface,zloc_on_xface, + is_eb_dirichlet,is_eb_inhomog, + on_x_face,domlo_x,domhi_x, + on_y_face,domlo_y,domhi_y, + on_z_face,domlo_z,domhi_z); + } else { + fxm = grad_x_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, + yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog); } - Real fxp = bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i,j,k,n)); - if ( (apxp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0)) ) { - Real yloc_on_xface = fcx(i+1,j,k,0); - Real zloc_on_xface = fcx(i+1,j,k,1); - - if(needs_bdry_stencil) { - fxp = grad_x_of_phi_on_centroids_extdir(i+1,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - yloc_on_xface,zloc_on_xface, - is_eb_dirichlet,is_eb_inhomog, - on_x_face,domlo_x,domhi_x, - on_y_face,domlo_y,domhi_y, - on_z_face,domlo_z,domhi_z); - } else { - fxp = grad_x_of_phi_on_centroids(i+1,j,k,n,x,phieb,flag,ccent,bcent, - yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog); - } + fxm *= bX(i,j,k,n); + } - fxp *= bX(i+1,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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0)) ) { + Real yloc_on_xface = fcx(i+1,j,k,0); + Real zloc_on_xface = fcx(i+1,j,k,1); + + if(needs_bdry_stencil) { + fxp = grad_x_of_phi_on_centroids_extdir(i+1,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + yloc_on_xface,zloc_on_xface, + is_eb_dirichlet,is_eb_inhomog, + on_x_face,domlo_x,domhi_x, + on_y_face,domlo_y,domhi_y, + on_z_face,domlo_z,domhi_z); + } else { + fxp = grad_x_of_phi_on_centroids(i+1,j,k,n,x,phieb,flag,ccent,bcent, + yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog); } - Real fym = bY(i,j,k,n)*(x(i,j,k,n) - x(i,j-1,k,n)); - if ( (apym != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0)) ) { - Real xloc_on_yface = fcy(i,j,k,0); - Real zloc_on_yface = fcy(i,j,k,1); - - if(needs_bdry_stencil) { - fym = grad_y_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - xloc_on_yface,zloc_on_yface, - is_eb_dirichlet,is_eb_inhomog, - on_x_face,domlo_x,domhi_x, - on_y_face,domlo_y,domhi_y, - on_z_face,domlo_z,domhi_z); - } else { - fym = grad_y_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, - xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog); - } + fxp *= bX(i+1,j,k,n); + } - fym *= bY(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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0)) ) { + Real xloc_on_yface = fcy(i,j,k,0); + Real zloc_on_yface = fcy(i,j,k,1); + + if(needs_bdry_stencil) { + fym = grad_y_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + xloc_on_yface,zloc_on_yface, + is_eb_dirichlet,is_eb_inhomog, + on_x_face,domlo_x,domhi_x, + on_y_face,domlo_y,domhi_y, + on_z_face,domlo_z,domhi_z); + } else { + fym = grad_y_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, + xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog); } - Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j,k,n)); - if ( (apyp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0)) ) { - Real xloc_on_yface = fcy(i,j+1,k,0); - Real zloc_on_yface = fcy(i,j+1,k,1); - - if(needs_bdry_stencil) { - fyp = grad_y_of_phi_on_centroids_extdir(i,j+1,k,n,x,phieb,flag,ccent,bcent,vfrc, - xloc_on_yface,zloc_on_yface, - is_eb_dirichlet,is_eb_inhomog, - on_x_face,domlo_x,domhi_x, - on_y_face,domlo_y,domhi_y, - on_z_face,domlo_z,domhi_z); - } else { - fyp = grad_y_of_phi_on_centroids(i,j+1,k,n,x,phieb,flag,ccent,bcent, - xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog); - } + fym *= bY(i,j,k,n); + } - fyp *= bY(i,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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0)) ) { + Real xloc_on_yface = fcy(i,j+1,k,0); + Real zloc_on_yface = fcy(i,j+1,k,1); + + if(needs_bdry_stencil) { + fyp = grad_y_of_phi_on_centroids_extdir(i,j+1,k,n,x,phieb,flag,ccent,bcent,vfrc, + xloc_on_yface,zloc_on_yface, + is_eb_dirichlet,is_eb_inhomog, + on_x_face,domlo_x,domhi_x, + on_y_face,domlo_y,domhi_y, + on_z_face,domlo_z,domhi_z); + } else { + fyp = grad_y_of_phi_on_centroids(i,j+1,k,n,x,phieb,flag,ccent,bcent, + xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog); } - Real fzm = bZ(i,j,k,n)*(x(i,j,k,n) - x(i,j,k-1,n)); - if ( (apzm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0)) ) { - Real xloc_on_zface = fcz(i,j,k,0); - Real yloc_on_zface = fcz(i,j,k,1); - - if(needs_bdry_stencil) { - fzm = grad_z_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - xloc_on_zface,yloc_on_zface, - is_eb_dirichlet,is_eb_inhomog, - on_x_face,domlo_x,domhi_x, - on_y_face,domlo_y,domhi_y, - on_z_face,domlo_z,domhi_z); - } else { - fzm = grad_z_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, - xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog); - } + fyp *= bY(i,j+1,k,n); + } - fzm *= bZ(i,j,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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0)) ) { + Real xloc_on_zface = fcz(i,j,k,0); + Real yloc_on_zface = fcz(i,j,k,1); + + if(needs_bdry_stencil) { + fzm = grad_z_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + xloc_on_zface,yloc_on_zface, + is_eb_dirichlet,is_eb_inhomog, + on_x_face,domlo_x,domhi_x, + on_y_face,domlo_y,domhi_y, + on_z_face,domlo_z,domhi_z); + } else { + fzm = grad_z_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent, + xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog); } - Real fzp = bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k,n)); - if ( (apzp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0)) ) { - Real xloc_on_zface = fcz(i,j,k+1,0); - Real yloc_on_zface = fcz(i,j,k+1,1); - - if(needs_bdry_stencil) { - fzp = grad_z_of_phi_on_centroids_extdir(i,j,k+1,n,x,phieb,flag,ccent,bcent,vfrc, - xloc_on_zface,yloc_on_zface, - is_eb_dirichlet,is_eb_inhomog, - on_x_face,domlo_x,domhi_x, - on_y_face,domlo_y,domhi_y, - on_z_face,domlo_z,domhi_z); - } else { - fzp = grad_z_of_phi_on_centroids(i,j,k+1,n,x,phieb,flag,ccent,bcent, - xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog); - } + fzm *= bZ(i,j,k,n); + } - fzp *= bZ(i,j,k+1,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)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0)) ) { + Real xloc_on_zface = fcz(i,j,k+1,0); + Real yloc_on_zface = fcz(i,j,k+1,1); + + if(needs_bdry_stencil) { + fzp = grad_z_of_phi_on_centroids_extdir(i,j,k+1,n,x,phieb,flag,ccent,bcent,vfrc, + xloc_on_zface,yloc_on_zface, + is_eb_dirichlet,is_eb_inhomog, + on_x_face,domlo_x,domhi_x, + on_y_face,domlo_y,domhi_y, + on_z_face,domlo_z,domhi_z); + } else { + fzp = grad_z_of_phi_on_centroids(i,j,k+1,n,x,phieb,flag,ccent,bcent, + xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog); } - Real feb = Real(0.0); - if (is_eb_dirichlet && flag(i,j,k).isSingleValued()) { - 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; - - feb = grad_eb_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, - anrmx,anrmy,anrmz,is_eb_inhomog, - on_x_face,domlo_x,domhi_x, - on_y_face,domlo_y,domhi_y, - on_z_face,domlo_z,domhi_z); - feb *= ba(i,j,k) * beb(i,j,k,n); - } + fzp *= bZ(i,j,k+1,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_eb_dirichlet && flag(i,j,k).isSingleValued()) { + 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; + + feb = grad_eb_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc, + anrmx,anrmy,anrmz,is_eb_inhomog, + on_x_face,domlo_x,domhi_x, + on_y_face,domlo_y,domhi_y, + on_z_face,domlo_z,domhi_z); + feb *= 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 -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 +234,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..5d67f3c72f0 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp @@ -23,14 +23,6 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c const auto *factory = dynamic_cast(m_factory[amrlev][mglev].get()); const FabArray* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr; - const MultiFab* vfrac = (factory) ? &(factory->getVolFrac()) : nullptr; - auto area = (factory) ? factory->getAreaFrac() - : Array{AMREX_D_DECL(nullptr,nullptr,nullptr)}; - auto fcent = (factory) ? factory->getFaceCent() - : Array{AMREX_D_DECL(nullptr,nullptr,nullptr)}; - const MultiCutFab* barea = (factory) ? &(factory->getBndryArea()) : nullptr; - const MultiCutFab* bcent = (factory) ? &(factory->getBndryCent()) : nullptr; - const auto *const ccent = (factory) ? &(factory->getCentroid()) : nullptr; const bool is_eb_dirichlet = isEBDirichlet(); const bool is_eb_inhomog = m_is_eb_inhomog && (!this->m_precond_mode); @@ -99,17 +91,7 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c } } else { Array4 const& ccmfab = ccmask.const_array(mfi); - Array4 const& flagfab = flags->const_array(mfi); - Array4 const& vfracfab = vfrac->const_array(mfi); - AMREX_D_TERM(Array4 const& apxfab = area[0]->const_array(mfi);, - Array4 const& apyfab = area[1]->const_array(mfi);, - Array4 const& apzfab = area[2]->const_array(mfi);); - AMREX_D_TERM(Array4 const& fcxfab = fcent[0]->const_array(mfi);, - Array4 const& fcyfab = fcent[1]->const_array(mfi);, - Array4 const& fczfab = fcent[2]->const_array(mfi);); - Array4 const& bafab = barea->const_array(mfi); - Array4 const& bcfab = bcent->const_array(mfi); - Array4 const& ccfab = ccent->const_array(mfi); + auto const& ebdata = factory->getEBData(mfi); Array4 const& bebfab = (is_eb_dirichlet) ? m_eb_b_coeffs[amrlev][mglev]->const_array(mfi) : foo; Array4 const& phiebfab = (is_eb_dirichlet && is_eb_inhomog) @@ -130,35 +112,30 @@ 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); #else - AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bx, tbx, - { - mlebabeclap_adotx_centroid(tbx, yfab, xfab, afab, AMREX_D_DECL(bxfab,byfab,bzfab), - flagfab, vfracfab, - AMREX_D_DECL(apxfab,apyfab,apzfab), - AMREX_D_DECL(fcxfab,fcyfab,fczfab), - ccfab, bafab, bcfab, bebfab, phiebfab, - 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), - is_eb_dirichlet, is_eb_inhomog, dxinvarr, - ascalar, bscalar, ncomp); - }); + AMREX_HOST_DEVICE_PARALLEL_FOR_4D( bx, ncomp, i, j, k, n, + { + mlebabeclap_adotx_centroid(i, j, k, n, yfab, xfab, afab, + AMREX_D_DECL(bxfab,byfab,bzfab), + ebdata, bebfab, phiebfab, + 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), + is_eb_dirichlet, is_eb_inhomog, dxinvarr, + ascalar, bscalar); + }); #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); - }); + 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);