EB: make face-centered EB data nodal in face_dir throughout

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));

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<G>::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<G>::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);
     }
 
     // **********************************************************************

Src/EB/AMReX_EB2_Level.cpp

@@ -1111,7 +1111,8 @@ Level::fillAreaFracFC (Array<MultiFab*,AMREX_SPACEDIM> const& areafrac, int face
     {
         const std::vector<IntVect>& 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())

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);

Tests/EB/FCFactory/main.cpp

@@ -15,6 +15,7 @@
 #include <AMReX_Reduce.H>
 
 #include <cmath>
+#include <string>
 
 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<ReduceOpSum,ReduceOpSum,ReduceOpSum,ReduceOpSum> reduce_op;
-        ReduceData<int,int,int,int> reduce_data(reduce_op);
-        using ReduceTuple = typename ReduceData<int,int,int,int>::Type;
+        ReduceOps<ReduceOpSum,ReduceOpSum,ReduceOpSum,ReduceOpSum,ReduceOpSum> reduce_op;
+        ReduceData<int,int,int,int,int> reduce_data(reduce_op);
+        using ReduceTuple = typename ReduceData<int,int,int,int,int>::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<Real const>{};,
                          auto const& apy = have_ap ? afrac[1]->const_array(mfi) : Array4<Real const>{};,
                          auto const& apz = have_ap ? afrac[2]->const_array(mfi) : Array4<Real const>{};);
+            const bool have_ec = ecent[0]->ok(mfi);
+            AMREX_D_TERM(auto const& ecx = have_ec ? ecent[0]->const_array(mfi) : Array4<Real const>{};,
+                         auto const& ecy = have_ec ? ecent[1]->const_array(mfi) : Array4<Real const>{};,
+                         auto const& ecz = have_ec ? ecent[2]->const_array(mfi) : Array4<Real const>{};);
+            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