Skip to content

ConvertCheckpointGrids::ConvertData: refine path takes interpolation slopes from ghost cells that still hold the setVal(10.) sentinel at coarse-fine and physical boundaries #247

Description

@WeiqunZhang

Location: Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp:686, Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp:720, Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp:730
Severity: Low — O(dx) error confined to the row of refined cells next to walls and coarse-fine boundaries, bounded by the slope limiter; utility code
Category: Correctness
Based on commit bb697bf5 (line numbers refer to that tree).

Problem

PR #205 (issue #191) added FillBoundary(cgeom.periodicity()) for the source copies, so
periodic ghosts are now valid. Two other kinds of ghost cell are still never filled:

// Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp:686-691, 720-721
      MultiFab * NewData_src = new MultiFab(save_grids_state,dm,ncomps,ngrow_loc);
      NewData_src -> setVal(10.);
      NewData_src -> copy(*(falRef_src.state[n].new_data),0,0,ncomps,0,ngrow_loc);
      ...
        NewData_src->FillBoundary(cgeom.periodicity());

copy with src_nghost=0 and FillBoundary only fill ghosts that overlap another box of
the same level. Ghost cells adjacent to coarser cells (every box edge of levels >= 1) and
ghosts outside a non-periodic domain boundary keep the value 10. cell_cons_interp
(CellConservativeLinear::CoarseBox grows by one) reads them for the slopes: with
u_{i+1}=10 next to u ~ 1 the MC-limited slope becomes 2*(u_i - u_{i-1}) or 0 instead
of the centred difference, so the refined children of every boundary cell are shifted by
O(dx). The tool already includes AMReX_Extrapolater.H (line 34) but never calls it;
Extrapolater::FirstOrderExtrap fills exactly these cells (crsebnd and physbnd share mask
value 0, AMReX_Extrapolater.H:20-21) from neighbouring valid data.

A second, latent hazard in the same loop: the source MultiFabs live on dm (line 596) and the
target ones on an independently built dm_trgt (line 638), yet both are indexed through one
MFIter over the target ((*NewData_src)[mfi], lines 720-740). This is safe only because the
tool's GNUmakefile pins USE_MPI=FALSE; a USE_MPI=TRUE build compiles and would misindex
whenever the two space-filling-curve maps differ. Building the targets on dm (or a
ParallelCopy onto dm_trgt before interpolating) removes the assumption.

Impact

Any checkpoint with walls/inflow/outflow (e.g. Exec/run2d lid-driven cavity, channel decks)
or with amr.max_level>0, converted with interp_kind=refine: all cell-centred state types
(velocity, density, tracer, gradp, divu, averages) carry an O(dx) error in the first row of
fine cells along those boundaries. Nodal pressure is unaffected (fine node at ic*ratio uses
only crse(ic); interior fine nodes use nodes the box owns).

Suggested fix

Extrapolate into the remaining ghost cells after the periodic fill, for cell-centred data.

--- a/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp
+++ b/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp
@@ -713,13 +713,22 @@
         // The copies above are non-periodic, so ghost cells at periodic
         // domain boundaries still hold the setVal(10.) sentinel. Fill them
         // from valid data before computing the interpolation slopes.
-        // NOTE: ghosts at non-periodic physical boundaries still hold the
-        //       sentinel -- the state descriptors (and hence the real BCs)
-        //       are never restored by this tool.
         //
         NewData_src->FillBoundary(cgeom.periodicity());
         OldData_src->FillBoundary(cgeom.periodicity());
 
+        //
+        // Ghost cells not covered by this level's grids (coarse-fine
+        // boundaries on levels > 0 and non-periodic domain boundaries)
+        // still hold the sentinel. Fill them by first-order extrapolation
+        // so the slope stencil sees real data. Nodal data (pressure) only
+        // reads the nodes it owns, so it needs no fill.
+        //
+        if (NewData_src->is_cell_centered()) {
+          Extrapolater::FirstOrderExtrap(*NewData_src, cgeom, 0, ncomps);
+          Extrapolater::FirstOrderExtrap(*OldData_src, cgeom, 0, ncomps);
+        }
+
         for (MFIter mfi(*NewData_trgt); mfi.isValid(); ++mfi)
         {
           FArrayBox& ffab = (*NewData_trgt)[mfi];

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

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