From 7284587ff6e96a0bd67e662a1bc006c62f929551 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 19 Nov 2024 18:07:33 -0800 Subject: [PATCH] EBData: New class containing Array4's for all EB data This can make the code cleaner because we can get access to all available EB data in one object. This also reduces many GPU kernels' register pressure because the GPU device lambdas are much smaller without the need to capture multiple Array4's. For example, the 3D EBABecLap's gsrb kernel's number of registers reduces from 198 to 130. Add a new function that returns a random point on EB for WarpX. --- Src/EB/AMReX_EBData.H | 56 ++++++++++ Src/EB/AMReX_EBDataCollection.H | 3 + Src/EB/AMReX_EBFabFactory.H | 5 + Src/EB/AMReX_EBFabFactory.cpp | 78 ++++++++++++- Src/EB/CMakeLists.txt | 1 + Src/EB/Make.package | 2 + .../MLMG/AMReX_MLEBABecLap_2D_K.H | 59 +++++----- .../MLMG/AMReX_MLEBABecLap_3D_K.H | 103 ++++++++++-------- .../MLMG/AMReX_MLEBABecLap_F.cpp | 24 +--- 9 files changed, 232 insertions(+), 99 deletions(-) create mode 100644 Src/EB/AMReX_EBData.H diff --git a/Src/EB/AMReX_EBData.H b/Src/EB/AMReX_EBData.H new file mode 100644 index 00000000000..ebfe03467d7 --- /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 6a8556dc9a5..bf2d8f15494 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 c9db5f39473..7d20e843547 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 f5154578d6f..4d4f9dcee76 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 6f17526fcd6..a6c82f743fd 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 13f9d2c5ae3..d6fde890556 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 557b14f7a4d..1c4e6153927 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 * ba * 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 174eef0af33..4c914b59655 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 89258b75fe9..f77cd8de900 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); });