Skip to content

Face-centered EB factories: follow-up items #5635

Description

@WeiqunZhang

Review of the face-centered (FC) EB feature at 548f9407f0 (#5548). Everything below is still open.

Reproductions use the in-tree test unless stated otherwise:

cmake -S . -B build -G Ninja -DAMReX_SPACEDIM=3 -DAMReX_EB=ON -DAMReX_ENABLE_TESTS=ON \
      -DAMReX_MPI=OFF -DAMReX_FORTRAN=OFF -DCMAKE_BUILD_TYPE=Release
ninja -C build -j8 Test_EB_FCFactory_3d
T=build/Tests/EB/FCFactory/3d/Test_EB_FCFactory_3d
I=Tests/EB/FCFactory/inputs

Items marked [verified] were reproduced by running. Items marked [read] are from source
inspection only. Nothing here was exercised under MPI, on GPU, or in 2D.


Blockers — silent wrong data

1. Covered boxes are reported as fluid [verified]

Src/EB/AMReX_EB2_Level.H:861, Src/EB/AMReX_EB2_Level.cpp:1062-1068, 1175-1194

const BoxArray fc_base_grids = amrex::coarsen(refined_cc_level->m_grids, 2);

m_grids holds only the cut boxes. The refined level's m_covered_grids is dropped and never
re-applied, and none of the fill*FC functions mask covered regions the way Level::fillVolFrac
(Level.cpp:499-518) and Level::fillEBCellFlag (Level.cpp:461-478) do. A fully covered box
therefore keeps the setVal(1.0) default: volfrac 1, flag regular, in the middle of solid.

$ $T $I                          # default eb2.max_grid_size
FCFactory test PASSED
$ $T $I eb2.max_grid_size=16     # domain now has covered boxes
ERROR: dir=0 vol_error=12.80154696 exceeds 0.001
ERROR: dir=1 vol_error=12.80154696 exceeds 0.001
ERROR: dir=2 vol_error=12.80154696 exceeds 0.001
FCFactory test FAILED with 3 errors

The stock test only passes because its geometry produces no covered boxes at the default grid size.
Any production-size domain has them.

Direction: carry amrex::coarsen(refined_cc_level->m_covered_grids, 2) in FCData and apply it in
fillVolFracFC / fillEBCellFlagFC / fillAreaFracFC as the CC versions do; port the
fine_covered_grids zeroing loop from Level.cpp:209-251.

2. FC cell flags claim connectivity through solid [verified]

Src/EB/AMReX_EB2_Level.H:1138-1140

buildFCData ends by calling buildCellFlag(), but Level::buildCellFlag() (Level.cpp:402-426)
operates on m_cellflag / m_areafrac — this level's cell-centered data, not the FC data.
coarsen_from_fine only calls setRegular / setCovered / setSingleValued, which touch the low
bits and leave all 27 neighbour bits at their default of "connected".

Counting cut cells that report isConnected() toward a covered neighbour (64³ sphere, r=0.25):

CC  cut=4760  cut-cells linked to a COVERED neighbour = 0
FC0 cut=4912  cut-cells linked to a COVERED neighbour = 33296

Anything consuming flag.isConnected() on FC data — slope reconstruction, redistribution, MLEB
stencils — is reading connectivity through the boundary.

Secondary: because the call lands on the CC level, it is a pure side effect repeated
AMREX_SPACEDIM × nlev times. It is idempotent, so the CC data is not corrupted, but it costs a
FillBoundary plus a full-domain kernel pass each time.

Direction: build the FC flags from the FC area fractions — a face-dir-aware build_cellflag_from_ap
over m_areafrac_fc — and drop the CC buildCellFlag() call.

3. No coarsenability check: overlapping coarse boxes and uninitialized data [verified]

Src/EB/AMReX_EB2_Level.H:861, versus the CC guard at Level.H:781-802

The CC coarse-level path checks fine_grids.coarsenable(2, /*min_width*/8) and otherwise routes
through prepareForCoarsening. buildFCData coarsens unconditionally. When the refined grids are
not coarsenable, amrex::coarsen yields an overlapping BoxArray:

n_cell=65, eb2.max_grid_size=33
  CC level BA:  n=8   isDisjoint=1   coarsenable(2,8)=0
  FC0 level BA: n=56  isDisjoint=0
  FC1 level BA: n=56  isDisjoint=0
  FC2 level BA: n=56  isDisjoint=0
  FC2 m_volfrac_fc.contains_nan(valid) = 1

NaN in the valid region of the level's FC data, with no abort and no warning. Which direction
surfaces it varies with the configuration; the non-disjoint BoxArray is present in all three.

Direction: assert refined_cc_level->m_grids.coarsenable(2,8), or make the transient level's grids
coarsenable by construction — derive its max_grid_size from the coarse level (even, doubled) rather
than reusing EB2::max_grid_size at 2× resolution.


Should-fix

4. Coarsening failure is computed, reduced, then discarded [read]

Src/EB/AMReX_EB2_Level.H:1132-1144

if (!error) { buildCellFlag(); }
// return error;                       <-- leftover
m_fc_data[face_dir]->m_built = true;   <-- unconditional

Compare Level.cpp:395-399 (return error;) and Level.H:793-801 (m_ok = (ierr == 0);, with a
rebuild-from-shop fallback at IndexSpaceI.H:53-64). A failed FC coarsening leaves m_built = true
and hasFCData() true, and the factory hands the user partially written data. m_ok is never
touched.

5. FC data is valid only 8 cells outside the domain, but the extension trusts the CC value [verified]

Src/EB/AMReX_EB2_Level.H:839-841, Src/EB/AMReX_EBDataCollection.cpp:170-172

The transient level is built with ngrow = m_ngrow[0] > 0 ? m_ngrow[0] : 4, rounded up to 16, so the
coarsened FC grids reach only 8 cells past the domain while the CC level reaches 16.
extendDataOutsideDomain is then handed a_level.nGrowVect() — the CC value — so it treats bands
9–16 as already valid and skips them.

Plane at z=0.5 with solid above, FC factory with ngrow=20, sampling volfrac(i,32,63):

i = -1 .. -8 :  CC 0   FC 0     correct (covered)
i = -9 .. -20:  CC 0   FC 1     wrong

Related: fillVolFracFC (Level.cpp:1062-1068) does setVal(1.0) then ParallelCopy with
src_nghost = 0, so any destination region the source's valid box does not cover keeps "regular
fluid". A caller requesting more ghost cells than the FC data has gets silently wrong values rather
than an error.

Direction: add Level::nGrowVectFC(face_dir) and pass that to extendDataOutsideDomain; assert or
clamp when the requested a_ngrow exceeds the FC valid extent.

6. BuildFC() costs ~7× Build() in time and ~4× in peak memory [verified]

Src/EB/AMReX_EB2_IndexSpaceI.H:211-219, Src/EB/AMReX_EB2_Level.H:838-846

The loop nest is face_dir outer, ilev inner, and each buildFCData builds its own transient
2×-refined GShopLevel from the shop — AMREX_SPACEDIM × nlev refined levels where one per level
would serve all directions. The face-dir-independent multi-cut scan (Level.H:871-928) is repeated
the same number of times.

n_cell=128, single level:  EB2::Build   0.50 s / VmHWM 1052 MB
                           EB2::BuildFC 3.50 s / VmHWM 4629 MB     (7.0x time, 4.4x peak RSS)
n_cell=128, mcl=3:         EB2::Build   0.62 s     BuildFC 4.87 s

Direction: swap the loop order and build the refined level once per level; hoist the multi-cut check
out of the per-direction path; consider letting callers request only the directions they need.

7. extend_domain_face and num_coarsen_opt are hard-coded [read]

Src/EB/AMReX_EB2_Level.H:841

max_grid_size, m_ngrow[0] > 0 ? m_ngrow[0] : 4, true, 0);
// xxxxx TODO: should we use IndexSpace's member variables extend_domain_face and num_crse_opt?

IndexSpaceImp already stores m_extend_domain_face and m_num_coarsen_opt (AMReX_EB2.H:121-122).
A user who built with eb2.extend_domain_face=0 gets FC data built from a different geometry than
their CC data, with no diagnostic. The 4 and the direction-0-only m_ngrow[0] probe are
undocumented.

8. BuildFC() is a silent no-op for STL and checkpoint index spaces [read]

Src/EB/AMReX_EB2.H:70, Src/EB/AMReX_EB2.cpp:325-330, Src/EB/AMReX_EBDataCollection.cpp:100

IndexSpace::buildAllFCData() defaults to {} and only IndexSpaceImp overrides it — neither
IndexSpaceSTL nor IndexSpaceChkptFile does. EB2::BuildFC() returns successfully and the user
finds out later through
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_level.hasFCData(face_dir), "EBDataCollection: FC data not available for face_dir").
STL is a mainstream EB source.

Direction: make the base method abort with a clear message, or implement it for STL (which has a
shop) and abort for the checkpoint case.

9. m_shift is ignored by every fill*FC [read]

Src/EB/AMReX_EB2_Level.cpp:1067, 1079, 1108, 1126, 1144, 1162, 1187, 1205, 1237

Every CC counterpart passes -m_shift to ParallelCopy (e.g. Level.cpp:488) and shifts the
periodic shift vectors. The FC versions pass neither, so after IndexSpace::setShift
(IndexSpaceI.H:198-207) the CC and FC data describe geometries offset from one another.

10. The FC factory reports an all-zero level set [read]

Src/EB/AMReX_EBDataCollection.cpp:115-118

A nodal MultiFab is allocated and setVal(0.0). Zero means "exactly on the surface" everywhere —
the worst available sentinel. It also feeds extendDataOutsideDomain
(EBDataCollection.cpp:245-262), where ls_a(...) >= 0 counts as covered, so a cut cell outside
level_domain is turned covered. EBFArrayBox::getLevelSetData() will hand it to users.

Direction: leave m_levelset null for FC and make getLevelSet() abort for an FC factory, or fill a
clearly invalid sentinel; guard the level-set branch of extendDataOutsideDomain for FC.

11. extendDataOutsideDomain assumes CC index types [read]

Src/EB/AMReX_EBDataCollection.cpp:171, 194-197, 290-300

lev_ap_domain[idim] = surroundingNodes(level_domain, idim) and the
apbx.smallEnd(idim) == nbx.smallEnd(idim) trimming assume m_areafrac[idim] is idim-nodal. Under
the FC convention they are cell-typed on a_ba (EBDataCollection.cpp:148-149), so the extension
applies the wrong boxes whenever it runs.

12. Test coverage [partly verified]

Tests/EB/FCFactory/main.cpp

  • The checks are NaN-blind. vf_min < 0.0 || vf_min > 1.0 (:95) and vol_error > vol_tol
    (:145) are both false for NaN. With eb2.max_grid_size=33 n_cell=65 sphere_radius=0.49 the test
    reports PASSED while the level's FC volfrac contains NaN over its valid region (item 3). Add an
    explicit contains_nan() check.
  • Only volfrac and "some cut cells exist" are validated. areafrac, facecent, edgecent,
    centroid, bndrycent, bndrynorm and flag connectivity are never checked — which is how items 1
    and 2 went unnoticed.
  • Not covered: covered boxes (any eb2.max_grid_size smaller than the domain), periodic domains,
    max_coarsening_level > 0, ngrow > 1, MPI (GNUmakefile:8 sets USE_MPI=FALSE), 2D,
    all-regular.
  • A strong oracle is available and cheap: build a second EB on the grid shifted by dx/2 in
    face_dir and compare index by index. Doing this confirms the current staggered offset and index
    conventions are correct (agreement to volfrac 1.5e-2, areafrac 2.2e-2, edgecent 9e-12), so it
    is worth locking in as a regression test.
  • main.cpp:129 uses std::numbers::pi_v but includes only <cmath>; add <numbers> (it compiles
    today by transitive include on libstdc++).

13. No documentation [verified: grep -rn "BuildFC\|faceDir" Docs/ is empty]

Nothing under Docs/sphinx_documentation/. Users cannot discover EB2::BuildFC() or the FC factory
constructor, and in particular nothing states the index-type convention: areafrac[idim] and
facecent[idim] live on the cell-centered base grid indexed by the staggered cell, rather than each
being nodal in its own direction as in the CC path. Also worth stating that FC data is a coarsening of
a 2×-refined EB, so it differs from a directly computed staggered EB at the O(10⁻²) level.


Nice-to-have

14. areafrac ghost count breaks the surrounding assertions [read]

Src/EB/AMReX_EBDataCollection.cpp:148facecent/edgecent use ng = m_ngrow[2], areafrac uses
m_ngrow[1]+1. The +1 is load-bearing (the topmost staggered cell's low face lands at domhi+1,
a ghost index once areafrac is cell-typed) but unexplained, and it can exceed m_ngrow[0], which
the AMREX_ALWAYS_ASSERTs at :123 and :134 exist to prevent — with the stock {1,1,1} it already
does (2 > 1).

15. ng = 3 on the FC MultiFabs, but only 2 layers are written and 0 read [read]

Src/EB/AMReX_EB2_Level.H:962coarsen_from_fine is driven over grow(bx,2) (Level.H:1040, 1107), so the third ghost layer is never written, and every fill*FC ParallelCopy uses
src_nghost = 0, so no ghost layer is ever read. The CC equivalent is 2 (Level.cpp:167).

16. Multi-cut geometries hard-abort inside BuildFC [read]

Src/EB/AMReX_EB2_Level.H:926-928amrex::Abort("GShopLevel::buildFCData: MVMC error"). The CC path
returns the error and lets the index space fall back to a direct shop build (IndexSpaceI.H:53-64),
and EB2::BuildMultiValuedMultiCut exists for these geometries. Here a CC build that already
succeeded is followed by an abort naming neither the level nor a remedy.

17. Ordering against addFineLevels / addRegularCoarseLevels is undocumented [read]

Src/EB/AMReX_EB2_IndexSpaceI.H:122-148, 150-195 — both prepend levels; existing levels keep their
m_fc_data, new ones have none, so hasFCData() is false and a later factory asserts. Nothing states
that BuildFC() must come last.

18. Dead code, stale comments, duplication [read]

  • Level.H:833, 846Level const* refined_cc_level = nullptr; and the STEP1: Get or build header
    (:830) are leftovers from a version that could reuse an existing level; it is now always the
    transient one.
  • Level.H:1142// return error;. Level.H:1145 — "If transient_level exists, it's automatically
    destroyed here"; it always exists.
  • Level.H:836-837 — the comment references an external project ("ERF EB-FC development").
  • Level.H:169-172Level::getFCData() is public, unused, and exposes internal MultiFabs.
  • EBDataCollection.H:99-100m_is_fc and m_face_dir are written and never read.
  • EBFabFactory.cpp:95-172 duplicates :14-93 almost verbatim, and the copy drops the comment
    (:26-28) explaining why const_arrays() is called, leaving what reads as dead code. Worth
    factoring into a shared helper.
  • AMReX_EB2.H:69, 279-280 and Level.H:243 say "all three face directions" / "[2]=z-face", which is
    wrong in 2D.

Checked, no problem found

Assertions plus AMReX_BOUND_CHECK=ON on the stock test and several off-nominal configurations
produced no Array4 bounds violations. The staggered offset (ii = 2i-1) and the FC index-type
conventions were validated against an independently constructed staggered EB and are correct.
build_cellflag_from_ap is idempotent, so the misplaced buildCellFlag() in item 2 wastes work but
does not corrupt the cell-centered data.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions