Compute only ∂ξ∂x[3,3] for the vertical WVector→Contravariant3 projection

[imreddyTeja]

Summary

Speeds up the vertical physical→contravariant projection used by the vertical finite-difference operators (e.g. divᵥ applied to a WVector).

The generic path

_to_components_type(::Contravariant, ::OrthonormalTensor, lg) = lg.∂ξ∂x * v

materializes the full 3×3 inverse ∂ξ∂x = inv(∂x∂ξ), even though for a purely-vertical vector v = (0, 0, w) only the (3,3) entry survives:

(∂ξ∂x · v)₃ = ∂ξ∂x[3,3] · w

This PR adds a specialized project(Components{Contravariant,(3,)}, ::WVector, lg) that computes that single entry as minor(3,3)/det(∂x∂ξ), avoiding the other eight cofactors of the inverse.

Both the CPU stencil (contravariant3/Jcontravariant3) and the GPU operator-matrix path (project_row2_for_mul → project(Contravariant3Axis, WVector, lg)) funnel through this method, so both benefit.

Why it's exact (no geometry guard)

∂ξ∂x[3,3] = minor(3,3)/det is the exact (3,3) entry of inv(∂x∂ξ) for any metric — including terrain-following, where the horizontal–vertical coupling enters through the determinant. The tempting shortcut 1/∂x∂ξ[3,3] is only valid for an orthogonal vertical (it gives e.g. 0.333 vs the true 0.362 on a coupled metric) and is deliberately not used, so no geometry-type guard is required.

Validation

• New unit test test/Geometry/contravariant3_projection.jl (registered in runtests.jl), 8009 assertions passing: machine-precision agreement with inv(∂x∂ξ)[3,3] across 2000 random well-conditioned and terrain-following metrics (Float32 + Float64), type stability, correct dispatch, and that non-vertical / already-contravariant sources are unchanged.
• CPU microbenchmark: 6.0 vs 8.3 ns/call (1.39×), allocation-free.

Remaining before merge

• GPU end-to-end benchmark of DivergenceC2F / DivergenceF2C on a 63-level VIJFH sphere (Nq=4). Expectation: the divergence kernels move from the projection-bound ~40-register regime toward the no-projection ~32-register regime (a probe showed WVector-projection divergence at 40 regs / 233 µs vs a Contravariant3 operand at 32 regs / 125 µs). (In progress — waiting on a GPU allocation on the cluster.)
• Full Operators test suite.

cc @imreddyTeja

---

Note: independent of the tr/cart-inds eager-FD work; this change benefits main directly (main pays the same full-inverse projection cost).

🤖 Generated with Claude Code

https://claude.ai/code/session_01DHA3sRpzBXgPDEJuQrNoBZ

[imreddyTeja]

GPU benchmark result (A/B) — the divergence win does not materialize on main

Ran the A/B on an A100 (z=63, helem=30, Nq=4, VIJFH, Float32), same base commit, only src/Geometry/conversions.jl differing:

op	without this change	with this change
DivergenceF2C(WVector)	30 regs / 201.6 µs	32 regs / 199.8 µs
DivergenceC2F(WVector, SetDivergence)	32 regs / 213.1 µs	32 regs / 207.4 µs
divᵥ(gradᵥ(·))	26 regs / 146.7 µs	26 regs / 146.7 µs
divᵥ(interpᵥ(·))	22 regs / 122.5 µs	22 regs / 122.6 µs
divᵥ(f·gradᵥ(uₕ))	30 regs / 241.0 µs	30 regs / 240.9 µs

Register counts are unchanged (±2) and times are within noise (~1–3%). On main, the compiler already avoids materializing the full inv(∂x∂ξ) in the vertical-divergence kernel, so eliminating it explicitly buys nothing on the GPU here.

The projection-bound behaviour I originally measured (vertical divergence at ~40 registers, ~2× slower than the no-projection case) was on a separate eager-FD refactor branch, not main; that register pressure does not reproduce on main.

Net: this change is still exact for all geometries and ~1.39× faster on the CPU in isolation (allocation-free), but it does not speed up main's GPU divergence kernels. Given it adds a specialized method for no measurable GPU benefit on main, it may be best to hold this until the eager-FD path where it actually helps is on main, or treat it purely as a small CPU/clarity improvement. @imreddyTeja — your call on whether to keep it open.