diff --git a/Docs/sphinx_doc/RegressionTests.rst b/Docs/sphinx_doc/RegressionTests.rst index bf9382263b..ad888d43e8 100644 --- a/Docs/sphinx_doc/RegressionTests.rst +++ b/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, diff --git a/Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/README.md b/Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/README.md index f671ba705e..5933dcaa2c 100644 --- a/Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/README.md +++ b/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 diff --git a/Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/inputs_hill3d b/Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/inputs_hill3d index 6e1d3ac07a..0f1f3eaf79 100644 --- a/Exec/CanonicalTests/Canonical_RANS/Neutral_Hill_3D/inputs_hill3d +++ b/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 diff --git a/Exec/CanonicalTests/Canonical_RANS/PLAN.md b/Exec/CanonicalTests/Canonical_RANS/PLAN.md index 40411c6365..468d1eaa02 100644 --- a/Exec/CanonicalTests/Canonical_RANS/PLAN.md +++ b/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 diff --git a/Exec/CanonicalTests/Canonical_RANS/README.md b/Exec/CanonicalTests/Canonical_RANS/README.md index 6a73fb5ebc..cec8ae1b3b 100644 --- a/Exec/CanonicalTests/Canonical_RANS/README.md +++ b/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_.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`). diff --git a/Exec/CanonicalTests/Canonical_RANS/RESULTS.md b/Exec/CanonicalTests/Canonical_RANS/RESULTS.md index b0280751e9..137d974323 100644 --- a/Exec/CanonicalTests/Canonical_RANS/RESULTS.md +++ b/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 diff --git a/Source/LinearSolvers/ERF_PoissonWallDist.cpp b/Source/LinearSolvers/ERF_PoissonWallDist.cpp index eb9093d96b..feff717fb1 100644 --- a/Source/LinearSolvers/ERF_PoissonWallDist.cpp +++ b/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 const& a) -> Real + { + Real m = std::numeric_limits::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; diff --git a/Tests/CTestList.cmake b/Tests/CTestList.cmake index 830d6c8c04..568be4bf44 100644 --- a/Tests/CTestList.cmake +++ b/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