ApplyPBLHSmoothing reads outside the PBLH FArrayBox and is decomposition-dependent

[asalmgren]

ApplyPBLHSmoothing (Source/PBL/ERF_PBLModels.H:494) reads outside the FArrayBox it was handed, and its results depend on the domain decomposition. This was found while reviewing #3972, which explicitly lists it under "Not changed here"; opening it so the problem is tracked independently of that PR.

Where

ApplyPBLHSmoothing is called from two places, both under if (turbChoice.enable_pblh_smoothing):

• Source/PBL/ERF_ComputeDiffusivityMRF.cpp:558
• Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp:982

Both call sites pass the per-tile pbl_height_corrector FAB together with xybx, where

const Box& gbx  = mfi.growntilebox(IntVect(1,1,0));
const Box xybx  = PerpendicularBox<ZDir>(gbx, IntVect{0, 0, 0});
FArrayBox pbl_height_corrector(xybx, 1, The_Async_Arena());

so the FAB is allocated on exactly xybx — the valid tile grown by one cell in x and y, with no further ghost cells.

The out-of-bounds read

The smoothing kernel iterates over all of xybx and clamps its 5-point stencil only at the domain edge:

const int i_xlo = (i > dom_lo.x) ? (i - 1) : i;

For any box or tile whose low corner is interior, xybx.smallEnd(0) == vbx.smallEnd(0) - 1, which is strictly greater than dom_lo.x. The clamp therefore does not fire, and the kernel reads pblh(xybx.smallEnd(0) - 1, j, 0) — one cell outside the FAB it was given. The same applies at the upper edge and in y. Only a box that touches the domain boundary is accidentally safe there, because the grown ring puts i at or below dom_lo.x.

Consequences:

• Under AMReX_BOUND_CHECK / DEBUG=TRUE, any run with erf.enable_pblh_smoothing = true and more than one box (or more than one OpenMP tile) aborts.
• In an optimized build it reads whatever is adjacent in the arena, so the smoothed PBLH — and through it K_turb and the countergradient terms — is nondeterministic.

Decomposition dependence, independent of the OOB read

Even with the stencil clamped to the FAB, passes > 1 cannot be correct as written. Pass 1 produces valid smoothed values only on the interior of xybx, i.e. on the valid tile; the one-cell ring is itself computed from data that is not there. Pass 2 needs smoothed neighbours at vbx ± 1, which is exactly that corrupted ring. So with erf.pblh_smoothing_passes at its default of 1 the valid region happens to be recoverable once the clamp is fixed, but any larger value gives answers that change with the box/tile layout.

Not covered by CI

enable_pblh_smoothing defaults to false (Source/DataStructs/ERF_TurbStruct.H:807) and no regression test sets it. The tiling-parity tests added in #3972 do not enable it either, so a regression here would not be caught.

Suggested fix

Either

1. give the smoothing a proper halo — grow the temporary's source region, or run the smoothing on a MultiFab over the whole level with FillBoundary between passes, so neighbour values are real; or
2. restrict the ParallelFor to the valid (ungrown) region and clamp the stencil to the FAB's own box in addition to the domain, and reject pblh_smoothing_passes > 1 until (1) is done.

Option 2 is the smaller change and makes single-pass smoothing correct and decomposition-independent. Whichever is chosen, a regression test with erf.enable_pblh_smoothing = true run at two different max_grid_size values would keep it honest.