dt has an accuracy limit BELOW its stability limit; between them near-surface resolved w is silently replaced by grid-scale noise

[tyatharva]

Summary

dt is a mandatory user constant that FastEddy never checks (SRC/TIME_INTEGRATION/time_integration.c; the only bound is FLT_MIN..FLT_MAX). Bisecting it on a production grid, we found two distinct thresholds, and the lower one is silent:

• above CFL_3d ≈ 1.79 the run goes non-finite and prints ****CORRUPTED***
• above CFL_3d ≈ 1.64 — but below 1.79 — the run completes, exits 0, prints no warning, and the resolved vertical velocity in the lowest ~3 levels is grid-scale acoustic noise instead of turbulence

Everything else looks fine in that window: u, v, theta, fricVel, and the whole profile above ~45 m are indistinguishable from a correct run. That is what makes it dangerous — for a footprint/dispersion application the near-surface w statistics are exactly the quantity of interest, and they are wrong by four orders of magnitude with no indication anything happened.

Here CFL_3d = dt * c * sqrt(1/dx² + 1/dy² + 1/dz²), with c ≈ 347.2 m/s at 300 K. FastEddy is fully compressible with RK3 (Wicker–Skamarock 2002) and no acoustic sub-stepping, so the acoustic CFL is the relevant number.

Reproduction

FastEddy v5.0.1 (e0cd2f3), single RTX 4080, sm_89, 1 rank, CUDA 11.8.
Case: dry neutral, flat, uniform z0 = 0.03 m, doubly periodic, U_g = 10 m/s, hydroBCs = 2,
advectionSelector = 3, turbulenceSelector = 1, TKESelector = 1,
filterSelector = 1, filter_6thdiff_vert = 1, filter_6thdiff_vert_coeff = 0.03
(i.e. the Example01_NBL numerics).

Diagnostic is the horizontal-plane variance of resolved w at the first interior level divided by the second (k0/k1). Physically this must be < 1 — the surface is impermeable, so w variance has to grow away from the wall.

Grid A — 434 x 146 x 122, dx = dy = 10 m, dz_surface = 10.0 m

Only dt was varied.

dt (s)	CFL_3d	k0/k1	<w'w'> at level 0 (m²/s²)	outcome
0.0250	1.503	0.16	3e-5	clean (production setting)
0.0272	1.636	0.27	6e-5	clean
0.0273	1.642	0.36	8e-5	degrading
0.0275	1.654	8.94	1.97	silent corruption — exit 0, no CORRUPTED
0.0300	1.794	—	—	NaN, ****CORRUPTED***

<w'w'> at the first level changes by ~33,000x across a 0.7% change in dt. The transition is sharp, not gradual. The lag-1 horizontal autocorrelation of w at 5 m collapses from +0.95 (resolved eddies) to +0.08 (grid noise) across the same step, while it stays at +0.95 at 85 m.

We eliminated the obvious alternatives by single-variable tests, all of which left k0/k1 = 8.94 unchanged: vertical stretching on/off, thread-block shape (4x4x16 vs 1x4x64), roughness length, spanwise domain width, damping-layer depth, and latitude. The only two changes that fixed it — coarsening dz_surface and coarsening dx,dy — both act by lowering CFL_3d, which is what pointed at dt.

Grid B — 146 x 50 x 90, dx = dy = 30 m, dz_surface = 20.0 m

A 3x coarser horizontal and 2x coarser vertical spacing, restarted from a 6-hour spun-up state of the same case and run 300 s at each dt:

dt (s)	CFL_3d	k0/k1	<w'w'> at level 0 (m²/s²)	outcome
0.0625	1.491	0.16	7e-4	clean (production setting)
0.0670	1.599	0.17	9e-4	clean
0.0713	1.701	7.44	577	silent corruption — exit 0, no CORRUPTED

The threshold sits between CFL_3d 1.60 and 1.70 on both grids, despite dx changing by 3x and dz by 2x. That is what identifies it as a property of the 3-D acoustic CFL rather than of a particular resolution — and it is why a user cannot avoid it by picking a "safe-looking" grid.

Why the tutorials don't show it

The shipped cases sit just below the threshold:

case	CFL_3d
Example01_NBL	1.603
Example03_SBL	0.904

Example01_NBL is within 2.5% of where we measure the onset. Its dt is therefore correct but has essentially no margin, and any user who scales a tutorial's dt to a new grid by 1-D reasoning (dt*c/min(d), which is the more intuitive form and gives a limit ~1.9x looser) will land above the accuracy threshold while remaining below the stability threshold.

Suggestion

Nothing here requires changing the numerics. Two cheap, non-breaking additions would close the trap:

1. At startup, compute and print dt * c * sqrt(1/dx² + 1/dy² + 1/dz_min²) alongside the echoed parameters, using the base-state sound speed. Right now a user has no way to see this number without deriving it.
2. Warn (not abort) above a conservative value. 1.6 would have flagged every corrupted case we produced; it would also flag Example01_NBL at 1.603, so ~1.65 may be the better line.

Happy to open a PR for the startup diagnostic if that would be useful. We are running v5.0.1 unmodified and would rather not carry a fork patch for this.

Context

This came out of setting up FastEddy as the target generator for a flux-footprint emulator (LES + backward Lagrangian dispersion), where the near-surface w statistics are the signal. We hit it as a physically implausible footprint before we understood it as a dt problem, hence the level of detail above.

Edit: Claude opened this issue on its own, so if something seems incoherent or wrong I apologize for that