EB is built in the incflo constructor from the inputs' prob_lo/prob_hi, but ReadCheckpointFile resets the physical domain from the checkpoint header with no check: a restart with edited prob_lo/prob_hi keeps the EB at the wrong physical location

[WeiqunZhang]

Severity: low · Category: io-restart
Locations: src/incflo.cpp:20-24, src/embedded_boundaries/embedded_boundaries.cpp:9-22, src/utilities/io.cpp:124-125, src/utilities/io.cpp:194-199, src/incflo.cpp:251-256
Based on commit 46de3367 (line numbers refer to that tree).

The defect

Order of operations on a restart:

1. incflo::incflo() calls ReadParameters() and then MakeEBGeometry() (incflo.cpp:20-24). Every builder
   evaluates its implicit function on geom.back(), i.e. the Geometry that AmrCore made from the
   inputs geometry.prob_lo/prob_hi (EB2::Build(gshop, geom.back(), ...)).
2. InitData → ReadCheckpointFile reads prob_lo/prob_hi from the checkpoint header (io.cpp:124-125,
   178-189) and installs them: Geometry::ResetDefaultProbDomain(rb); SetGeometry(lev, Geometry(Geom(lev).Domain(), rb, ...))
   (io.cpp:194-199).
3. MakeNewLevelFromScratch builds the level factories with makeEBFabFactory(geom[lev], ...)
   (incflo.cpp:251-256). AMReX finds the EB level by the index-space Domain() box
   (IndexSpaceImp::getLevel, AMReX_EB2_IndexSpaceI.H:97-104), which is unchanged, so no error is raised.

If the inputs used for the restart carry a different geometry.prob_lo/prob_hi than the run that wrote the
checkpoint, the cut cells are those of the inputs' physical domain (the implicit function was evaluated
there) while the state, the BCs and every coordinate-based routine now use the checkpoint's domain. Nothing
compares the two. The only EB path with a check is chkptfile (AMReX asserts the geom_chk prob domain
against geom, AMReX_EB_chkpt_file.cpp:169-172), and that compares against the inputs geometry too, so it
does not help here.

Why it matters

A user who restarts with, say, a doubled prob_hi (a common way to rescale a case) gets a run whose EB body
has silently moved and been rescaled in index space relative to the checkpointed flow, with no diagnostic.
The non-EB fields are protected because the checkpoint wins; the EB is the one object that is not rebuilt
from the checkpoint. Defensive gap, not a bug in a shipped deck.

How to reach it

EB build. Run test_2d/benchmark.channel_cylinder with amr.check_int = 5 max_step = 5, then restart with
amr.restart = chk00005 geometry.prob_hi = 2.4 0.8. The sphere (center 0.15, 0.2, radius 0.05 in inputs
coordinates) is now cut at index cells 6/16 instead of 12/16 of the domain; the restarted plotfile shows the
obstacle at physical (0.075, 0.1) inside the checkpoint's [0,1.2]x[0,0.4] domain.

Suggested fix

Refuse the restart when the header's prob domain differs from the inputs' and an EB was built. Diff checked
with git apply --check (audit/notes/T7-scratch/027.diff). EBFactory(0) does not exist yet at this point,
so the check goes through the EB2 index space:

--- a/src/utilities/io.cpp
+++ b/src/utilities/io.cpp
@@ -2,6 +2,9 @@
 #include <AMReX_PlotFileUtil.H>
 #include <AMReX_buildInfo.H>
 #include <incflo.H>
+#ifdef AMREX_USE_EB
+#include <AMReX_EB2.H>
+#endif
 
 using namespace amrex;
 
@@ -192,6 +195,20 @@
 
     // Set up problem domain
     RealBox rb(prob_lo, prob_hi);
+#ifdef AMREX_USE_EB
+    // The EB index space was built in the incflo constructor from the inputs'
+    // geometry.prob_lo/prob_hi; the physical domain is about to be reset from the
+    // checkpoint header, so refuse a mismatch instead of silently moving the EB.
+    if (!EB2::IndexSpace::top().getLevel(Geom(0)).isAllRegular()) {
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            if (std::abs(prob_lo[idim] - Geom(0).ProbLo(idim)) > Real(1.e-12) ||
+                std::abs(prob_hi[idim] - Geom(0).ProbHi(idim)) > Real(1.e-12)) {
+                amrex::Abort("ReadCheckpointFile: geometry.prob_lo/prob_hi differ from the checkpoint, "
+                             "but the EB geometry was already built from the inputs");
+            }
+        }
+    }
+#endif
     Geometry::ResetDefaultProbDomain(rb);
     for (int lev = 0; lev <= max_level; ++lev) {
         SetGeometry(lev, Geometry(Geom(lev).Domain(), rb, Geom(lev).CoordInt(),

An alternative is to read the header before MakeEBGeometry and build the EB from the checkpoint's domain;
that is a larger restructuring of the constructor and is left to the maintainers.

[asalmgren]

fixed by #267