Canonical RANS CTests: run without FFT in CI; terrain-height wall distance on z-split layouts

Docs/sphinx_doc/RegressionTests.rst

@@ -560,7 +560,10 @@ dt = 10 s with the implicit column solve of momentum, theta and KE),
 ``RANS_Stable_ABL_Flat``, ``RANS_Convective_ABL_Flat``,
 ``RANS_Neutral_Hill_2D`` (and ``_Poisson``),
 ``RANS_Flat_Fitted_2D`` (and ``_Poisson``), ``RANS_Neutral_Hill_3D`` (and
-``_Poisson``); label ``rans``.
+``_Poisson``); label ``rans``. The flat cases run with the MLMG projection
+(``erf.use_fft=false``); the terrain-fitted cases need the FFT-preconditioned
+projection and are registered only when the build enables FFT
+(``ERF_ENABLE_FFT``).
 
 A sixth case, ``Timestep_Limits``, measures the largest stable time step of
 the vertical eddy diffusion on a neutral 4 x 4 x 200 column (dx = 800 m,

Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/README.md

@@ -18,24 +18,20 @@ the crest speed-up against the axisymmetric estimate of about 1.6 h/L
 | physics run | 4 h (9600 steps), about 10 min on 2 ranks |
 | smoke run | 40 steps (`ctest -R RANS_Neutral_Hill_3D`) |
 
-The mesh has unit aspect ratio on purpose. With dz different from dx
-(20 or 80 m at dx = 40 m) this deck, and the same mesh flattened, produce
-a divergence of order 1e139 before the very first projection, the wall
-distance Poisson solve then diverges, and the run aborts. The failure is
-not the aspect ratio as such (the Askervein deck runs at dx/dz of 0.5, 1
-and 2), not the lateral boundaries, the terrain source, grid stretching
-or the box layout, and not the RANS code (it appears in the initial
-projection, before any turbulence call, and with the closure switched to
-Smagorinsky). It is deterministic: under `amrex.init_snan = 1` with the
-invalid-operation trap armed, initialisation completes without a trap
-and the same 1.788e139 divergence, so it is an arithmetic error tied to
-dz relative to dx, not a memory read; a second, separate uninitialised
-read in the w boundary fill (`ERFPhysBCFunct_w`) then trips the trap in
-the first advance. Reproducer: `inputs_hill3d amr.n_cell="64 64 40"
-prob.hmax=1e-6 max_step=0 erf.mg_v=2 erf.v=1`. The
-wall distance is `erf.wall_dist_type = terrain_height` (exact to 1e-10 on
-a flat mesh, 0.01 % mean error on this hill); the Poisson distance is
-exercised by the `_Poisson` CTest variant.
+The mesh has unit aspect ratio. The 1.788e139 divergence before the
+first projection that this deck showed with dz different from dx
+(`amr.n_cell = 64 64 40` at `amr.max_grid_size = 32`) was not the aspect
+ratio: that mesh splits the BoxArray in z, and the initial projection
+read the unfilled momenta ghost faces at the internal box faces while
+the planar surface-layer arrays were duplicated across the stacked boxes
+(fixed in erf-model/ERF#3970). What does depend on dz relative to dx is
+the Poisson wall-distance solve, whose multigrid diverges at dx = 2 dz
+with or without the split, so the `_Poisson` CTest variant needs dz = dx;
+the `terrain_height` distance the deck uses has no such limit and is
+gathered from the surface boxes, so a z-split layout is fine with it.
+The wall distance is `erf.wall_dist_type = terrain_height` (exact to
+1e-10 on a flat mesh, 0.01 % mean error on this hill); the Poisson
+distance is exercised by the `_Poisson` CTest variant.
 
 ## Running
 

Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/inputs_hill3d

@@ -17,9 +17,10 @@ amrex.fpe_trap_invalid = 0
 fabarray.mfiter_tile_size = 1024 1024 1024
 
 # PROBLEM SIZE & GEOMETRY
-# Unit aspect ratio: on a 3D fitted mesh both the projection and the
-# wall-distance Poisson solves (MLMG) fail at dx = 2 dz (the 2D ridge and
-# Askervein, dx = dz, are fine), so dz = dx here.
+# Unit aspect ratio: the Poisson wall-distance solve (the _Poisson CTest
+# variant) diverges at dx = 2 dz; the terrain_height default has no such
+# limit. (The projection failure once seen at dz != dx on this deck was a
+# BoxArray split in z, fixed in erf-model/ERF#3970, not the aspect ratio.)
 geometry.prob_lo     = -2560. -2560.   0.
 geometry.prob_hi     =  2560.  2560. 800.
 amr.n_cell           =    64     64   20     # dx = dy = dz = 40 m

Exec/CanonicalTests/Canonical_RANS/PLAN.md

@@ -229,26 +229,26 @@ hill decks use it; the Poisson path keeps three `_Poisson` CTest
 variants. The Poisson solve is posed positive definite now (same
 iterates as before).
 
-Finding, outside this plan: on a 3D terrain-fitted mesh with dz different
-from dx (20 or 80 m at dx = 40 m, flat or hill, periodic or inflow,
-custom or file terrain, stretched or not, any box layout) the divergence
-of the initial field is of order 1e139 before the first projection, the
-wall-distance MLMG then diverges (residual 18x after one cycle, 1e10 by
-iteration 100, unchanged by the sign convention or by semi-coarsening),
-and the run aborts within a step. Askervein is unaffected at dx/dz of
-0.5, 1 and 2, so the aspect ratio is not the cause; the RANS code is not
-involved (same with Smagorinsky). Two separate defects: (1) the
-initialisation one is deterministic, since under `amrex.init_snan = 1`
-with the invalid-operation trap armed initialisation completes with no
-trap and the bit-identical 1.788e139 divergence, so it is an arithmetic
-error that depends on dz relative to dx, not a memory read; (2) the trap
-then fires in `ERFPhysBCFunct_w` during the first advance, an
-uninitialised read in the w boundary fill that the unit-aspect mesh does
-not trigger. Reproducer for (1): `inputs_hill3d amr.n_cell="64 64 40"
-prob.hmax=1e-6 max_step=0 erf.mg_v=2 erf.v=1` and read the divergence
-before the solve. Harish's rule for the decks: mass inflow and
-pressure outflow instead of periodic if periodic turns out to be the
-issue; it did not, so the hill decks stay periodic at unit aspect ratio.
+Finding, outside this plan (root cause found after this phase): with
+`amr.n_cell = 64 64 40` at `amr.max_grid_size = 32` the initial field's
+divergence was 1.788e139 before the first projection and the first
+advance tripped the invalid-operation trap in `ERFPhysBCFunct_w`. It
+was never the aspect ratio: that mesh splits the BoxArray in z, and (1)
+`project_initial_velocity` converted velocities to momenta on the valid
+faces only while `OmegaFromW` read one z-ghost face at every box face,
+which held the 1e150 allocation placeholder at the internal faces, and
+(2) the planar surface-layer and MOST arrays hold one 2D box per 3D box,
+so the split duplicated them and `FillBoundary` copied the uncomputed
+duplicate into the surface copy. Both are fixed in erf-model/ERF#3970
+(momenta ghost fill in `project_momenta`, `FillPlanarBoundary`). A third
+instance was in this branch: the `terrain_height` wall distance read the
+surface nodes `z_nd(:,:,klo)` from boxes above the split, which do not
+hold them (220 m error); it now gathers the surface slab onto every box.
+What does depend on dz relative to dx is the wall-distance MLMG solve
+(residual 18x after one cycle, 1e10 by iteration 100, unchanged by the
+sign convention or semi-coarsening, with or without the split), so the
+`_Poisson` variants keep dz = dx and `terrain_height` is the default
+for the hill decks.
 
 ## Phase 7: documentation and diagnostics
 

Exec/CanonicalTests/Canonical_RANS/README.md

@@ -22,8 +22,10 @@ record is `PLAN.md` and the numbers per phase are in `RESULTS.md`.
   (`erf_plotfile.py`, standard library only) and prints one row per check:
   measured value, target, tolerance, pass. Its exit code is the verdict.
   `--smoke` runs the structural checks that must hold after a few steps
-  (the CTest entries, `ctest -L rans`); `--physics` adds the checks that
-  need the converged run in the table above.
+  (the CTest entries, `ctest -L rans`; the flat entries pass
+  `erf.use_fft=false` and the terrain entries exist only in FFT builds);
+  `--physics` adds the checks that need the converged run in the table
+  above.
 - `Timestep_Limits` has no `check_<case>.py`: its `sweep_dt.py` runs ERF
   itself (spin-up, then a restart per time step) and applies the checks
   listed in its README (`ctest -L dt_sweep`, not part of `regression`).

Exec/CanonicalTests/Canonical_RANS/RESULTS.md

@@ -177,12 +177,13 @@ every field identical to 1e-14.
 Askervein (20 steps, 4 ranks, 26 s): walldist 7.9 to 707 m, KE up to
 6.9 m2/s2, Kmv up to 9.4 kg/m/s, all finite.
 
-Mesh finding (not a RANS matter, see PLAN phase 6): any 3D fitted mesh
-here with dz != dx yields a deterministic pre-projection divergence of
-1.788e139 and aborts, independent of the closure and not a memory read
-(no trap under `amrex.init_snan` at initialisation); a separate
-uninitialised read in the w boundary fill trips the trap in the first
-advance.
+Mesh finding (not a RANS matter, see PLAN phase 6): the 1.788e139
+pre-projection divergence first seen at dz != dx was a BoxArray split in
+z (unfilled momenta ghost faces in the initial projection and duplicated
+planar surface-layer arrays), fixed in erf-model/ERF#3970; the
+`terrain_height` wall distance had the same class of read and now
+gathers the surface nodes onto every box. The Poisson wall-distance
+multigrid does diverge at dx = 2 dz, with or without the split.
 
 
 ## Phase 9: implicit vertical diffusion of scalars under anelastic

Source/LinearSolvers/ERF_PoissonWallDist.cpp

@@ -90,11 +90,36 @@ void ERF::poisson_wall_dist (int lev)
         // terrain height used by the amr-wind immersed terrain and Kynema).
         Print() << "Calculating wall distance from the terrain height (normal-projected)" << std::endl;
         const int klo = geomdata.Domain().smallEnd(2);
+
+        // The surface nodes z_nd(:,:,klo) live only in the boxes that touch the surface;
+        // when the BoxArray is split in z the boxes above hold nodes from their own k
+        // range, so gather the surface slab onto every box (the same 2D footprint, at klo).
+        BoxList bl_surf = z_phys_nd[lev]->boxArray().boxList();
+        for (auto& b : bl_surf) { b.setRange(2,klo); }
+        BoxArray ba_surf(std::move(bl_surf));
+        IntVect ng_surf = z_phys_nd[lev]->nGrowVect(); ng_surf[2] = 0;
+        MultiFab znd_surf(ba_surf, z_phys_nd[lev]->DistributionMap(), 1, ng_surf);
+        znd_surf.setVal(bogus_large_value);
+        znd_surf.ParallelCopy(*z_phys_nd[lev], 0, 0, 1, ng_surf, ng_surf, geom[lev].periodicity());
+        // Every node the stencil below reads (valid plus the x/y ghosts) must
+        // have been gathered; the slab has no z ghosts, so reduce over ng_surf
+        // rather than a scalar ghost count.
+        Real znd_max = ReduceMax(znd_surf, ng_surf,
+            [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4<Real const> const& a) -> Real
+            {
+                Real m = std::numeric_limits<Real>::lowest();
+                amrex::Loop(bx, [&] (int i, int j, int k) { m = amrex::max(m, a(i,j,k)); });
+                return m;
+            });
+        ParallelAllReduce::Max(znd_max, ParallelContext::CommunicatorSub());
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(znd_max < bogus_large_value,
+            "poisson_wall_dist: the surface nodes were not gathered onto every box");
+
         for (MFIter mfi(*walldist[lev]); mfi.isValid(); ++mfi) {
             const Box& bx = mfi.validbox();
             auto dist_arr = walldist[lev]->array(mfi);
             const auto zcc_arr = z_phys_cc[lev]->const_array(mfi);
-            const auto znd_arr = z_phys_nd[lev]->const_array(mfi);
+            const auto znd_arr = znd_surf.const_array(mfi);
             ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real z_surf = fourth * ( znd_arr(i,j,klo) + znd_arr(i+1,j,klo)
                                        + znd_arr(i,j+1,klo) + znd_arr(i+1,j+1,klo) );
@@ -300,9 +325,9 @@ void ERF::poisson_wall_dist (int lev)
     // definite form MLABecLaplacian documents (B = -1 with f = -h_zeta is
     // the same equation and gave the same iterates). Note: on a 3D fitted
     // mesh with dx != dz this multigrid diverges (residual 18x after the
-    // first cycle, 1e10 by iteration 100) with or without semi-coarsening;
-    // the anelastic projection on such a mesh fails the same way. Use
-    // erf.wall_dist_type = terrain_height there.
+    // first cycle, 1e10 by iteration 100) with or without semi-coarsening
+    // and independent of the box layout. Use erf.wall_dist_type =
+    // terrain_height there.
     // ****************************************************************************
     constexpr Real constA = zero;
     constexpr Real constB = one;

Tests/CTestList.cmake

@@ -1061,6 +1061,12 @@ endif()
 # Each case runs a short smoke deck and then its Python check script, which
 # compares planar-averaged numbers against stated targets with tolerances.
 # A clean exit alone is never the pass criterion.
+#
+# The decks run the anelastic projection with the FFT solver (erf.use_fft),
+# which no CI configuration builds. The flat decks are therefore run here
+# with the MLMG projection (erf.use_fft=false), and the terrain-fitted decks,
+# whose general-terrain projection has no non-FFT path, are registered only
+# when the build enables FFT (ERF_ENABLE_FFT).
 #=============================================================================
 find_package(Python3 COMPONENTS Interpreter QUIET)
 if(Python3_Interpreter_FOUND)
@@ -1111,8 +1117,9 @@ function(add_test_rans TEST_NAME CASE_DIR INPUT_FILE NSTEPS CHECK_SCRIPT)
     set(test_log "${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.log")
     set(check_log "${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.check.log")
     # The check script's exit code is the verdict; its table is echoed into
-    # the ctest output so a failure shows the measured numbers.
-    set(test_command sh -c "${MPI_COMMANDS} ${TEST_EXE} ${CURRENT_TEST_BINARY_DIR}/${INPUT_FILE} ${RUNTIME_OPTIONS} > ${test_log} 2>&1 && rm -f ${CURRENT_TEST_BINARY_DIR}/CHECK_FAILED && ( ${ERF_RANS_PYTHON} ${CURRENT_TEST_BINARY_DIR}/${CHECK_SCRIPT} --smoke ${CURRENT_TEST_BINARY_DIR}/${PLTFILE} > ${check_log} 2>&1 || touch ${CURRENT_TEST_BINARY_DIR}/CHECK_FAILED ) && cat ${check_log} && test ! -f ${CURRENT_TEST_BINARY_DIR}/CHECK_FAILED")
+    # the ctest output so a failure shows the measured numbers, and the tail
+    # of the run log is echoed when the executable itself exits non-zero.
+    set(test_command sh -c "${MPI_COMMANDS} ${TEST_EXE} ${CURRENT_TEST_BINARY_DIR}/${INPUT_FILE} ${RUNTIME_OPTIONS} > ${test_log} 2>&1 || ( tail -n 60 ${test_log} && false ) && rm -f ${CURRENT_TEST_BINARY_DIR}/CHECK_FAILED && ( ${ERF_RANS_PYTHON} ${CURRENT_TEST_BINARY_DIR}/${CHECK_SCRIPT} --smoke ${CURRENT_TEST_BINARY_DIR}/${PLTFILE} > ${check_log} 2>&1 || touch ${CURRENT_TEST_BINARY_DIR}/CHECK_FAILED ) && cat ${check_log} && test ! -f ${CURRENT_TEST_BINARY_DIR}/CHECK_FAILED")
 
     add_test(${TEST_NAME} ${test_command})
     set_tests_properties(${TEST_NAME}
@@ -1125,18 +1132,22 @@ function(add_test_rans TEST_NAME CASE_DIR INPUT_FILE NSTEPS CHECK_SCRIPT)
     )
 endfunction(add_test_rans)
 
-add_test_rans(RANS_Neutral_ABL_Flat     Neutral_ABL_Flat     inputs_neutral     40  check_neutral.py)
-add_test_rans(RANS_Neutral_ABL_Flat_Implicit Neutral_ABL_Flat  inputs_neutral     40  check_neutral.py RUNTIME_OPTIONS "erf.vert_implicit=true erf.fixed_dt=10")
-add_test_rans(RANS_Stable_ABL_Flat      Stable_ABL_Flat      inputs_stable      40  check_stable.py)
-add_test_rans(RANS_Convective_ABL_Flat  Convective_ABL_Flat  inputs_convective  40  check_convective.py)
-# terrain-fitted mesh: Poisson wall distance against the exact ridge distance,
-# and the same deck flattened (prob.hmax = 1e-6) against the analytic height
-add_test_rans(RANS_Neutral_Hill_2D        Neutral_Hill_2D      inputs_hill        40  check_hill.py)
-add_test_rans(RANS_Neutral_Hill_2D_Poisson Neutral_Hill_2D     inputs_hill        40  check_hill.py RUNTIME_OPTIONS "erf.wall_dist_type=poisson")
-add_test_rans(RANS_Flat_Fitted_2D         Neutral_Hill_2D      inputs_hill        40  check_flat_fitted.py RUNTIME_OPTIONS "prob.hmax=1e-6")
-add_test_rans(RANS_Flat_Fitted_2D_Poisson Neutral_Hill_2D      inputs_hill        40  check_flat_fitted.py RUNTIME_OPTIONS "prob.hmax=1e-6 erf.wall_dist_type=poisson")
-add_test_rans(RANS_Neutral_Hill_3D        Neutral_Hill_3D      inputs_hill3d      40  check_hill3d.py)
-add_test_rans(RANS_Neutral_Hill_3D_Poisson Neutral_Hill_3D     inputs_hill3d      40  check_hill3d.py RUNTIME_OPTIONS "erf.wall_dist_type=poisson")
+# flat meshes: MLMG projection so the tests run in every build
+add_test_rans(RANS_Neutral_ABL_Flat     Neutral_ABL_Flat     inputs_neutral     40  check_neutral.py    RUNTIME_OPTIONS "erf.use_fft=false")
+add_test_rans(RANS_Neutral_ABL_Flat_Implicit Neutral_ABL_Flat  inputs_neutral     40  check_neutral.py    RUNTIME_OPTIONS "erf.use_fft=false erf.vert_implicit=true erf.fixed_dt=10")
+add_test_rans(RANS_Stable_ABL_Flat      Stable_ABL_Flat      inputs_stable      40  check_stable.py     RUNTIME_OPTIONS "erf.use_fft=false")
+add_test_rans(RANS_Convective_ABL_Flat  Convective_ABL_Flat  inputs_convective  40  check_convective.py RUNTIME_OPTIONS "erf.use_fft=false")
+if(ERF_ENABLE_FFT)
+    # terrain-fitted mesh (FFT-preconditioned projection): wall distance against
+    # the exact ridge distance, and the same deck flattened (prob.hmax = 1e-6)
+    # against the analytic height, each with the terrain_height and Poisson paths
+    add_test_rans(RANS_Neutral_Hill_2D        Neutral_Hill_2D      inputs_hill        40  check_hill.py)
+    add_test_rans(RANS_Neutral_Hill_2D_Poisson Neutral_Hill_2D     inputs_hill        40  check_hill.py RUNTIME_OPTIONS "erf.wall_dist_type=poisson")
+    add_test_rans(RANS_Flat_Fitted_2D         Neutral_Hill_2D      inputs_hill        40  check_flat_fitted.py RUNTIME_OPTIONS "prob.hmax=1e-6")
+    add_test_rans(RANS_Flat_Fitted_2D_Poisson Neutral_Hill_2D      inputs_hill        40  check_flat_fitted.py RUNTIME_OPTIONS "prob.hmax=1e-6 erf.wall_dist_type=poisson")
+    add_test_rans(RANS_Neutral_Hill_3D        Neutral_Hill_3D      inputs_hill3d      40  check_hill3d.py)
+    add_test_rans(RANS_Neutral_Hill_3D_Poisson Neutral_Hill_3D     inputs_hill3d      40  check_hill3d.py RUNTIME_OPTIONS "erf.wall_dist_type=poisson")
+endif()
 
 # Largest stable time step of one closure under explicit anelastic, implicit
 # anelastic and implicit compressible integration (Timestep_Limits). The
@@ -1178,9 +1189,13 @@ function(add_test_rans_dt TEST_NAME CLOSURE)
     )
 endfunction(add_test_rans_dt)
 
-add_test_rans_dt(RANS_Timestep_Limits_kEqn      kEqn)
-add_test_rans_dt(RANS_Timestep_Limits_Deardorff Deardorff)
-add_test_rans_dt(RANS_Timestep_Limits_MRF       MRF)
+# The sweep's spin-up deck uses the FFT projection (the limits in its README
+# were measured with it), so the entries exist only in FFT builds.
+if(ERF_ENABLE_FFT)
+    add_test_rans_dt(RANS_Timestep_Limits_kEqn      kEqn)
+    add_test_rans_dt(RANS_Timestep_Limits_Deardorff Deardorff)
+    add_test_rans_dt(RANS_Timestep_Limits_MRF       MRF)
+endif()
 
 #=============================================================================
 # MOST reference height on flat stretched meshes