Fix momenta ghost faces and planar surface-layer arrays on z-split BoxArrays (terrain-fitted meshes) - #3970
Fix momenta ghost faces and planar surface-layer arrays on z-split BoxArrays (terrain-fitted meshes)#3970hgopalan wants to merge 10 commits into
Conversation
…ZSplit On a terrain-fitted mesh whose BoxArray is split in z (more cells in z than amr.max_grid_size) two things went wrong, independently of the closure, the boundaries, the terrain source, stretching or the solver: 1. project_initial_velocity converted velocity to momentum on the valid faces only, and project_momenta's OmegaFromW averages (rho0 u, rho0 v) across each w-face, reading one z-ghost face at the top and bottom face of every box; at a box face inside the domain that ghost still held the bogus_large_value placeholder, so the divergence handed to the solver was of order 1e139 (the domain top reads the same placeholder times an exactly zero slope, which is why single-box columns never showed it). Fill the momenta ghost faces after the conversion. 2. The planar surface-layer arrays (u*, t*, q*, L, t_surf, q_surf, pblh, diagnostic source) and the MOST region averages hold one 2D box per 3D box, so the split gave duplicate 2D boxes of which only the surface copy is computed, and FillBoundary filled the surface copy's ghost cells from the uncomputed duplicate: 1e150 in the surface stress at box-boundary faces, -inf momenta, then 0 * inf in WFromOmega and an invalid-operation trap in ERFPhysBCFunct_w during the first advance. New FillPlanarBoundary (Source/Utils/ERF_PlanarBoundary) gathers the surface copies on a BoxArray without duplicates and fills every copy from that; plain FillBoundary is kept when there are no duplicates or with EB terrain. Both are silent under amrex.init_snan with the trap armed because the leaked value is the deterministic placeholder. Release-active asserts now abort if a placeholder would be used at either site. ABL_MOST_WOA_ZSplit: anelastic neutral flow over a radial Witch of Agnesi hill on a 32 x 32 x 40 fitted mesh with amr.max_grid_size = 32, MOST and Smagorinsky, 10 steps. On development it traps in the first advance; the gold plotfile is from this branch. Existing MOST and terrain regression tests are unchanged. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
|
@asalmgren @AMLattanzi I do not understand the bug per-se. I was trying to clean-up the one equation RANS model which Eliot put added and this bug came along the way. |
Review: ERF PR #3970Substantive1. The
|
…ed planar fill, unit test) Momenta ghost faces. The FillBoundary in project_initial_velocity ran before project_momenta's ConvertForProjection, which rescales (rho u) to (rho0 u) on the valid faces only, so with erf.fixed_density = 0 the ghost faces OmegaFromW reads at a split kept (rho u) while the valid faces held (rho0 u). The fill now sits in project_momenta, after the conversion and just before OmegaFromW, for every caller, and covers only what OmegaFromW reads: one z-ghost layer of x- and y-momentum. The per-step release-active reduce that checked those faces is gone. The comment now gives the real reasons the domain boundaries are safe: the w-face k = 0 is set to zero and nothing below it is read, and the top ghost face holds the momentum VelocityToMomentum extrapolates there (an exactly zero metric slope was not the reason). Planar fill. MakeSurfaceBoxes/FillPlanarBoundary become a PlanarBoundary class that keeps a gather buffer per index type and component count, so a fill no longer allocates a MultiFab each call. It takes the planar BoxArray itself, so the surface copies are the planar boxes as built rather than a second collapse to k = 0. The header states that the fill also overwrites the valid region of the uncomputed copies, how shared faces of face-centered averages are handled, and that a level with no box at the surface is left as it is. w_star is now filled with the other surface-layer fields when it is computed (erf.most.include_wstar). The doxygen blocks the new declarations had displaced are back on their functions. Tests. ABL_MOST_WOA_ZSplit runs with amr.max_grid_size = 16 (4 columns of 3 boxes, split in z at k = 16 and 32; the old comment's k = 32 split was wrong, it was 20 + 20), so the surface fill crosses box and rank boundaries in the plane; gold regenerated. New parallel gtest PlanarBoundaryParallel covers the surface-copy selection, cell- and face-centered fields with 1 and 3 components, a refill with new values, a BoxArray without a split, and a level above the surface. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
rebalance_columns integrated every box of the level from its own lowest cell. A box stacked on another box of the level (a BoxArray split in z, e.g. amr.max_grid_size below n_cell z) therefore started from its lowest cell's density as interpolated from the sounding, not from the rebalanced column below, and the hydrostatic base state depended on the split. On ABL_MOST_WOA_ZSplit (32x32x40 hill, input_sounding) at step 0, one box against 12 boxes (split at k = 16 and 32): dens_hse and density first differ at k = 16 by 2.1e-7 and again at k = 32, pres_hse by 0.025 Pa; every cell below k = 16 is bitwise identical. With this change both are bitwise identical in every layout, and a single-box run is unchanged. The boxes are now integrated in bands of equal lowest index, bottom up. The state each column reaches in a cell (pressure, theta, qv, qt and the cell-centre height) is kept in a MultiFab whose z ghost cells are filled from the bands already done, and a box continues from the cell below it where that cell belongs to the level. A column with nothing of the level below it (the domain bottom, a refined patch aloft) starts from its own lowest cell as before. With a single band nothing is kept or filled. The use_sfc assertion now accepts a split column whose cells below are covered. New parallel gtest RebalanceColumnsParallel: 4-column, 3-box, uneven and staggered-footprint layouts match one box (maintain_Th true and false, with and without use_sfc), and so does a patch aloft. On the old code the split cases fail and the use_sfc cases abort. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
After the terrain Poisson solve, project_momenta adds the fluxes to the valid faces of (rho0 u, rho0 v, Omega) and converts Omega back to rho0 w with WFromOmega, which averages (rho0 u) and (rho0 v) over the faces below and above each w-face. At the lowest and highest w-face of a box inside the domain it reads a z-ghost face, and that face still held the momentum from before the projection (filled for OmegaFromW). With a BoxArray split in z, w at every split face came out wrong. ABL_MOST_WOA_ZSplit, one box against 12 boxes on one rank: after the initial projection w differed by 0.0093 m/s at the faces k = 16 and 32 only, with identical GMRES residual histories; after 10 steps u differed by 0.0091 m/s and theta by 0.18 K. With this fill and the previous commit, every split layout (12 boxes on 1 and 2 ranks, 2 boxes on 2 ranks) agrees with one box to round-off (u 6.6e-13, theta 2.8e-13), as the x/y-only split already did. The fill alone leaves u 9.7e-7 and density 2.1e-7 from the base state. Single-box runs are bitwise unchanged. The ABL_MOST_WOA_ZSplit gold held the split-dependent result (u 0.011 and theta 0.18 K from the new one) and is regenerated. A new box-parity CTest, ABL_MOST_WOA_ZSplit_BoxParity (add_test_box_parity, RunBoxParity.cmake), runs the deck on one box against its 12 boxes and compares the plotfiles at a relative tolerance of 1e-9; the unfixed binary fails it with the numbers above. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
|
Thanks, @AMLattanzi, this was a useful review. Everything is addressed in 3202eaa. Two further commits, 6d5a6bd and de4ecf3, fix a split dependence that came up while checking the result. Point by point: 1. Fill before 2. Release-active assert on the per-step path. Removed. With the fill inside 3. Why the domain boundaries are safe. You're right; the zero-slope explanation was wrong. The comment now says that the w-face k = 0 is set to zero and nothing below it is read, and that the top ghost face holds the momentum One related thing I left alone:
Face-centered averages. The header notes that neighbouring surface boxes share a face in the buffer, and that ParallelCopy picks which box supplies it, as FillBoundary did; both should hold the same value. The index-order assumption for the xvel/yvel BoxArrays is still caught by the box-equality assert, which now also checks the owning rank.
Minor:
Beyond the review: the result still depended on the split. Comparing the deck on one box with its split layouts showed two more defects:
With both commits, every split layout agrees with one box to about 6e-13 (density identical), and single-box runs are bitwise unchanged. The All 76 tests of the |
erf_init_dens_hse_dry (terrain path) and erf_enforce_hse integrated every box of a level from the ghost cell below it. On level 0 that never matters: initHSE copies a z-split base state onto full-height ChopGrids2D columns first. On a refined level it does: a box stacked on another fine box started from the density and pressure interpolated from the coarse level, not from the column below, so the fine base state depended on where the fine boxes were split. 32x32x40 WoA hill, Isentropic init, level 1 over the centre (ref 2 2 1, terrain_smoothing 2 so both runs have the same fine mesh), fine level as full columns against split at k = 16 and 32: at step 0 dens_hse differed from k = 16 up by 3.9e-5 (3.6e-7 from k = 32), pres_hse by 4.7 Pa, and still did after 5 steps. With this change the split fine level is bitwise identical at step 0 and within 1e-15 (density) and 2e-13 (theta) after 5 steps. Level 0 split against one box, and every layout with a single band, are bitwise unchanged; Terrain2Lev_STF_interp and _transform have one band per level on 1 and 2 ranks, so their gold stands. Both loops now visit the boxes in bands of equal lowest z index, bottom up. Before each band but the first, fill_below_band (new, Utils/ERF_ColumnBands) fills the cells just below the band, widened laterally by one cell, from the bands already done: the lateral ghost cells of the boxes below first, then their valid cells, which take precedence. A cell with nothing of the level below it (under a patch aloft) keeps its interpolated value. The lateral ghost cells matter because FillPatchTwoLevels adds the fine base state into the coarse-fine ghost cells of the state. New parallel gtest ColumnBandsParallel: the fill precedence, and erf_init_dens_hse_dry on a full periodic domain, a patch and a patch aloft in five split layouts, compared exactly with one box on valid and lateral ghost cells. Without the fill all four tests fail; with only the valid-cell copy (what a FillBoundary would give) the two patch tests fail. Not addressed: with terrain_smoothing = 0 (BTF) the fine mesh itself depends on the z split (z_phys differs by 0.54 m above k = 16). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… runner CI on every platform that runs ctest failed ABL_MOST_WOA_ZSplit and its box-parity test: the deck is anelastic on a terrain-fitted mesh, so its projection is the FFT-preconditioned GMRES, and no CI job builds with ERF_ENABLE_FFT (amrex::Abort "Rebuild with USE_FFT = TRUE"). On WIN32 the box-parity runner also received the unexpanded exe glob, since it was registered outside the MPI block its sibling add_test_cloud_chamber_parity lives in. - The anelastic tests are registered only with ERF_ENABLE_FFT, and the cmake -P box-parity tests only with ERF_ENABLE_MPI AND NOT WIN32. - New ABL_MOST_WOA_ZSplit_NoSub: the same hill, MOST surface layer, sounding and 4 x 3 box split, run compressible with erf.substepping_type = None and fixed_dt = 0.02 (the acoustic limit), so it runs in every build. It covers the planar surface-layer fill and rebalance_columns; the acoustic substep is off because it has its own per-box column assumption. The e279ca5 binary aborts on the split at start-up ("Base state is too far out of HSE", extrema at the split faces k = 16 and 32) and runs on one box; with this branch the split agrees with one box to 1e-14 after 10 steps, over which u moves by 0.76 m/s and w by 1.9 m/s. Gold and a box-parity test added. - RunBoxParity.cmake tests its arguments for a value rather than DEFINED (-DX= defines X as empty), requires MPIEXEC_NUMPROC_FLAG only with MPIEXEC, and fails unless the split run has more level-0 boxes than the reference, so a deck change cannot turn the comparison trivial. - rebalance_columns uses column_bands() from ERF_ColumnBands instead of its own copy; results are bitwise unchanged. Checked: configure with MPI off registers only ABL_MOST_WOA_ZSplit_NoSub; with MPI on and FFT off, _NoSub and _NoSub_BoxParity only. The runner fails as intended with the split option removed, an empty TEST_EXE, a missing numproc flag, and the e279ca5 binary. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
|
The CI reds on this PR were both mine, in the tests rather than the fixes, and are addressed in d3e29c9: The two new tests aborted on every platform that runs ctest.
Checks, in the spirit of the #3972 review: configuring with MPI off registers only I also went through the rest of |
…; correct the dz != dx notes The terrain_height wall distance read z_nd(:,:,klo) from every box, but boxes above a z split of the BoxArray do not hold the surface nodes, so their distance was garbage (778 m error on the 3D hill deck at amr.max_grid_size_z = 10, clamped to epsilon). Gather the surface slab onto every box with a ParallelCopy first, and guard the gathered slab with a reduction over its own ghost width (a scalar-ghost MultiFab::max would read past the slab, which has no z ghosts). With the fix the split and unsplit layouts give bit-identical wall distances, and the unsplit distance is unchanged. The 1.788e139 pre-projection divergence recorded in PLAN.md phase 6, the 3D hill README, RESULTS.md, the deck and the Poisson solver comment as a dz != dx failure was the same z-split layout (unfilled momenta ghost faces in the initial projection and duplicated planar surface-layer arrays, fixed in erf-model#3970), not the aspect ratio; the text now says so. The Poisson wall-distance multigrid does still diverge at dx = 2 dz. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
3 participants
Summary
On a terrain-fitted mesh whose BoxArray is split in z (more cells in z than
amr.max_grid_size, so a column is more than one box) five things went wrong. The first two made the run blow up. The other three made a run that did complete depend on where the columns were split.Projection, before the solve.
VelocityToMomentumandConvertForProjectionwrite the valid momentum faces only.project_momentathen convertsrho0 wto Omega withOmegaFromW, which averagesrho0 uandrho0 vover the faces below and above each w-face. At the lowest and highest w-face of a box it therefore reads one face in the z-ghost layer. Where that box face lies inside the domain, the ghost face heldbogus_large_value(1e150) from allocation in the initial projection, and the divergence handed to the solver was of order 1e139. The domain boundaries never showed it, for other reasons: the w-face k = 0 is set to zero and nothing below it is read, and the top ghost face holds the momentumVelocityToMomentumextrapolates there. Fix:project_momentafills the z-ghost layer of x- and y-momentum after therho0conversion, just beforeOmegaFromW, for every caller. In the time stepapply_bcshas already filled those faces, so this fill changes no value there.Surface layer and MOST averages. The planar 2D MultiFabs are built on the z-collapse of the 3D BoxArray, one 2D box per 3D box. They are
u_star,w_starwhen computed,t_star,q_star,olen,t_surf,q_surf,pblh, the diagnostic source, and the region averages inMOSTAverage.FillBoundaryon such a BoxArray is ill-defined: it fills the surface copy's ghost cells from the uncomputed duplicate. That puts 1e150 into the surface stress at box-boundary faces, then-infmomenta, then0 * infinWFromOmega, and an invalid-operation trap inERFPhysBCFunct_wduring the first advance.PlanarBoundary(Source/Utils/ERF_PlanarBoundary.{H,cpp}) records the surface copies once per level. It gathers them into a cached buffer without duplicates, one per index type and component count, and fills every copy from that withParallelCopy, valid region and ghost cells alike. UnlikeFillBoundaryit also overwrites the valid region of the uncomputed copies, so any copy can be read afterwards. PlainFillBoundaryis kept when there are no duplicates, and with EB terrain (3D fields).Column rebalance (
input_sounding).rebalance_columnsintegrated every box of the level from its own lowest cell. A box stacked on another box therefore started from its lowest cell's density as interpolated from the sounding, not from the rebalanced column below.dens_hseand the density first differed at the split (k = 16) by 2.1e-7, andpres_hseby 0.025 Pa.use_sfcassertion now accepts a split column whose cells below are covered.Projection, after the solve.
WFromOmegareads the same z-ghost faces asOmegaFromW, but the fluxes are added to the valid faces only, so those ghost faces still held the momentum from before the projection. GMRES residual histories were identical between layouts; only w at the split faces differed, by 0.0093 m/s after the initial projection. Fix: fill the z-ghost layer of x- and y-momentum again beforeWFromOmega.Fine-level hydrostatic initialisation over terrain.
erf_init_dens_hse_dry(terrain path) anderf_enforce_hseintegrated every box of a level from the ghost cell below it.initHSEcopies a z-split base state onto full-height columns first.dens_hsediffered from k = 16 up by 3.9e-5, andpres_hseby 4.7 Pa.fill_below_band(Source/Utils/ERF_ColumnBands.{H,cpp}) fills the cells just below it, widened laterally by one cell, from the bands already done. It copies the lateral ghost cells of the boxes below first, then their valid cells, which take precedence. The lateral ghost cells matter becauseFillPatchTwoLevelsadds the fine base state into the coarse-fine ghost cells.Items 1 and 2 are silent under
amrex.init_snan = 1with the invalid-operation trap armed, because the leaked value is the deterministic placeholder, not a NaN.PlanarBoundaryaborts with a message if the surface boxes overlap in the plane, or if a surface box does not match its planar box.Numbers
Mesh that exposed items 1 and 2: a periodic 64 x 64 x 40 flat fitted mesh (dx = dy = 40 m, dz = 20 m), uniform 10 m/s wind, anelastic. With
amr.max_grid_size = 32, each column is boxes k = 0..23 and k = 24..39.developmentu_starghost cells at box-boundary columns afterupdate_fluxesamrex.fpe_trap_invalid = 1ERFPhysBCFunct_wItems 3 and 4 on
ABL_MOST_WOA_ZSplit(below). The reference is the same deck on one box on one rank (amr.max_grid_size = 64). The table gives the max-norm difference after 10 steps;developmentcannot run the split layouts at all.Single-box runs are bitwise unchanged by all five items. So is the 64 x 64 x 20 single-box-per-column mesh above: 40 steps, all 12 plotfile variables at 0 difference against
development.Tests
ABL_MOST_WOA_ZSplit(regression, FFT builds only).amr.max_grid_size = 16gives 4 columns of 3 boxes, split in z at k = 16 and 32. The surface fill therefore crosses box boundaries in the plane, and rank boundaries on 2 ranks.development. There it trips the invalid-operation trap in the first advance. With the trap disabled, the first step's projection sees a NaN divergence and the step aborts on a negative theta (RhoTheta is negative ... -inf).ERF_ENABLE_FFT, so this test and its box-parity test are registered only underif(ERF_ENABLE_FFT). They run in FFT builds locally; the first push of this PR had them unguarded, and every CI platform that runs ctest failed them withRebuild with USE_FFT = TRUE.ABL_MOST_WOA_ZSplit_NoSub(regression, every build). The same hill, surface layer, sounding and 4 x 3 split, run compressible witherf.substepping_type = Noneandfixed_dt = 0.02(the acoustic limit). It covers items 2 and 3; the acoustic substep stays off because it has its own per-box column assumption (Dirichlet rows at every box's lowest and highest face), which is not addressed here. Over its 10 steps u moves by 0.76 m/s and w by 1.9 m/s, so the comparison is on a live field. Ondevelopmentand on the first commit of this PR the split run aborts at start-up withBase state is too far out of HSE, with the extrema at the split faces k = 16 and 32 (item 3), while one box runs; with this branch the split agrees with one box to 1e-14. Gold from this branch, 2.7 MB. On Windows (no MPI) it runs on one rank asadd_test_rtests do.ABL_MOST_WOA_ZSplit_BoxParityand_NoSub_BoxParity(regression, box-parity label). Newadd_test_box_parity/Tests/RunBoxParity.cmake. Each runs its deck on one box on one rank and on its 12 boxes on 2 ranks, then compares the plotfiles at a relative tolerance of 1e-9, and fails unless the split run has more level-0 boxes than the reference, so a deck change cannot make the comparison trivial. They take about 2 s each. Without items 3 and 4 the anelastic one fails with the numbers above; the_NoSubone fails on the first commit of this PR because its split run aborts. The runner is acmake -Pscript, so both are registered underif(ERF_ENABLE_MPI AND NOT WIN32); it checks its arguments for a value rather thanDEFINED(-DX=defines X as empty). Checked by configuring with MPI off (only_NoSubis registered) and with MPI on but FFT off (_NoSuband_NoSub_BoxParity), and by running the script with the split option removed, an emptyTEST_EXE, a missing numproc flag and the old binary, each of which fails with its own message.PlanarBoundaryParallel(parallel gtest). It checks that:FillBoundary;It passes on 1 and 2 ranks, and fails when
fillis replaced by a plainFillBoundary: the 1e30 placeholder is left in the uncomputed copies.RebalanceColumnsParallel(parallel gtest). (rebalance_columnssharescolumn_bands()with item 5; results are bitwise unchanged.) Layouts of 4 columns, 3 boxes per column, uneven boxes and staggered footprints all match one box. This holds withmaintain_Thtrue and false, with and withoutuse_sfc, and for a patch aloft. On the old code the split cases fail and theuse_sfccases abort.ColumnBandsParallel(parallel gtest). It covers:fill_below_band;erf_init_dens_hse_dryon a full periodic domain, a refined patch and a patch aloft, in five split layouts, compared exactly with one box on valid and lateral ghost cells.Without the fill all four tests fail. With only the valid-cell copy (what a
FillBoundarywould give), the two patch tests fail.Full runs (Release, MPI, FFT on, macOS arm64, on the branch merged with
developmentat dab7826):regressionlabel pass on 2 ranks, including the four ZSplit tests, BTF fine level: make the terrain-fitted mesh independent of the z split #3976'sTerrain2Lev_BTF_ZSplit, and the existing MOST and terrain cases (ABL_MOST,ABL_MOST_IMP_DIFF_WOA,MovingTerrain_nosub, and the two-levelTerrain2Lev_STF_interp/_transform, which have one band per level) against their unchanged gold files;Not run locally: Debug builds, single precision, GPU backends and Windows; CI covers those, and the first CI round showed the FFT and Windows-glob failures described above, both now fixed by the registration guards and the
_NoSubvariant.Item 5 changes no level-0 result: the
ABL_MOST_WOA_ZSplitlayouts are bit-identical with and without it.Notes
Found while running a 3D hill deck at dz = dx/2 in the ERF-Hazard RANS work, where the failure was first attributed to the aspect ratio. The same class of bug is worth keeping in mind for any new per-column operation or planar field: per-column data or integrations indexed by the 3D MFIter. The new box-parity test is a cheap way to catch the next one.
Same class, not changed here. The remaining per-box column integrations on
developmentfail loudly rather than silently:init_base_state_from_wrfinput, the NetCDF initialisation and the column diagnostics inERF_Deriveassert that boxes reach the bottom of the domain. Two are silent: the compressible acoustic substep (MakeFastCoeffsand the w solve set Dirichlet rows at every box's lowest and highest face) and the implicit vertical diffusion (vert_implicit_fac, box-localklo/khi), which is why the_NoSubdeck runs without substepping and every ZSplit test passeserf.vert_implicit=false.Not changed here:
ConvertForProjectionalso skips the momentum faces in the ghost layer above the domain top. So witherf.fixed_density = 0, the extrapolated top ghost face keepsrho uwhile the valid faces holdrho0 u. That is ondevelopmentalready, independent of the box split.terrain_smoothing = 0(BTF), the fine-level mesh itself still depends on the z split:z_physdiffers by 0.54 m above k = 16. That is being looked at separately.🤖 Generated with Claude Code