MRF, YSUNew and YSU: stay inside MFIter tiles (tiled runs abort or store a garbage PBLH); restore #3961's HOL_abs bases

[hgopalan]

Summary

MRF, YSUNew and YSU break when MFIter tiling splits a box. On CPU that is the default, because AMReX's CPU fabarray.mfiter_tile_size is 1024000 8 8.

• MRF and YSUNew compute the corrected PBL height in a per-tile planar FArrayBox (PerpendicularBox<ZDir>(mfi.growntilebox(IntVect(1,1,0))) under TileNoZ()), then copy it into pblh_mf for SurfLayer->set_pblh with a ParallelFor over mfi.validbox(). Once a box is 16 or more cells wide in y, the valid box extends past the tile, and the copy reads outside the tile's array.
  • Debug and AMReX_BOUND_CHECK builds abort at step 1, e.g. (0,8,0,0) is out of bound (-1:16,-1:7,0:0,0:0) in ComputeDiffusivityMRF.
  • Release builds do not abort. Each tile writes the whole valid box, so SurfaceLayer's pblh holds out-of-bounds values for every tile but the last. That field feeds the 2D plotfile pblh and the Beljaars w* term (erf.most.include_wstar). The eddy diffusivities are unaffected.
• YSU iterates with TilingIfNotGPU(), which on CPU also tiles in z. Its full-column AMREX_ALWAYS_ASSERT (sbx must span the domain in z) therefore aborts at step 1 whenever nz ≥ 16.

Most MRF/YSUNew canonical decks set fabarray.mfiter_tile_size = 1024 1024 1024, which is why this has not shown up.

Changes

• ERF_ComputeDiffusivityMRF.cpp, ERF_ComputeDiffusivityYSUNew.cpp: loop over mfi.tilebox() for the pblh_mf copy. Under TileNoZ the tile spans the full column, and pblh_mf has no ghost cells. The tiles therefore partition the valid box, and every cell is written exactly once before set_pblh copies the valid cells.
• ERF_ComputeDiffusivityYSU.cpp: TileNoZ() instead of TilingIfNotGPU(), like the other column schemes.
• ERF_ComputeDiffusivityMRF.cpp: restore the 1 + 16*HOL_abs bases from Keep the MRF scheme from raising invalid floating-point flags it never uses #3961 in the Pass 4 wstar recompute and in the K-profile phiM/phit. PBLH Corrections #3486 replaced them with pow(max(1 - 16*HOL_bounded, 0.01), ...), which left HOL_abs unused in Pass 4. development therefore builds with warning: unused variable 'HOL_abs', which the CI warning check counts. It also brings back the pow-of-a-negative-base form that Keep the MRF scheme from raising invalid floating-point flags it never uses #3961 removed for amrex.fpe_trap_invalid runs. On the unstable arm HOL ≤ 0, so max(1 - 16*HOL, 0.01) equals 1 - 16*HOL = 1 + 16*|HOL| and the results do not change.
• New regression tests: ABL_MRF_Tiling, ABL_YSUNew_Tiling, ABL_YSU_Tiling.
  • They use add_test_tiling_parity in Tests/CTestList.cmake, Tests/RunTilingParity.cmake, and the deck Tests/test_files/ABL_MRF_Tiling. They are registered only when ERF_ENABLE_MPI is on and not on Windows, since the script runs mpiexec and the plotfile tools through execute_process.
  • The deck is 32×32×32 (dx = dy = 200 m, dz = 32 m) on four 16×16×32 boxes, with a surface layer and deterministic sinusoidal u/v/θ perturbations, run for 10 steps. The sounding has a 6 K inversion at 150–250 m and a 5 m/s wind, and MOST uses local averaging, so the PBL height, u*, t* and L differ from column to column (see below).
  • Each test runs it with fabarray.mfiter_tile_size = 1024000 8 8 (two tiles per box in y) and again with tiling off.
  • amrex_fcompare then requires the 3D plotfiles (ρ, u, v, w, p, θ, Kmv, Khv, Lturb) and the 2D plotfiles (pblh, u_star, t_star, Olen) to agree. No gold files are needed.
  • VARYING_3D / VARYING_2D name fields that must take more than one value in the untiled run (amrex_fextrema, min < max), so the test cannot pass by comparing a constant with itself. MRF and YSUNew check Lturb, Kmv, pblh and u_star.
  • YSUNew and YSU reuse the deck through RUNTIME_OPTIONS. YSU gets a cooled surface, since it aborts in unstable conditions. Legacy YSU never writes the surface-layer pblh, so its run leaves pblh out of the 2D plotfile and covers only the full-column assert (plus Lturb, Kmv, u_star).

Results

Built on development, macOS, AppleClang, Release with -DAMReX_BOUND_CHECK=ON -DERF_ENABLE_MPI=ON -DERF_ENABLE_ALL_WARNINGS=ON -DERF_ENABLE_TESTS=ON, 2 MPI ranks. The "development" column uses the three PBL sources from 489a2aa with the final deck and harness.

test	development PBL sources	this PR
ABL_MRF_Tiling	tiled run aborts at step 1: (0,8,0,0) is out of bound (-1:16,-1:7,0:0,0:0) on rank 0, (0,24,0,0) is out of bound (-1:16,15:23,0:0,0:0) on rank 1, in ComputeDiffusivityMRF	passed, every field bit-identical (all norms exactly 0)
ABL_YSUNew_Tiling	tiled run aborts at step 1 with the same two out-of-bound reads, in ComputeDiffusivityYSUNew	passed, bit-identical
ABL_YSU_Tiling	tiled run fails Assertion 'sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2)'	passed, bit-identical

Ranges of the fields the tests require to vary (untiled run, step 10):

	Lturb (m)	Kmv (m²/s)	pblh (m)	u_star (m/s)
MRF	168.0 – 203.9	0.10 – 19.1	168.0 – 203.9	0.457 – 0.487
YSUNew	180.3 – 200.7	0.10 – 18.4	180.3 – 200.7	0.457 – 0.487
YSU	152.2 – 181.9	0.034 – 2.25	not written	0.352 – 0.388

On development the untiled run of the same deck completes all 10 steps, so the failures come from tiling alone.

With the HOL_abs restoration, erf_exec builds with no warnings outside Submodules (development has one, the unused HOL_abs). On the MRF deck (heated, so the unstable arm is taken) the plotfiles with and without the restoration are bit-identical (fcompare, zero tolerance, 3D and 2D).

The same change was first made and tested in the hgopalan/ERF fork (hgopalan#393). There, a Release build of the old code (no bound checks) failed ABL_MRF_Tiling and ABL_YSUNew_Tiling on the 2D pblh: absolute difference 16 m, relative 1.0. The 3D fields were identical.

Review changes

• 53f5ad7 (from @asalmgren's review): register the tests only with ERF_ENABLE_MPI AND NOT WIN32; treat an empty -D value as missing in the script; remove the step-suffixed plotfiles (<run>_plt*) before each run; note that legacy YSU does not write pblh.
• Follow-up from auditing those findings: amrex_fextrema showed the original deck's pblh/Lturb were 16 m (the 0.5·dz floor) in every column, because with a 15 m/s wind and the inversion above the domain top the bulk-Richardson search never crossed Ribcr. The comparison was between two copies of one constant. The deck now puts the inversion inside the perturbed layer and the harness requires the named fields to vary; on the step-0 2D plotfile (pblh still 1e150) it stops with pblh is uniform (1e+150).

Not changed here

ApplyPBLHSmoothing (ERF_PBLModels.H, opt-in erf.enable_pblh_smoothing, default off) has a related problem. It smooths over the grown tile but clamps its stencil only at the domain edge. So it reads one column outside the array at every interior box or tile edge, even without tiling. With bound checks and several boxes, MRF and YSUNew abort there. A decomposition-independent fix needs neighbouring columns from outside the tile, so it is left for a separate PR (#3978).

🤖 Generated with Claude Code

[asalmgren]

Code review: correctness pass over the diff

What I verified as correct in the three PBL source changes (no findings against them):

• pblh_mf is built on eddyViscosity.boxArray() with 0 ghost cells, and the loop is an MFIter(..., TileNoZ()), whose tileboxes exactly partition the valid box. So switching the copy from mfi.validbox() to mfi.tilebox() in ERF_ComputeDiffusivityMRF.cpp and ERF_ComputeDiffusivityYSUNew.cpp both (a) removes the out-of-bounds read of the per-tile pbl_height_corrector FAB, which is allocated only on PerpendicularBox<ZDir>(mfi.growntilebox(IntVect(1,1,0))), and (b) removes a real OpenMP write race — every tile previously wrote the whole valid box, so with #pragma omp parallel multiple threads stored to the same pblh_mf cells. Coverage of pblh_mf is unchanged (union of tileboxes == valid box), so set_pblh's MultiFab::Copy still sees every valid cell.
• The copy point is still the right one: nothing rewrites pblh_corr_arr after it. MRF writes it only at lines 545/548, YSUNew only at 581/583 and 971/973 — both strictly before the new tbx loop.
• ERF_ComputeDiffusivityYSU.cpp: TileNoZ() does satisfy the full-column AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && ...), and the remaining x/y tiling is safe — growntilebox only grows where a tile abuts the valid-box boundary, so tiles never overlap and the K_turb writes stay race-free while the union still covers the full one-cell ghost ring exactly as the untiled path did. On GPU TileNoZ() returns the zero vector, which AMReX treats as no tiling (AMReX_MFIter.cpp:311), so behaviour there is identical to the old TilingIfNotGPU().
• The "16+ cells" wording in the new YSU comment is accurate, not off by a factor: AMReX's buildTileArray uses max(ncells/tilesize, 1), so a 9-cell-deep box stays one tile and 16 is indeed the first split.

All my findings are in the new test harness, posted inline, plus one pre-existing item below.

---

Source/PBL/ERF_PBLModels.H:517 — low (pre-existing; filed as #3978, not introduced by this PR).

Noting it here rather than inline because the file is not part of this diff. ApplyPBLHSmoothing allocates pblh_temp on xybx and clamps its 5-point stencil only against the domain edge:

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

xybx is the tile grown by one, so for any box or tile whose lo is interior to the domain, i == xybx.smallEnd(0) still satisfies i > dom_lo.x and the kernel reads pblh(i-1, j, 0) one column outside the FArrayBox. Concretely: with erf.enable_pblh_smoothing=true (default off) and more than one box in x or y, MRF and YSUNew abort under -DAMReX_BOUND_CHECK=ON and otherwise produce decomposition-dependent PBL heights. This matches the "Not changed here" section of the PR description. The new parity tests do not enable smoothing, so nothing in CI will catch a regression here — tracked in #3978 so the follow-up PR is not lost.

[asalmgren]

high — these three tests are registered unconditionally, unlike every other execute_process-driven test in this file (add_test_anelastic_wall_diffusion, add_test_cloud_chamber*), which all sit inside the if(ERF_ENABLE_MPI) block starting at line 491.

Concrete failure: .github/workflows/windows.yml configures with -DERF_ENABLE_MPI:BOOL=OFF -DERF_ENABLE_TESTS:BOOL=ON and then runs ctest -C Debug -L regression, which selects these by their regression label. With MPI off, MPIEXEC_EXECUTABLE and MPIEXEC_NUMPROC_FLAG are undefined, so add_test passes -DMPIEXEC= and -DMPIEXEC_NUMPROC_FLAG= (empty). In RunTilingParity.cmake those become empty elements of _mpi_run, and unquoted expansion in COMMAND ${_mpi_run} ... drops empty list elements — I confirmed with a standalone cmake -P that the list ;;1 expands to just 1, so execute_process tries to exec a program named 1 and returns "No such file or directory", tripping the FATAL_ERROR. All three tests fail on that job.

Independently, on Windows resolve_test_exe returns a glob (${CMAKE_BINARY_DIR}/Exec/*/erf_exec.exe) and FCOMPARE_EXE is also a glob (CMakeLists.txt:238). add_test_r gets away with those because it runs through sh -c, which expands them; execute_process does no globbing, so even with MPI enabled on Windows the exe and fcompare paths would not resolve.

Suggest gating the three registrations on ERF_ENABLE_MPI (and excluding WIN32), the same way the sibling parity tests are.

[hgopalan]

Fixed in 53f5ad7. The three registrations are now inside if(ERF_ENABLE_MPI AND NOT WIN32), with a comment giving both reasons: execute_process drops the empty MPIEXEC, and it does not expand the Windows erf_exec/fcompare globs.

Checked by configuring with -DERF_ENABLE_MPI=OFF -DERF_ENABLE_TESTS=ON: ctest -N -R Tiling lists 0 tests. The MPI build still registers all three, and they pass with every field bit-identical.

[asalmgren]

low — this case advertises a pblh parity check that cannot fail.

Only ComputeDiffusivityMRF and ComputeDiffusivityYSUNew call SurfLayer->set_pblh; legacy ComputeDiffusivityYSU never does. With erf.most.pblh_calc=YSU, SurfaceLayer::computes_pblh() is true, so ERF_Plotfile2D.cpp takes pblh_source = m_SurfaceLayer->get_pblh(lev) — but that MultiFab still holds the bogus_large_value (1e18) it was initialised with, in both the tiled and untiled run. The pblh column of the 2D comparison is therefore vacuous for this test (the other three 2D fields and all 3D fields are still meaningful, and the test does still cover the sbx full-column abort, which is the point of the case).

Worth a comment on the line so a future reader does not treat ABL_YSU_Tiling as pblh coverage.

[hgopalan]

Agreed, and confirmed: SurfaceLayer::update_pblh does nothing for the YSU type, and only MRF and YSUNew call set_pblh. The 2D pblh therefore stays at bogus_large_value in both runs.

53f5ad7 adds a comment on the registration. It says ABL_YSU_Tiling covers only the full-column assert, and that its pblh comparison is vacuous.

[asalmgren]

medium — if(NOT DEFINED ${_required}) is satisfied by an empty value, so this guard does not catch the case it exists for.

add_test_tiling_parity always emits -DMPIEXEC=${MPIEXEC_EXECUTABLE} and -DFCOMPARE=${FCOMPARE_EXE}; when those CMake variables are unset the arguments still arrive as -DMPIEXEC= / -DFCOMPARE=, which defines the script variable as the empty string. DEFINED then returns true, the intended FATAL_ERROR "RunTilingParity.cmake requires MPIEXEC" never fires, and the failure surfaces much later as an opaque execute_process error (see the non-MPI Windows scenario on the registration line).

if(NOT DEFINED ${_required} OR "${${_required}}" STREQUAL "")

[hgopalan]

Fixed in 53f5ad7 with your condition, plus a comment explaining why an empty value has to be rejected too. cmake -DMPIEXEC= -DMPIEXEC_NUMPROC_FLAG= ... -P Tests/RunTilingParity.cmake now stops straight away with RunTilingParity.cmake requires MPIEXEC.

[asalmgren]

low — this REMOVE_RECURSE never matches anything, so the intended cleanup does not happen.

AMReX names plotfiles with the step suffix, so the directories actually written are tiled_plt00000, tiled_plt00010, tiled_plt2d00010, … — there is never a directory called tiled_plt or tiled_plt2d. Stale plotfiles from an earlier configuration of the deck therefore survive into the next run. Since PLTFILE/PLT2DFILE are hard-coded to 00010, a later edit to max_step or erf.plot_int_1 that stops producing a step-10 plotfile would leave both runs' step-10 directories from a previous build in place, and fcompare would compare two stale trees and report a pass.

file(GLOB _stale "${WORKING_DIRECTORY}/${_run}_plt*")
if(_stale)
  file(REMOVE_RECURSE ${_stale})
endif()

[hgopalan]

Fixed in 53f5ad7. Each run now globs <run>_plt* and removes what matches, so the step-suffixed 3D and 2D plotfiles are both covered.

Checked with a stale tiled_plt00010 and tiled_plt2d00010 in the working directory: both were gone before the run started.

[hgopalan]

Pushed 32d86ec, a follow-up from auditing the review findings above.

Running amrex_fextrema on the untiled plotfiles showed that the MRF and YSUNew pblh (2D) and Lturb (3D) were 16 m in every column. That is the floor 0.5*dz (dz is 32 m in this deck): with a 15 m/s wind and the inversion above the domain top, the bulk-Richardson search never crossed Ribcr. So the parity check on pblh was comparing one constant with itself, which catches an out-of-bounds read but not a tile bug that permutes columns.

• The deck now has a 6 K inversion at 150–250 m, a 5 m/s wind and local MOST averaging. pblh ranges over 168–204 m (MRF) and 180–201 m (YSUNew), and u_star, t_star, Olen differ between columns too.
• RunTilingParity.cmake takes VARYING_3D / VARYING_2D lists and fails unless amrex_fextrema reports min < max for each named field in the untiled run. On the step-0 2D plotfile (pblh still 1e150) it stops with pblh is uniform (1e+150).
• Legacy YSU drops pblh from erf.plot2d_vars_1 instead of comparing the constant, and its cooling is -0.02 (with the 5 m/s wind, -0.05 stopped the MOST iteration converging).

Re-checked with the three PBL sources from 489a2aa and this deck: the tiled MRF and YSUNew runs abort at step 1 ((0,8,0,0) is out of bound (-1:16,-1:7,0:0,0:0) on rank 0, (0,24,0,0) ... (-1:16,15:23,0:0,0:0) on rank 1), YSU fails its full-column assert; with the fix all three pass bit-identical. The PR description is updated with the field ranges.