Merge the k-equation RANS branches (claude-RANS, claude-RANS-implicit) into ERF-Hazard - #405
Merged
Conversation
…CTest hook Adds Exec/CanonicalTests/Canonical_RANS with the eight-phase plan for a minimal working one-equation k RANS (Axell & Liungman 2001), a neutral Ekman-layer deck on flat ground under MOST, and a check script that reads the plotfile with a standard-library-only AMReX reader and compares planar-averaged numbers against stated targets with tolerances. add_test_rans in Tests/CTestList.cmake runs a 40-step smoke deck and then the script; the script's exit code is the verdict. RANS_Neutral_ABL_Flat passes on 1 and 2 ranks, which agree to 2e-15. The 12 h run records the "before" state in RESULTS.md: the log law holds within 7 %, but the wall cell keeps only 80 % of its Dirichlet k through a step and KE(0)/u*^2 is 2.60 against the AL01 value 3.23. Phase 2 fixes the wall condition. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…e step
With erf.dirichlet_k the surface layer writes the AL01 Eq. 16 value of k
into the first cell of S_old once per step. The logical BC for RhoKE was
then switched to ext_dir, so the diffusion applied a one-sided stencil
against a zero ghost cell and drained a fifth of the wall k before the
step ended; the converged wall k sat 20 % below the AL01 value and the
implicit KE solve was silently switched off by the BC sanity check.
Keep foextrap for RhoKE at a surface-layer wall (ghost cell equals the
first cell, zero flux through the wall face), restore the first-cell
value from S_old after every RK stage in erf_slow_rhs_post, and pin the
bottom row of ImplicitDiffForStateLU_{N,S,T} for RhoKE so the implicit
solve treats it as a Dirichlet neighbour.
Neutral_ABL_Flat, 12 h: KE(0)/u*^2 = 3.232 against 3.2325 (was 2.60),
k_start/k_end = 1.000 (was 1.244), interior profile unchanged, restart
bit-exact, and the wall-cell check now runs in the CTest smoke mode. The
implicit path was exercised with a compressible variant since anelastic
disables it. erf.dirichlet_k keeps its default of false.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Without the Dirichlet wall value the first cell cannot resolve the near-wall shear production and k there settles at about half the Axell & Liungman equilibrium u*^2/Cmu0^2 (12 h neutral run: 1.59 against 3.23), while the mean wind still follows the log law because MOST supplies the stress. Warn at startup in that configuration; the default stays false. Every deck that enables kEqn now sets the flag (Askervein added; the canonical deck already had it). Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…table length, unit tests Factor the Axell & Liungman closure relations into ERF_RANSClosure.H (namespace AL01) so ComputeTurbulentViscosityRANS and the unit tests evaluate the same expressions; the refactor is bit-identical to phase 2. Validate Cmu0, Cb, max_geom_lscale, Rt_crit <= 0, Rt_min < Rt_crit and Rt_min > -3.6 (the poles of Eqs. 31-32) at input time. Add erf.tke_floor (runtime floor on k; distinct from erf.tke_min, the initial value) and the opt-in erf.implicit_tke_dissipation, which drops the explicit sink from the source and divides the update by 1 + dt c with c = diss_old / (rho k)_old, half-weighted on the anelastic stage 1. Two robustness changes came out of the new gtests and the review against the paper and the Kynema KLAxell implementation: the Burchard & Petersen smoothing is rewritten as Rt_crit + a x / (x + a), the same function without the cancellation that returned -2 for Rt = -1e16 and +1.4e14 for Rt = -1e30; and the unstable length evaluates Eq. 28 once with the smoothed Rt from the geometric length, bounded by about 1.31 l_g, instead of the two-pass corrector that iterated the fixed-point map of Eq. 26, which has no fixed point in strong convection (AL01 p. 78). Also: fill the derived diffusivity components on the tilebox only (the wrapper fills ghost cells), drop unused parameters, correct the surface heat-flux index in the TKE-source comment, and use the level's own periodicity in the wall-distance Poisson solve. Neutral_ABL_Flat: all checks unchanged to the digits shown; the 4x-dt run showed the anelastic integrator's explicit vertical diffusion is the real time-step limit (dz^2/2K about 19 s), not the dissipation; a compressible dt = 60 s pair with implicit vertical diffusion passes with the dissipation explicit and implicit. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…er anelastic Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…ap, stable and convective decks Two opt-in closure options: erf.rans_consistent_diffusivities makes the horizontal heat and every scalar and moisture diffusivity follow the AL01 scalar stability function (rho cmu' sqrt(k) L) like the vertical heat diffusivity, instead of Pr_t and Sc_t times the eddy viscosity; and erf.rans_lscale_from_pblh caps the geometric length at kappa 0.1 zi from the surface layer's MYNN25 PBL-height diagnostic, clamped to [erf.rans_lscale_min, erf.max_geom_lscale]. The neutral deck is bit-identical with both off. Two new Canonical_RANS cases with numeric check scripts and CTest smoke entries: Stable_ABL_Flat (GABLS1, 9 h) passes every check first time, u* 0.244, jet 1.23 Ug at 154 m, depth 134 m; Convective_ABL_Flat (MS94-B sounding, 0.24 K m/s, 4 h) closes the column heat budget to 0.02 %, keeps the inversion at 1020 m and carries the AL01 buoyancy term at the wall. It runs at dt = 2 s: at 5 s the eddy viscosity reaches the explicit vertical-diffusion limit of the anelastic integrator (40 m2/s) and the run aborts, the same limit phases 9 and 10 address. The mixed-layer spread tolerance (2 K) reflects the superadiabatic lapse a local-K closure keeps without countergradient transport; Kynema's KLAxell has none either and its Prandtl function matches cmu'/cmu. Shared check code moves to rans_checks.py; the dissipation-lag check allows 10 % in smoke mode (early transient) and 5 % in physics mode. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
… deck The Poisson wall distance (Tucker 2003) took the cell's low-face fluxes as its gradient, half a cell off centre in every direction, which overstates |grad phi| by dz/2 and shortens every distance by z dz/(2H) on a flat mesh, 0.78 % for 64 cells. Form a cell-centred gradient instead from centred differences with the cell-centre terrain metrics (chain rule for a mesh deformed in z); it needs one ghost cell of phi and the cell's own nodes. On a flat fitted mesh the distance is now exact to 1e-6 m above the first cell (1.5 cm there, the odd-reflection Dirichlet ghost); on the ridge the mean error drops from 1.5 % to 1.0 %. New Canonical_RANS case Neutral_Hill_2D: periodic Witch of Agnesi ridge (h 100 m, L 500 m) on a terrain-fitted mesh under MOST, one cell wide in y, anelastic with the MLMG projection. Its check script compares the wall distance with the exact distance to the curve, and after 6 h measures the crest speed-up (0.55, 0.45, 0.39 in the lowest cells against the Jackson-Hunt estimate 0.4) and the upstream log law (within 4 %). A flat-fitted variant of the same deck (prob.hmax = 1e-6) is the second CTest entry; both pass. erf_plotfile.py gains a full-field reader for the terrain checks. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…estart on terrain Add erf.wall_dist_type = terrain_height: the height above the local surface projected on its normal, no linear solve (after the amr-wind immersed terrain and Kynema). It is exact to 1e-10 on a flat fitted mesh and closer to the true distance than the Poisson solve on both hills (mean 0.02 % and 0.01 % against 1.0 % and 0.3 %); the hill decks use it and three _Poisson CTest variants keep the Poisson path covered. The Poisson solve is now posed in the positive-definite form MLABecLaplacian documents (same iterates as before). New case Neutral_Hill_3D: a periodic radial Witch of Agnesi hill at dx = dy = dz = 40 m; after 4 h the crest speed-up is 0.41, 0.28, 0.22 in the lowest cells against the axisymmetric estimate 1.6 h/L = 0.32, and the upstream log law holds within 10 %. Restart is bit-exact on the 2D and 3D terrain decks; the Askervein deck runs clean. Recorded in the plan, outside its scope: on a 3D fitted mesh with dz different from dx the initial field's divergence is a deterministic 1.788e139 before the first projection and the run aborts, independent of the closure, the boundaries, the terrain source, stretching, the box layout and the projection solver, while Askervein runs at the same aspect ratios; a separate uninitialised read in the w boundary fill trips the invalid-operation trap in the first advance. Reproducer in PLAN.md. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…agnostics Add Docs/sphinx_doc/theory/RANS.rst, describing the Axell & Liungman one-equation closure as implemented: the eddy viscosity and diffusivity, the k equation and its dissipation, the geometric and stratified length scales with the bounded unstable form, the stability functions with the cancellation-free smoothing, the wall condition on k, the two wall distances on a fitted mesh, the limitations and an input table. Link it in the THEORY toctree, document the Canonical RANS regression cases in RegressionTests.rst, and add a top-level README for the case directory. Three new EddyDiff components, written by the closure, expose the smoothed turbulent Richardson number and the two stability functions as the plot variables Rt, cmu and cmu_prime; the neutral smoke check verifies them against AL01 Eqs. 31-32 to round-off and Kmv against rho cmu sqrt(k) L. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…istency check Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…stic integrator Opt-in: erf.vert_implicit = true (or an explicit erf.vert_implicit_fac) in the inputs turns the column solve on for the scalars under anelastic; without it the run stays fully explicit as before, and a banner says how to enable it. Momentum stays explicit (its solve lives in the substepping path). The anelastic second stage is trapezoidal and rebuilds the first-stage tendency from the state difference, which already carries the first stage's implicit increment, so the pre-stage theta solve and the post-stage KE/moisture solves now act with half the step on that stage. With the full step the convective column gained 1.45 times the surface heat flux per step (stage decomposition with vert_implicit_fac 1 0 0 / 0 1 0 / 0 0.5 0); it is 0.975 for every combination now. The column solves take each box's z extent as the column, so level 0 is no longer decomposed in z when the solve is on and init aborts if any box is cut in z (this also closes a silent error in the compressible path). check_for_negative_theta now trips on NaN, which it let through. The convective RANS deck runs at dt = 5 s with the solve on (it aborted at 5 s explicit); a CTest entry runs the neutral deck at dt = 10 s implicit. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…explicit face flux The k-equation closures (Deardorff, Axell-Liungman) write -K_h dtheta/dz from the cell-centred gradient into Hfx3 at the start of the step and the TKE source reads it there (AL01 Eq. 15; Kynema forms it the same way). Since "Optimize Post" (erf-model#2310) the explicit diffusion operator overwrote that value above the first cell with the face flux times the explicit fraction, so with vert_implicit_fac = 1 the buoyancy production and destruction of k vanished everywhere but the first cell (the convective RANS case carried a third of its KE and the mean theta was 1.3 K off), and even fully explicit the term sat half a cell low. The overwrite is removed; the surface layer still writes Hfx3(i,j,0). Verified by restarting the neutral RANS deck at 6 h and taking one step implicit and explicit: KE differed by 1e-6, growing linearly in time and independent of dt, only when the buoyancy term was present; after the fix the two agree to 8e-9 in one step and to 5e-7 after 10 s, and over the full runs to 3e-4 m/s (neutral 12 h), 4e-5 m/s (stable 9 h) and 3e-3 m/s (convective 4 h, dt 5 s implicit against dt 2 s explicit). The Deardorff and implicit-diffusion gold tests pass. Theory page, plan and results updated for phase 9. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Phase 9 left momentum explicit under anelastic on the grounds that its solve lives in the substepping path. That was wrong: with implicit_before_substep, which is forced true when there is no substepping, ERF_ImplicitPre.H folds the momentum solve into the slow tendency before the stage update, so under anelastic the momenta are diffused and then projected, which is the order the divergence constraint needs. Removing the flag is the whole change; the half-step on the trapezoidal second stage from phase 9 applies to the momentum solve unchanged. erf.vert_implicit = true under anelastic now covers u, v, theta, k and moisture. w stays explicit unless built with ERF_IMPLICIT_W, as in the compressible path, and the solve remains opt-in under anelastic. Neutral deck, 12 h: dt = 5, 10, 20, 30 and 60 s all pass every physics check, and at 60 s (6 times the explicit limit) the profiles are within 2.9e-3 m/s of the explicit dt = 5 s run. Compressible at dt = 5, 20 and 60 s lands within 1.8e-3 m/s of the same reference. The convective deck reaches dt = 20 s, ten times its original explicit step, and the stable deck dt = 8 s. On the 2D hill the advective Courant number binds first (it fails at dt = 3 s identically with the solve on and off), and there the terrain-fitted momentum solve reproduces the explicit answer to 1.2e-4 m/s. Restart at dt = 20 s is bit-exact. Recorded in PLAN.md and RESULTS.md, with one caveat found while verifying: the implicit theta solve makes the answer sensitive to the box decomposition at the 1e-5 level in wind, because the column tridiagonal spreads a one-ulp difference over the whole column in a single step and the near-neutral buoyancy term amplifies it. The k and momentum solves are decomposition invariant to round-off on their own, and with Smagorinsky instead of the k equation the theta solve stays at 1e-12. No physics check moves. plot_dt_overlay.py draws the comparison. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…y and Deardorff too The column solve is in the dycore, not the closure. On the neutral geometry under anelastic, explicit Smagorinsky fails at dt = 20 s and explicit Deardorff already at dt = 5 s, while both run at 5, 20, 60 and 120 s with erf.vert_implicit = true. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ressible Neutral deck, 12 h: anelastic passes every check up to dt = 240 s and fails at 480 s; compressible passes to dt = 120 s and loses seven checks at 240 s. At their best passing steps the anelastic run costs 1 s of wall time against 433 s, since it pays one FFT solve per step rather than thousands of acoustic substeps. Before phases 9 and 10 the ordering was inverted, the compressible path having had the column solve all along. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The phase 10 note blamed the k equation's buoyancy term. That is wrong. Chasing it properly: it reproduces on a single rank with four boxes, so MPI reduction order is not involved; it survives local MOST averaging, a nulled buoyancy term, dirichlet_k off, a TKE floor, and Deardorff in place of the k equation; two identical runs agree bitwise; and fab.init_snan with the invalid trap armed finds no uninitialised read. What holds up: the deck stays exactly uniform horizontally, so the two decompositions differ only through the round-off of the projection, which leaves a w of order 1e-19 where w should be zero. The state first differs at step 3 by one ulp of rho theta and grows before saturating. The growth needs both the implicit solve and the anelastic projection (the compressible path with the same solve stays at 5e-12) and scales with the implicit increment: 9e-7 in wind at dt 5 s, 1.3e-5 at 20 s and 3.6e-4 at 60 s, while theta alone at vert_implicit_fac = 1 0 0 is clean at one ulp. A scalar amplification analysis of the fold-in gives a factor in (0,1] for every mode, so it is not the theta update alone. Not fixed. The mechanisms ruled out are recorded so the search does not have to be repeated. It moves no physics check; at dt 60 s the spread is 3e-5 relative against an implicit-explicit difference of 2.5e-4. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
New case Timestep_Limits: a neutral column on a RANS-like grid (4 x 4 x 200,
dx = 800 m, dz = 5 m) that measures the largest stable step of kEqn,
Deardorff and MRF under explicit anelastic, implicit anelastic and implicit
compressible integration. sweep_dt.py spins each closure up for 1 h at
dt = 5 s with the implicit solve (anelastic and compressible checkpoints),
restarts over a ladder from 0.125 s to 1024 s for 200 steps per rung, and
reports the largest rung below the first failure.
Largest passing step (first failure):
closure explicit anel. implicit anel. implicit comp. dz^2/(2K)
kEqn 2 (4) 256 (512) 64 (128) 2.13
Deardorff 0.25 (0.5) 512 (1024) 64 (128) 0.34
MRF 0.5 (1) 256 (512) 8 (16) 0.58
The checks require the explicit step within a factor 2 of dz^2/(2K) (0.94,
0.74, 0.86) and both implicit integrators at 8 times it or more. With the
solve switched off in the implicit anelastic mode the MRF sweep fails that
check. MRF under implicit compressible stops at 8 s; the 16 s run keeps
going with |u| in the thousands of m/s, with ERF's own substep count as
with pinned substeps. Not investigated.
Two things the driver has to do: lift erf.change_max, which caps step growth
at 10 % per step even with erf.fixed_dt, so a restart would otherwise ramp
up from 5 s; and pin erf.fixed_mri_dt_ratio, because the substep count sized
from a state that is going unstable reaches billions and the run hangs.
Restarts from the spin-up checkpoints are bit-exact for all six.
CTest entries RANS_Timestep_Limits_{kEqn,Deardorff,MRF}, 60-90 s each on one
rank in Release, labelled rans and dt_sweep but not regression, which the CI
runs in Debug.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…and Canonical_RANS cases Brings in the twelve commits of claude-RANS (open upstream as erf-model#3971) ahead of its upstream merge; both branches share the development base 489a2aa. Tests/CTestList.cmake conflicted only because both sides appended a section at the same place; the fire and dust smoke tests and the Canonical RANS cases are both kept. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…on under anelastic and the time-step sweep Brings in the seven commits of claude-RANS-implicit on top of claude-RANS: the implicit column solve of theta, KE, moisture and momentum under the anelastic integrator (opt-in with erf.vert_implicit = true), the TKE buoyancy term using the closure's cell-centred heat flux, NaN detection in check_for_negative_theta, and the Timestep_Limits dt sweep. It merged without conflicts. Its level-0 full-column abort in update_diffusive_arrays duplicates ERF-Hazard's check_stacked_boxes_in_z, which runs first, and it does not touch fine levels, so the whole-column joining of fine boxes from #404 is unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
This brings the k-equation RANS work into ERF-Hazard now, ahead of its upstream review. It adds two merge commits, in this order:
23285ff2cmergesclaude-RANS, 12 commits, open upstream as One-equation k RANS: wall condition, closure fixes, terrain wall distance, Canonical_RANS regression cases erf-model/ERF#3971. It contains the one-equation k RANS closure fixes, theExec/CanonicalTests/Canonical_RANScases with their Python check scripts, seven gtests and the theory pageDocs/sphinx_doc/theory/RANS.rst.25bffd236mergesclaude-RANS-implicit, 7 commits on top ofclaude-RANS. It contains:check_for_negative_theta;Canonical_RANS/Timestep_Limitsstep sweep.Both branches start from
489a2aa64, the development commit ERF-Hazard last synced, so nothing else from development comes in.Conflicts and overlaps
Tests/CTestList.cmake, where both sides appended a section at the same place. The fire and dust smoke tests and the RANS cases are both kept.check_stacked_boxes_in_z(Stop at start-up when boxes stacked in z meet the column solves or the surface layer #400) runs first and gives the more detailed message.ABL_ZSplit_*_aborttests and theBubble2D/Terrain2Levfine z-split parity tests pass.ERF_TurbStruct.H,ERF_SetupDiff.H,ERF_ComputeTurbulentViscosity.cpp,ERF_IndexDefines.H,ERF_PlotfileSelection.H,ERF_SlowRhsPost.cppandTests/Unit/CMakeLists.txt. ERF-Hazard's changes there are fire, dust and MRF additions, the PBL full-column check and thetheta_refFPE guard. They sit beside the RANS changes, with no shared lines.Behaviour changes on ERF-Hazard
-K_h dtheta/dz. The explicit diffusion operator used to overwrite it with the face flux: with the implicit solve that zeroed it above the first cell, and in explicit runs it sat half a cell low.ABL_MOSTandDeardorff_stationarystill pass.check_for_negative_thetaaborts on a NaN rho or rho thetaerf.vert_implicit = true, or an expliciterf.vert_implicit_fac, now turns the column solve on under anelastic. The default stays explicit.amr.max_grid_size_zbelow nz and x/y uncapped now run on whole columns instead of stopping at the #400 checkNew options from erf-model#3971, all off by default:
erf.tke_floor,erf.implicit_tke_dissipation,erf.rans_consistent_diffusivities,erf.rans_lscale_from_pblhwitherf.rans_lscale_min, anderf.wall_dist_type = poisson | terrain_height. There are also new plot variablesRt,cmuandcmu_prime.The largest stable step on the RANS-like
Timestep_Limitscolumn (dx = 800 m, dz = 5 m), with the first failing step in brackets:MRF under implicit compressible stopping at 8 s is not investigated.
Tests
The build is Release with MPI, fire, dust and FFT, all warnings on, and has no warnings outside
Submodules.ctest -L unit -LE parallel(536) andctest -L parallel(16) pass.ctest -L regressionpasses 122 of 122. That covers the 10 RANS cases, fire smoke and abort tests, the tiling, box and fine z-split parity tests, and the gold tests includingABL_MOST*,ABL_MOST_IMP_DIFF*andDeardorff_stationary.erf_exechad been built:amrex_fcompareand the checker executables were missing. After building every target, those 80 pass.ctest -L dt_sweeppasses 3 of 3. It is kept out ofregressionbecause the CI runs that label in Debug.git diff --check. The exceptions are blank lines at the end ofPLAN.mdanderf_plotfile.py, which the style job does not check.Later, upstream
erf-model/ERF squash-merges pull requests, so erf-model#3971 will land in development as one new commit, not these SHAs.
git rebase --onto upstream/development claude-RANS claude-RANS-implicit.🤖 Generated with Claude Code