diff --git a/Src/EB/AMReX_EB2_3D_C.H b/Src/EB/AMReX_EB2_3D_C.H index 8baf25d07b..1ffb7743ec 100644 --- a/Src/EB/AMReX_EB2_3D_C.H +++ b/Src/EB/AMReX_EB2_3D_C.H @@ -336,9 +336,19 @@ int coarsen_from_fine (int i, int j, int k, Box const& bx, int ngrow, const Box& xgbx = amrex::surroundingNodes(gbx,0); const Box& ygbx = amrex::surroundingNodes(gbx,1); const Box& zgbx = amrex::surroundingNodes(gbx,2); - const Box& exbx = amrex::convert(bx,IntVect(0,1,1)); // edge boxes - const Box& eybx = amrex::convert(bx,IntVect(1,0,1)); - const Box& ezbx = amrex::convert(bx,IntVect(1,1,0)); + Box exbx = amrex::convert(bx,IntVect(0,1,1)); // edge boxes + Box eybx = amrex::convert(bx,IntVect(1,0,1)); + Box ezbx = amrex::convert(bx,IntVect(1,1,0)); + // For face-centered data bx is nodal in face_dir, and the edge running along face_dir is + // indexed by the staggered cell it spans. There is one more of those than the conversion + // above leaves after turning face_dir from nodal into cell, so grow it back. + if (face_dir == 0) { + exbx.growHi(0,1); + } else if (face_dir == 1) { + eybx.growHi(1,1); + } else if (face_dir == 2) { + ezbx.growHi(2,1); + } const Box& exgbx = amrex::convert(gbx,IntVect(0,1,1)); const Box& eygbx = amrex::convert(gbx,IntVect(1,0,1)); const Box& ezgbx = amrex::convert(gbx,IntVect(1,1,0)); diff --git a/Src/EB/AMReX_EB2_Level.H b/Src/EB/AMReX_EB2_Level.H index a2b3bcc5f8..7612adfaf0 100644 --- a/Src/EB/AMReX_EB2_Level.H +++ b/Src/EB/AMReX_EB2_Level.H @@ -173,6 +173,20 @@ public: bool m_built = false; }; + /** + * \brief Index type of the face-centered edge centroids for edge direction \p idim. + * + * The edge is cell-typed in the direction it runs along and nodal in the others, as in the + * cell-centered path, except that \p face_dir is always nodal: that is how a box spans the + * staggered cells it owns. Used for both the level's data and the factory's, which have to + * agree - fillEdgeCentFC copies one into the other. + */ + [[nodiscard]] static IntVect fcEdgeType (int idim, int face_dir) noexcept { + IntVect edge_type{1}; + if (idim != face_dir) { edge_type[idim] = 0; } + return edge_type; + } + //! Check if face-centered data is available for face_dir (0=x, 1=y, 2=z). [[nodiscard]] bool hasFCData (int face_dir) const noexcept { return face_dir >= 0 && face_dir < AMREX_SPACEDIM && @@ -960,20 +974,31 @@ GShopLevel::buildFCData (G const& gshop, int face_dir, int max_grid_size) auto& fc_facecent = m_fc_data[face_dir]->m_facecent_fc; auto& fc_edgecent = m_fc_data[face_dir]->m_edgecent_fc; - // Index types of the face-centered data: + // Index types of the face-centered data. One rule covers all of it: face_dir is nodal, + // because that is how a box spans staggered cells - there is one more of them than there + // are cell-centered cells, and nodal is how AMReX says so. The remaining directions carry + // whatever type the cell-centered path would give. + // + // cellflag / volfrac / centroid / bndry* nodal in face_dir alone, i.e. (N,C,C) for + // face_dir=0, (C,N,C) for 1, (C,C,N) for 2. + // These are the staggered cells themselves. + // areafrac[idim] / facecent[idim] the same, for every idim. + // edgecent[idim] cell in idim and nodal elsewhere, as in the + // cell-centered path, with face_dir forced nodal. // - // cellflag / volfrac / centroid / bndry* nodal in face_dir, i.e. (N,C,C) for face_dir=0, - // (C,N,C) for 1, (C,C,N) for 2. These are the - // staggered cells themselves. - // areafrac[idim] / facecent[idim] the cell-centered base grids, the same for - // every idim. - // edgecent[idim] set by which two directions the edge spans. + // areafrac and facecent are kept uniform across idim on purpose. The cell-centered path + // gives each direction its own nodal index type (areafrac[0] is (N,C,C), areafrac[1] is + // (C,N,C), ...); carrying that over to a grid already staggered in face_dir would put three + // different index types in play at once for no real gain. Here all three are simply indexed + // by the staggered cell they belong to, and areafrac[idim](i,j,k) is that cell's low face in + // direction idim, with areafrac[idim](i+1,j,k) its high one. // - // areafrac and facecent are kept uniform on purpose. The cell-centered path gives each - // direction its own nodal index type (areafrac[0] is (N,C,C), areafrac[1] is (C,N,C), ...); - // carrying that over to a grid already staggered in face_dir would put three different index - // types in play at once for no real gain. Here all three are simply indexed by the staggered - // cell they belong to. + // Indexing by the staggered cell is why face_dir has to be nodal even here. The valid box + // then holds every staggered cell the box owns, which is what coarsen_from_fine below fills + // and what the fill*FC ParallelCopy can hand out. A cell-typed face_dir would stop one short + // and leave the last staggered cell of every box in the ghost region, where the copy cannot + // reach it: the covered-region mask in fill*FC would clear that plane with nothing able to + // write it back. // Define FC MultiFabs BoxArray fc_grids = amrex::convert(fc_base_grids, IntVect::TheDimensionVector(face_dir)); @@ -986,14 +1011,10 @@ GShopLevel::buildFCData (G const& gshop, int face_dir, int max_grid_size) fc_bndrynorm.define(fc_grids, fc_dmap, AMREX_SPACEDIM, ng); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - fc_areafrac[idim].define(fc_base_grids, fc_dmap, 1, ng); - fc_facecent[idim].define(fc_base_grids, fc_dmap, AMREX_SPACEDIM-1, ng); - IntVect edge_type{1}; - if (idim != face_dir) { - edge_type[face_dir] = 0; - edge_type[idim] = 0; - } - fc_edgecent[idim].define(amrex::convert(fc_base_grids, edge_type), fc_dmap, 1, ng); + fc_areafrac[idim].define(fc_grids, fc_dmap, 1, ng); + fc_facecent[idim].define(fc_grids, fc_dmap, AMREX_SPACEDIM-1, ng); + fc_edgecent[idim].define(amrex::convert(fc_base_grids, Level::fcEdgeType(idim,face_dir)), + fc_dmap, 1, ng); } // ********************************************************************** diff --git a/Src/EB/AMReX_EB2_Level.cpp b/Src/EB/AMReX_EB2_Level.cpp index fabbe5f3ee..a191d2a8b1 100644 --- a/Src/EB/AMReX_EB2_Level.cpp +++ b/Src/EB/AMReX_EB2_Level.cpp @@ -1111,7 +1111,8 @@ Level::fillAreaFracFC (Array const& areafrac, int face { const std::vector& pshifts = geom.periodicity().shiftIntVect(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - // FC areafrac is cell-typed on the base grids, indexed by the staggered cell + // the area fractions are indexed by the staggered cell, so the conversion reaches + // one further in face_dir, over the cells that straddle the covered interface BoxArray cov_fc = amrex::convert(cov, areafrac[idim]->ixType()); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Src/EB/AMReX_EBDataCollection.cpp b/Src/EB/AMReX_EBDataCollection.cpp index df257b94f4..be3ee23a80 100644 --- a/Src/EB/AMReX_EBDataCollection.cpp +++ b/Src/EB/AMReX_EBDataCollection.cpp @@ -144,16 +144,18 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level, m_bndrynorm = new MultiCutFab(fc_grids, a_dm, AMREX_SPACEDIM, ng, *m_cellflags); a_level.fillBndryNormFC(*m_bndrynorm, face_dir, m_geom); + // Index types follow the level's face-centered convention: face_dir is nodal, so that a + // box spans the staggered cells it owns, and the other directions carry the type the + // cell-centered path would give. See the comment in GShopLevel::buildFCData. + // + // The extra ghost layer on the area fractions is for callers that sweep a grown box and + // ask about the far face of its last cell, which is one index further still. for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - m_areafrac[idim] = new MultiCutFab(a_ba, a_dm, 1, m_ngrow[1]+1, *m_cellflags); - m_facecent[idim] = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM-1, ng, *m_cellflags); - IntVect edge_type{1}; - if (idim != face_dir) { - edge_type[face_dir] = 0; - edge_type[idim] = 0; - } - m_edgecent[idim] = new MultiCutFab(amrex::convert(a_ba, edge_type), a_dm, - 1, ng, *m_cellflags); + m_areafrac[idim] = new MultiCutFab(fc_grids, a_dm, 1, m_ngrow[1]+1, *m_cellflags); + m_facecent[idim] = new MultiCutFab(fc_grids, a_dm, AMREX_SPACEDIM-1, ng, *m_cellflags); + m_edgecent[idim] = new MultiCutFab( + amrex::convert(a_ba, EB2::Level::fcEdgeType(idim,face_dir)), a_dm, + 1, ng, *m_cellflags); } a_level.fillAreaFracFC(m_areafrac, face_dir, m_geom); diff --git a/Tests/EB/FCFactory/main.cpp b/Tests/EB/FCFactory/main.cpp index f960c3359d..2aad904909 100644 --- a/Tests/EB/FCFactory/main.cpp +++ b/Tests/EB/FCFactory/main.cpp @@ -15,6 +15,7 @@ #include #include +#include using namespace amrex; @@ -149,6 +150,34 @@ void main_main() } } + // Index types. The convention is one rule - face_dir nodal, the other directions as in the + // cell-centered path - and most of it has no observable effect on the values, because the + // plane it puts back inside the valid box holds what the ghost region held before. A value + // check cannot see that, so state the contract directly. The rule is written out here rather + // than taken from Level::fcEdgeType, so that the test is an independent statement of it. + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + const IndexType fc_type(IntVect::TheDimensionVector(dir)); + auto check = [&] (IndexType t, IndexType expect, const std::string& what) { + if (t != expect) { + ++nerrors; + amrex::Print() << "ERROR: dir=" << dir << " " << what << " index type is " + << t << ", expected " << expect << "\n"; + } + }; + check(fc_factories[dir]->getVolFrac().ixType(), fc_type, "volfrac"); + check(fc_factories[dir]->getMultiEBCellFlagFab().ixType(), fc_type, "cellflag"); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + check(fc_factories[dir]->getAreaFrac()[idim]->boxArray().ixType(), fc_type, + "areafrac[" + std::to_string(idim) + "]"); + check(fc_factories[dir]->getFaceCent()[idim]->boxArray().ixType(), fc_type, + "facecent[" + std::to_string(idim) + "]"); + IntVect edge_type{1}; + if (idim != dir) { edge_type[idim] = 0; } + check(fc_factories[dir]->getEdgeCent()[idim]->boxArray().ixType(), + IndexType(edge_type), "edgecent[" + std::to_string(idim) + "]"); + } + } + // Covered regions. With eb2.max_grid_size small enough that whole boxes fall inside the // body, the EB level has covered boxes, and the face-centered data must report them as // covered rather than as regular fluid. A face-centered cell is covered exactly when both @@ -160,30 +189,35 @@ void main_main() const MultiFab& volfrac_fc = fc_factories[dir]->getVolFrac(); const auto& ebflag = fc_factories[dir]->getMultiEBCellFlagFab(); auto const& afrac = fc_factories[dir]->getAreaFrac(); + auto const& ecent = fc_factories[dir]->getEdgeCent(); const IntVect fdir = IntVect::TheDimensionVector(dir); - ReduceOps reduce_op; - ReduceData reduce_data(reduce_op); - using ReduceTuple = typename ReduceData::Type; + ReduceOps reduce_op; + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename ReduceData::Type; for (MFIter mfi(volfrac_fc); mfi.isValid(); ++mfi) { auto const& vf = volfrac_fc.const_array(mfi); auto const& flag = ebflag[mfi].const_array(); auto const& cflag = ccflag[mfi].const_array(); - // the area fractions carry the cell-centered type, so enclosedCells of the - // face-centered valid box is exactly their valid range, and those indices are - // valid in the face-centered arrays too - const Box abx = amrex::enclosedCells(mfi.validbox()); + // every face-centered array is indexed by the staggered cell, so one box covers + // all of them + const Box abx = mfi.validbox(); const bool have_ap = afrac[0]->ok(mfi); AMREX_D_TERM(auto const& apx = have_ap ? afrac[0]->const_array(mfi) : Array4{};, auto const& apy = have_ap ? afrac[1]->const_array(mfi) : Array4{};, auto const& apz = have_ap ? afrac[2]->const_array(mfi) : Array4{};); + const bool have_ec = ecent[0]->ok(mfi); + AMREX_D_TERM(auto const& ecx = have_ec ? ecent[0]->const_array(mfi) : Array4{};, + auto const& ecy = have_ec ? ecent[1]->const_array(mfi) : Array4{};, + auto const& ecz = have_ec ? ecent[2]->const_array(mfi) : Array4{};); + const int ndj = (AMREX_SPACEDIM == 3) ? 2 : 1; reduce_op.eval(abx, reduce_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { const IntVect iv{AMREX_D_DECL(i,j,k)}; if (!(cflag(iv-fdir).isCovered() && cflag(iv).isCovered())) { - return {0,0,0,0}; + return {0,0,0,0,0}; } int nbad_vf = (vf(i,j,k) != Real(0.0)) ? 1 : 0; int nbad_flag = (!flag(i,j,k).isCovered()) ? 1 : 0; @@ -193,7 +227,22 @@ void main_main() || apy(i,j,k) != Real(0.0), || apz(i,j,k) != Real(0.0))) { nbad_ap = 1; } } - return {1, nbad_vf, nbad_flag, nbad_ap}; + // Every edge of a covered cell lies in the body and has to read -1. The + // edges running along a direction sit at the corners of the other two, and + // in the staggered direction the next index is this cell's far side: the + // low edge of the next staggered cell is the high edge of this one. Those + // +1 offsets reach one layer into the ghost region, which the factory's + // ghost cells cover and the fill's ParallelCopy writes. + int nbad_ec = 0; + if (have_ec) { + for (int di = 0; di < 2; ++di) { + for (int dj = 0; dj < ndj; ++dj) { + AMREX_D_TERM(if (ecx(i,j+di,k+dj) != Real(-1.0)) { nbad_ec = 1; }, + if (ecy(i+di,j,k+dj) != Real(-1.0)) { nbad_ec = 1; }, + if (ecz(i+di,j+dj,k) != Real(-1.0)) { nbad_ec = 1; }) + }} + } + return {1, nbad_vf, nbad_flag, nbad_ap, nbad_ec}; }); } @@ -202,10 +251,12 @@ void main_main() int nbad_vf = amrex::get<1>(rv); int nbad_flag = amrex::get<2>(rv); int nbad_ap = amrex::get<3>(rv); + int nbad_ec = amrex::get<4>(rv); ParallelDescriptor::ReduceIntSum(ncovered); ParallelDescriptor::ReduceIntSum(nbad_vf); ParallelDescriptor::ReduceIntSum(nbad_flag); ParallelDescriptor::ReduceIntSum(nbad_ap); + ParallelDescriptor::ReduceIntSum(nbad_ec); if (ncovered == 0) { ++nerrors; @@ -227,6 +278,11 @@ void main_main() amrex::Print() << "ERROR: dir=" << dir << " " << nbad_ap << " of " << ncovered << " covered cells have a nonzero area fraction\n"; } + if (nbad_ec > 0) { + ++nerrors; + amrex::Print() << "ERROR: dir=" << dir << " " << nbad_ec << " of " << ncovered + << " covered cells have an edge that is not covered\n"; + } } // Report