diff --git a/Src/EB/AMReX_EBData.H b/Src/EB/AMReX_EBData.H new file mode 100644 index 0000000000..ebfe03467d --- /dev/null +++ b/Src/EB/AMReX_EBData.H @@ -0,0 +1,56 @@ +#ifndef AMREX_EB_DATA_H_ +#define AMREX_EB_DATA_H_ +#include + +#include + +namespace amrex +{ + +enum struct EBData_t : int +{ + levelset, // level set + volfrac, // volume fraction + centroid, // volume centroid + bndrycent, // boundary centroid + bndrynorm, // boundary normal + bndryarea, // boundary area + AMREX_D_DECL(apx, apy, apz), // area fraction + AMREX_D_DECL(fcx, fcy, fcz), // face centroid + AMREX_D_DECL(ecx, ecy, ecz), // edge centroid + cellflag // EBCellFlag +}; + +struct EBData +{ + template + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + auto get (int i, int j, int k) const noexcept + { + if constexpr (T == EBData_t::cellflag) { + return (*m_cell_flag)(i,j,k); + } else { + return m_real_data[static_cast(T)](i,j,k); + } + } + + template = 0> + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + auto get (int i, int j, int k, int n) const noexcept + { + return m_real_data[static_cast(T)](i,j,k,n); + } + + static constexpr int real_data_size = static_cast(EBData_t::cellflag); + + Array4 const* m_cell_flag = nullptr; + Array4 const* m_real_data = nullptr; +}; + +} +#endif diff --git a/Src/EB/AMReX_EBDataCollection.H b/Src/EB/AMReX_EBDataCollection.H index 6a8556dc9a..bf2d8f1549 100644 --- a/Src/EB/AMReX_EBDataCollection.H +++ b/Src/EB/AMReX_EBDataCollection.H @@ -9,6 +9,7 @@ namespace amrex { +class EBFArrayBoxFactory; template class FabArray; class MultiFab; class iMultiFab; @@ -19,6 +20,8 @@ class EBDataCollection { public: + friend class EBFArrayBoxFactory; + EBDataCollection (const EB2::Level& a_level, const Geometry& a_geom, const BoxArray& a_ba, const DistributionMapping& a_dm, Vector a_ngrow, EBSupport a_support); diff --git a/Src/EB/AMReX_EBFabFactory.H b/Src/EB/AMReX_EBFabFactory.H index c9db5f3947..7d20e84354 100644 --- a/Src/EB/AMReX_EBFabFactory.H +++ b/Src/EB/AMReX_EBFabFactory.H @@ -4,10 +4,12 @@ #include +#include #include #include #include #include +#include namespace amrex { @@ -89,12 +91,15 @@ public: //! One should use getMultiEBCellFlagFab for normal levels. [[nodiscard]] iMultiFab const* getCutCellMask () const noexcept { return m_ebdc->getCutCellMask(); } + [[nodiscard]] EBData getEBData (MFIter const& mfi) const noexcept; + private: EBSupport m_support; Geometry m_geom; std::shared_ptr m_ebdc; EB2::Level const* m_parent = nullptr; + Gpu::DeviceVector> m_eb_data; }; std::unique_ptr diff --git a/Src/EB/AMReX_EBFabFactory.cpp b/Src/EB/AMReX_EBFabFactory.cpp index f5154578d6..4d4f9dcee7 100644 --- a/Src/EB/AMReX_EBFabFactory.cpp +++ b/Src/EB/AMReX_EBFabFactory.cpp @@ -20,7 +20,70 @@ EBFArrayBoxFactory::EBFArrayBoxFactory (const EB2::Level& a_level, m_geom(a_geom), m_ebdc(std::make_shared(a_level,a_geom,a_ba,a_dm,a_ngrow,a_support)), m_parent(&a_level) -{} +{ + auto const& ebflags = getMultiEBCellFlagFab(); +#ifdef AMREX_USE_GPU + m_eb_data.resize(EBData::real_data_size*ebflags.local_size()); + Gpu::PinnedVector> eb_data_hv; +#else + auto& eb_data_hv = m_eb_data; +#endif + + eb_data_hv.reserve(EBData::real_data_size*ebflags.local_size()); + + for (MFIter mfi(ebflags,MFItInfo{}.DisableDeviceSync()); mfi.isValid(); ++mfi) { + Array4 a{}; + + bool cutfab_is_ok = ebflags[mfi].getType() == FabType::singlevalued; + + a = ( m_ebdc->m_levelset ) + ? m_ebdc->m_levelset->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + + a = ( m_ebdc->m_volfrac ) + ? m_ebdc->m_volfrac->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + + a = ( m_ebdc->m_centroid && cutfab_is_ok ) + ? m_ebdc->m_centroid->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + + a = ( m_ebdc->m_bndrycent && cutfab_is_ok ) + ? m_ebdc->m_bndrycent->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + + a = ( m_ebdc->m_bndrynorm && cutfab_is_ok ) + ? m_ebdc->m_bndrynorm->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + + a = ( m_ebdc->m_bndryarea && cutfab_is_ok ) + ? m_ebdc->m_bndryarea->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + a = ( m_ebdc->m_areafrac[idim] && cutfab_is_ok ) + ? m_ebdc->m_areafrac[idim]->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + } + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + a = ( m_ebdc->m_facecent[idim] && cutfab_is_ok ) + ? m_ebdc->m_facecent[idim]->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + } + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + a = ( m_ebdc->m_edgecent[idim] && cutfab_is_ok ) + ? m_ebdc->m_edgecent[idim]->const_array(mfi) : Array4{}; + eb_data_hv.push_back(a); + } + } + +#ifdef AMREX_USE_GPU + Gpu::copyAsync(Gpu::hostToDevice, eb_data_hv.begin(), eb_data_hv.end(), m_eb_data.begin()); + Gpu::streamSynchronize(); +#endif +} AMREX_NODISCARD FArrayBox* @@ -115,6 +178,19 @@ EBFArrayBoxFactory::hasEBInfo () const noexcept return m_parent->hasEBInfo(); } +EBData +EBFArrayBoxFactory::getEBData (MFIter const& mfi) const noexcept +{ + int const li = mfi.LocalIndex(); + auto const& ebflags_ma = this->getMultiEBCellFlagFab().const_arrays(); +#ifdef AMREX_USE_GPU + auto const* pebflag = ebflags_ma.dp + li; +#else + auto const* pebflag = ebflags_ma.hp + li; +#endif + return EBData{pebflag, m_eb_data.data()+EBData::real_data_size*li}; +} + std::unique_ptr makeEBFabFactory (const Geometry& a_geom, const BoxArray& a_ba, diff --git a/Src/EB/CMakeLists.txt b/Src/EB/CMakeLists.txt index 6f17526fcd..a6c82f743f 100644 --- a/Src/EB/CMakeLists.txt +++ b/Src/EB/CMakeLists.txt @@ -20,6 +20,7 @@ foreach(D IN LISTS AMReX_SPACEDIM) AMReX_EBMultiFabUtil.cpp AMReX_EBCellFlag.H AMReX_EBCellFlag.cpp + AMReX_EBData.H AMReX_EBDataCollection.H AMReX_EBDataCollection.cpp AMReX_MultiCutFab.H diff --git a/Src/EB/Make.package b/Src/EB/Make.package index 13f9d2c5ae..d6fde89055 100644 --- a/Src/EB/Make.package +++ b/Src/EB/Make.package @@ -12,6 +12,8 @@ CEXE_sources += AMReX_EBMultiFabUtil.cpp CEXE_headers += AMReX_EBCellFlag.H CEXE_sources += AMReX_EBCellFlag.cpp +CEXE_headers += AMReX_EBData.H + CEXE_headers += AMReX_EBDataCollection.H CEXE_sources += AMReX_EBDataCollection.cpp diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H index 557b14f7a4..4f107593fe 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H @@ -378,12 +378,8 @@ void mlebabeclap_gsrb (Box const& box, Array4 const& m1, Array4 const& m3, Array4 const& f0, Array4 const& f2, Array4 const& f1, Array4 const& f3, - 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, Array4 const& beb, + EBData const& ebdata, bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid, Box const& vbox, int redblack, int ncomp) noexcept { @@ -394,7 +390,8 @@ void mlebabeclap_gsrb (Box const& box, { if ((i+j+k+redblack) % 2 == 0) { - if (flag(i,j,k).isCovered()) + auto const flag = ebdata.get(i,j,k); + if (flag.isCovered()) { phi(i,j,k,n) = Real(0.0); } @@ -409,7 +406,7 @@ void mlebabeclap_gsrb (Box const& box, Real cf3 = (j == vhi.y && m3(i,vhi.y+1,k) > 0) ? f3(i,vhi.y,k,n) : Real(0.0); - if (flag(i,j,k).isRegular()) + if (flag.isRegular()) { Real gamma = alpha*a(i,j,k) + dhx * (bX(i+1,j,k,n) + bX(i,j,k,n)) @@ -428,19 +425,19 @@ void mlebabeclap_gsrb (Box const& box, } 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 kappa = ebdata.get(i,j,k); + Real apxm = ebdata.get(i ,j ,k); + Real apxp = ebdata.get(i+1,j ,k); + Real apym = ebdata.get(i ,j ,k); + Real apyp = ebdata.get(i ,j+1,k); Real fxm = -bX(i,j,k,n)*phi(i-1,j,k,n); Real oxm = -bX(i,j,k,n)*cf0; Real sxm = bX(i,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); + Real fcx = ebdata.get(i,j,k); + int jj = j + static_cast(std::copysign(Real(1.0),fcx)); + Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fxm = (Real(1.0)-fracy)*fxm + @@ -460,9 +457,9 @@ void mlebabeclap_gsrb (Box const& box, Real oxp = bX(i+1,j,k,n)*cf2; Real sxp = -bX(i+1,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); + Real fcx = ebdata.get(i+1,j,k); + int jj = j + static_cast(std::copysign(Real(1.0),fcx)); + Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fxp = (Real(1.0)-fracy)*fxp + @@ -482,9 +479,9 @@ void mlebabeclap_gsrb (Box const& box, Real oym = -bY(i,j,k,n)*cf1; Real sym = bY(i,j,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); + Real fcy = ebdata.get(i,j,k); + int ii = i + static_cast(std::copysign(Real(1.0),fcy)); + Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fym = (Real(1.0)-fracx)*fym + @@ -505,9 +502,9 @@ void mlebabeclap_gsrb (Box const& box, Real oyp = bY(i,j+1,k,n)*cf3; Real syp = -bY(i,j+1,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); + Real fcy = ebdata.get(i,j+1,k); + int ii = i + static_cast(std::copysign(Real(1.0),fcy)); + Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fyp = (Real(1.0)-fracx)*fyp + @@ -543,9 +540,9 @@ void mlebabeclap_gsrb (Box const& box, 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(vfrc(i,j,k)); + Real bctx = ebdata.get(i,j,k,0); + Real bcty = ebdata.get(i,j,k,1); + Real dx_eb = get_dx_eb(kappa); Real dg, gx, gy, sx, sy; if (std::abs(anrmx) > std::abs(anrmy)) { @@ -569,10 +566,12 @@ void mlebabeclap_gsrb (Box const& box, // In gsrb we are always in residual-correction form so phib = 0 Real dphidn = ( -phig)/dg; - Real feb = dphidn * ba(i,j,k) * beb(i,j,k,n); + Real ba = ebdata.get(i,j,k); + + Real feb = dphidn * ba * beb(i,j,k,n); rho += -vfrcinv*(-dh)*feb; - Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n); + Real feb_gamma = -phig_gamma/dg * beb(i,j,k,n); gamma += vfrcinv*(-dh)*feb_gamma; } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H index 174eef0af3..4c914b5965 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_3D_K.H @@ -557,14 +557,8 @@ void mlebabeclap_gsrb (Box const& box, Array4 const& f4, Array4 const& f1, Array4 const& f3, Array4 const& f5, - 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, + Array4 const& ccm, Array4 const& beb, + EBData const& ebdata, bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid, Box const& vbox, int redblack, int ncomp) noexcept { @@ -584,7 +578,8 @@ void mlebabeclap_gsrb (Box const& box, { if ((i+j+k+redblack) % 2 == 0) { - if (flag(i,j,k).isCovered()) + auto const flag = ebdata.get(i,j,k); + if (flag.isCovered()) { phi(i,j,k,n) = Real(0.0); } @@ -603,7 +598,7 @@ void mlebabeclap_gsrb (Box const& box, Real cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0) ? f5(i,j,vhi.z,n) : Real(0.0); - if (flag(i,j,k).isRegular()) + if (flag.isRegular()) { Real gamma = alpha*a(i,j,k) + dhx*(bX(i+1,j,k,n) + bX(i,j,k,n)) @@ -626,24 +621,26 @@ void mlebabeclap_gsrb (Box const& box, } 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 kappa = ebdata.get(i,j,k); + Real apxm = ebdata.get(i ,j ,k ); + Real apxp = ebdata.get(i+1,j ,k ); + Real apym = ebdata.get(i ,j ,k ); + Real apyp = ebdata.get(i ,j+1,k ); + Real apzm = ebdata.get(i ,j ,k ); + Real apzp = ebdata.get(i ,j ,k+1); Real fxm = -bX(i,j,k,n)*phi(i-1,j,k,n); Real oxm = -bX(i,j,k,n)*cf0; Real sxm = bX(i,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))); + auto fcx0 = ebdata.get(i,j,k,0); + auto fcx1 = ebdata.get(i,j,k,1); + int jj = j + static_cast(std::copysign(Real(1.0), fcx0)); + int kk = k + static_cast(std::copysign(Real(1.0), fcx1)); Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) - ? std::abs(fcx(i,j,k,0)) : Real(0.0); + ? std::abs(fcx0) : Real(0.0); Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) - ? std::abs(fcx(i,j,k,1)) : Real(0.0); + ? std::abs(fcx1) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm @@ -668,12 +665,14 @@ void mlebabeclap_gsrb (Box const& box, Real oxp = bX(i+1,j,k,n)*cf3; Real sxp = -bX(i+1,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))); + auto fcx0 = ebdata.get(i+1,j,k,0); + auto fcx1 = ebdata.get(i+1,j,k,1); + int jj = j + static_cast(std::copysign(Real(1.0),fcx0)); + int kk = k + static_cast(std::copysign(Real(1.0),fcx1)); Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) - ? std::abs(fcx(i+1,j,k,0)) : Real(0.0); + ? std::abs(fcx0) : 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); + ? std::abs(fcx1) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp @@ -699,12 +698,14 @@ void mlebabeclap_gsrb (Box const& box, Real oym = -bY(i,j,k,n)*cf1; Real sym = bY(i,j,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))); + auto fcy0 = ebdata.get(i,j,k,0); + auto fcy1 = ebdata.get(i,j,k,1); + int ii = i + static_cast(std::copysign(Real(1.0),fcy0)); + int kk = k + static_cast(std::copysign(Real(1.0),fcy1)); Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) - ? std::abs(fcy(i,j,k,0)) : Real(0.0); + ? std::abs(fcy0) : Real(0.0); Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) - ? std::abs(fcy(i,j,k,1)) : Real(0.0); + ? std::abs(fcy1) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym @@ -729,12 +730,14 @@ void mlebabeclap_gsrb (Box const& box, Real oyp = bY(i,j+1,k,n)*cf4; Real syp = -bY(i,j+1,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))); + auto fcy0 = ebdata.get(i,j+1,k,0); + auto fcy1 = ebdata.get(i,j+1,k,1); + int ii = i + static_cast(std::copysign(Real(1.0),fcy0)); + int kk = k + static_cast(std::copysign(Real(1.0),fcy1)); Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) - ? std::abs(fcy(i,j+1,k,0)) : Real(0.0); + ? std::abs(fcy0) : 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); + ? std::abs(fcy1) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp @@ -759,12 +762,14 @@ void mlebabeclap_gsrb (Box const& box, Real ozm = -bZ(i,j,k,n)*cf2; Real szm = bZ(i,j,k,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))); + auto fcz0 = ebdata.get(i,j,k,0); + auto fcz1 = ebdata.get(i,j,k,1); + int ii = i + static_cast(std::copysign(Real(1.0),fcz0)); + int jj = j + static_cast(std::copysign(Real(1.0),fcz1)); Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) - ? std::abs(fcz(i,j,k,0)) : Real(0.0); + ? std::abs(fcz0) : Real(0.0); Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) - ? std::abs(fcz(i,j,k,1)) : Real(0.0); + ? std::abs(fcz1) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm @@ -789,12 +794,14 @@ void mlebabeclap_gsrb (Box const& box, Real ozp = bZ(i,j,k+1,n)*cf5; Real szp = -bZ(i,j,k+1,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))); + auto fcz0 = ebdata.get(i,j,k+1,0); + auto fcz1 = ebdata.get(i,j,k+1,1); + int ii = i + static_cast(std::copysign(Real(1.0),fcz0)); + int jj = j + static_cast(std::copysign(Real(1.0),fcz1)); Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1)) - ? std::abs(fcz(i,j,k+1,0)) : Real(0.0); + ? std::abs(fcz0) : 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); + ? std::abs(fcz1) : Real(0.0); if (!beta_on_centroid && !phi_on_centroid) { fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp @@ -840,10 +847,10 @@ void mlebabeclap_gsrb (Box const& box, Real anrmx = dapx * anorminv; Real anrmy = dapy * anorminv; Real anrmz = dapz * anorminv; - 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(vfrc(i,j,k)); + Real bctx = ebdata.get(i,j,k,0); + Real bcty = ebdata.get(i,j,k,1); + Real bctz = ebdata.get(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)); @@ -874,10 +881,12 @@ void mlebabeclap_gsrb (Box const& box, + (gxy + gxyz) * phi(ii,jj,k,n) + (-gxyz) * phi(ii,jj,kk,n); + Real ba = ebdata.get(i,j,k); + Real dphidn = ( -phig)/dg; - Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n); + Real feb_gamma = -phig_gamma/dg * ba * beb(i,j,k,n); gamma += vfrcinv*(-dhx)*feb_gamma; - Real feb = dphidn * ba(i,j,k) * beb(i,j,k,n); + Real feb = dphidn * ba * beb(i,j,k,n); rho += -vfrcinv*(-dhx)*feb; } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp index 89258b75fe..f77cd8de90 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp @@ -203,13 +203,6 @@ MLEBABecLap::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs, 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; bool is_eb_dirichlet = isEBDirichlet(); @@ -279,19 +272,11 @@ MLEBABecLap::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs, 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& bebfab = (is_eb_dirichlet) ? m_eb_b_coeffs[amrlev][mglev]->const_array(mfi) : foo; + auto const& ebdata = factory->getEBData(mfi); + bool beta_on_centroid = (m_beta_loc == Location::FaceCentroid); bool phi_on_centroid = (m_phi_loc == Location::CellCentroid); @@ -307,10 +292,7 @@ MLEBABecLap::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs, AMREX_D_DECL(m1,m3,m5), AMREX_D_DECL(f0fab,f2fab,f4fab), AMREX_D_DECL(f1fab,f3fab,f5fab), - ccmfab, flagfab, vfracfab, - AMREX_D_DECL(apxfab,apyfab,apzfab), - AMREX_D_DECL(fcxfab,fcyfab,fczfab), - bafab, bcfab, bebfab, + ccmfab, bebfab, ebdata, is_eb_dirichlet, beta_on_centroid, phi_on_centroid, vbx, redblack, nc); });