Guide 07 · September 10, 2026

Design the mesh and clock together

This dated design record describes the uniform-grid PhiFlow investigation. Fixed full-domain refinement is now implemented in AMReX / incflo; see the implementation and sensitivity checks and completed Oliver run. Completion is not a new convergence certificate.

The prescribed finite hierarchy gives us a candidate spatial metric and a predetermined logarithmic time grid. Neither is an accuracy certificate: the design must resolve the forcing, retain the full domain, and pass independent refinement checks.

The temporal launch gate failed

The full-interval 32³ pilot integrated from exact rest to t = 0.985 with source 14ea695. It halved both the maximum timestep and forcing-phase ceiling twice, holding the force-difference half-window at 10−6. Successive final-velocity differences were 1.22316 × 10−4 and 6.47590 × 10−5 in relative L². Their ratio was 1.889, corresponding to observed order 0.917, rather than the factor of approximately four expected in a second-order regime.

The launcher stopped without starting the 128³/192³ production pair. The pilot summary preserves the measurements; elapsed times for cached levels are reload times, not original integration costs. Small differences alone do not justify a Richardson error estimate when the expected order is absent. The current criterion remains a ratio between 2.5 and 6 and a finest-pair difference below 10−3; even passing this coarse-grid criterion would not establish spatial accuracy.

Remove an unintended source of temporal roughness

The old profile lookup used piecewise-linear interpolation. Moving similarity coordinates cross table knots, where the interpolated profile has slope jumps. The manufactured force differentiates this time-dependent profile. That gives a plausible numerical source of the reduced convergence order, independently of the intended pulse oscillations.

Source 9d2afb8 adds versioned cubic interpolation using SciPy CubicSpline, which is C² between table endpoints. Inputs are clamped to the table range to preserve the old exterior extension; global smoothness across clamped endpoints is not implied. This changes the numerical interpolant, not the analytical construction. Linear interpolation remains available for controlled comparisons, and old checkpoints explicitly retain it.

Four from-rest pilots were dispatched with immutable source snapshots: 32³ cubic interpolation; 64³ cubic interpolation; 32³ cubic interpolation with a five-times-smaller force-difference window; and a 32³ linear-table, background-only ablation. Each has three timestep levels. The no-pulse ablation also disables the pulse-phase ceiling, so it is not an otherwise identical-clock comparison. These are diagnostic jobs, not new production results.

The linear background-only ablation completed with successive differences 5.18569 × 10−4 and 3.05940 × 10−4, ratio 1.695. Reduced order therefore persists without the added pulses. A paired cubic background-only run was dispatched next to isolate interpolation under the same no-pulse clock. The remaining cubic comparisons were still in progress when this record was written.

A predetermined time grid

Let τ = 1 − t and s = −log τ. The implemented phase-rate bound, including quadratic combinations in the force, is C ≈ 42.693 for the default hierarchy. Set a maximum phase advance φ = 0.15:

Δs = φ/C ≈ 0.0035134;   tₙ = 1 − exp(−sₙ);   Δtₙ = τₙ [1 − exp(−Δs)].

This schedule is known before integration. At t = 0.985 it gives Δt ≤ 5.26088 × 10−5. The solver must still subdivide for velocity-dependent CFL, viscous stability, the early activation ramp, and exact output events. For the present uniform-grid explicit method, the additional diffusion cap is 0.12 Δx²/ν; a future nonuniform discretization needs its own stability bound.

For a single constant-frequency source, midpoint quadrature at phase advance 0.15 has relative integral-amplitude error φ/[2 sin(φ/2)] − 1 ≈ 0.094%. This is a useful scale estimate, not the global Navier–Stokes error. The separate display clock uses phase spacing 0.3; render cadence does not determine integration accuracy.

Spatial accuracy is stricter than avoiding aliasing

At the midplane, the finest retained radial phase is α log X, with α = 9.2, X = r²/(2 R²τ), and R = 0.52. Thus its local radial wavenumber is kᵣ = 2α/r, and wavelength λᵣ = πr/α. At the annulus midpoint X = 1.025 and t = 0.985, this gives λᵣ ≈ 0.03114.

For a resolved Fourier wave sampled with p cells per wavelength, the centered first derivative has amplitude ratio sin(2π/p)/(2π/p). The three-point Laplacian has ratio [sin(π/p)/(π/p)]². Eight cells leave about 10% first-derivative amplitude error; sixteen leave about 2.6%. These local stencil estimates are not solution-error bounds.

Full-box gridMidpoint cells/wavelengthFirst-derivative amplitude error
128³1.99Below Nyquist; aliased
192³2.9959.0%
512³7.9710.0%
1024³15.942.6%

Eight cells across this velocity phase require an even uniform resolution of at least 514³ at the midpoint, or 600³ at the inner half-height of the radial bump. Quadratic advection can double phase wavenumbers: the corresponding eight-cell requirements become 1028³ and 1200³. These are design estimates for selected locations, not claims that all such modes have substantial amplitude at the endpoint. Gates change their activity; envelopes, axial structure, correction terms, and the background need separate analysis.

The measured 192³ capacity preflight reached a 10.31 GiB macOS peak footprint. A crude cubic extrapolation gives about 195 GiB at 512³, before revisiting storage and solver design. Four machines running this solver do not pool their RAM. We therefore need a more efficient discretization, not just a larger uniform run.

A fixed pre-refined mesh, not a smaller box

September 12 update: the full-domain AMReX projection pilot now passes manufactured-solution checks on four fixed Cartesian levels. It reaches 1024³-equivalent core spacing with measured operator memory and spatial convergence. The forcing-informed mesh and full momentum evolution remain unvalidated; the design constraints below still apply.

A geometric radial mesh follows directly from the phase: requiring at least p cells per cycle gives Δlog r ≤ π/(αp). Accounting for doubled radial phase halves that limit. With p = 8, the doubled-phase rule permits approximately 2.16% radial growth per cell in the relevant annulus.

The mesh must cover the union of active supports over the entire start-from-rest trajectory, not just the final core. It must also resolve axial and angular phase gradients, the axis, amplitude envelopes, the background field, and the original exterior taper. Those requirements define a spatial metric to evaluate before selecting cell counts. A finite near-axis patch and compatible pressure solve are necessary; geometric spacing cannot continue to r = 0.

This is a design proposal, not an implemented mesh. The current solver uses one constant Δx and a uniform-grid FFT projection. A stretched grid needs a compatible nonuniform discretization and projection; changing coordinate arrays alone would be incorrect. The existing core-refined preset shrinks the box and changes the taper, so it does not implement this proposal.

Reproduce and test the design

uv run python -m navier_stokes_sim.numerical_budget \
  --t-end 0.985 --points-per-wavelength 8 \
  --output artifacts/numerical-budget.json
uv run python -m scripts.check_forcing_timestep \
  --resolution 32 --profile-interpolation cubic \
  --derivative-epsilon 0.000001 --output outputs/cubic-pilot

The generated design record contains the configuration, wavelength estimates, stencil responses, candidate radial spacings, and complete phase-only time grid. Regression tests compare the symbol formula with the actual centered stencil and the precomputed clock with the runtime phase ceiling.

Next: establish temporal order at fixed spatial resolution and force-difference window; check that window independently; then refine space with a sufficiently tight time clock. Compare velocity, vorticity, force components, spectra, and core/taper quantities at matched times. Manufactured-target tracking alone can hide spatial underresolution because the discrete force is built from the discrete target. No finite grid can demonstrate an infinite singular value.