Honour phi_scale end to end, guard truncating rotations, drop the dead FFT cache

.github/workflows/ci.yml

@@ -204,6 +204,7 @@ jobs:
           mpiexec -n 4 julia --project=. test/parallel/test_mpi_plan_preflight.jl
           mpiexec -n 2 julia --project=. test/parallel/test_mpi_transpose_operand_preflight.jl
           mpiexec -n 2 julia --project=. test/parallel/test_mpi_comm_cleanup.jl
+          mpiexec -n 4 julia --project=. test/parallel/test_mpi_2d_alignment.jl
           mpiexec -n 2 julia --project=. test/parallel/test_mpi_ad_tangent_spaces.jl
           mpiexec -n 8 julia --project=. test/parallel/test_mpi_audit_fixes.jl
         env:

CHANGELOG.md

@@ -37,8 +37,91 @@ coefficient was `1/2π` too small — they inverted neither `synthesis_axisym` n
 the m=0 column of the full `analysis`. They now agree with both. Anything that
 compensated for the old scale downstream must drop that compensation.
 
+**`analysis` and the direct evaluators now honour `phi_scale`; an unset
+`phi_scale` resolves to `:dft`.** Two halves of one convention had drifted apart.
+
+*`analysis` ignored `phi_scale` entirely.* `synthesis` scales its Fourier bins by
+`phi_inv_scale(cfg)`, but `analysis` applied a fixed `cfg.cphi`, so under `:quad`
+the pair was not mutually inverse: `analysis(cfg, synthesis(cfg, alm))` came back
+as `alm / 2π` exactly. Analysis now applies `cphi · nlon / phi_inv_scale(cfg)`,
+which is `cphi` under `:dft` and restores the inverse property under `:quad`. The
+same factor was threaded through the batch, complex-packed, planned, distributed
+and adjoint analysis paths so they all agree.
+
+*The direct evaluators applied no φ factor at all.* `synthesis_point`,
+`synthesis_point_cplx`, `synthesis_axisym`, `synthesis_axisym_l`, `SH_to_lat`,
+`SH_to_lat_cplx`, `SHqst_to_point`, `SH_to_grad_point`, `SHqst_to_lat` and the
+PencilArray local evaluations each disagreed with the grid they claim to sample by
+exactly 2π under `:quad`. They now carry the same `phi_inv_scale(cfg)/nlon` factor
+`synthesis` does.
+
+*An unset `phi_scale` (`:auto`) used to fall back to a grid-type guess* —
+`grid_type == :gauss ? nlon : nlon/2π` — so a regular grid built through the
+exported `SHTConfig(; …)` keyword constructor disagreed by 2π with the identical
+grid from `create_regular_config`, which sets `:dft` explicitly. Every constructor
+emits `:dft`, so an unset value now means `:dft` too.
+
+**Default `:dft` behaviour is unchanged in all three cases**; only `:quad` and
+hand-built `:auto` configurations move, and they move to the values that make
+`analysis` and `synthesis` inverses.
+
+**Order-mixing rotations now reject `mmax < lmax` instead of truncating.**
+A Wigner-d rotation through a general `β` couples `Y_l^m` to every `Y_l^{m'}` with
+`|m'| ≤ l`. When storage stopped at `mmax < lmax`, the `|m'| > mmax` components
+were silently dropped — measured at `lmax = 8`, that discarded **14.8 %** of the
+field's energy at `mmax = 5` and **24.0 %** at `mmax = 3`, with no error and no
+warning. `SH_Yrotate`, `SH_Yrotate90`, `SH_Xrotate90` and the Euler-angle API now
+raise an `ArgumentError` on such a configuration. Pure Z-rotations are unaffected:
+`β ≡ 0` is diagonal and `β ≡ π` is anti-diagonal (`m' = -m`), so both still work
+at any `mmax`.
+
+*Porting:* use `mmax == lmax` for anything but a Z-rotation. Results that appeared
+to work before were missing the truncated energy.
+
+**Structural `SHTConfig` fields are no longer silently inconsistent.**
+Assigning `cfg.lmax = 10` left `size(cfg.Nlm) == (7, 7)` while the transforms index
+it as `(lmax+1, mmax+1)` under `@inbounds` — an out-of-bounds read of a live array.
+`lmax`, `mmax` and `mres` now rebuild the derived spectral layout (`Nlm`, `nlm`,
+`li`, `mi`, cached scale matrix and m-ordering) and drop the now-stale Legendre
+tables; `nlat`, `nlon`, `grid_type`, `nlm`, `li`, `mi` and `nspat` raise an
+`ArgumentError` pointing at the `create_*_config` constructors, because there is no
+grid-type-independent way to regenerate the quadrature in place.
+
+**The exported `SHTConfig(; …)` keyword constructor validates its invariants.**
+It previously checked nothing, so a hand-built configuration could violate
+`nlon ≥ 2*mmax+1` and then silently synthesize an all-zero field for any mode it
+could not resolve, or hand `use_rfft=true` a raw `BoundsError`. It now enforces the
+same constraints the `create_*_config` helpers always have. The exported keyword
+signature is otherwise unchanged.
+
 ### Fixed
 
+- **`analysis_turbo` / `synthesis_turbo` ignored `mres`.** Both walked a bare
+  `0:mmax` instead of `0:mres:mmax`, so `analysis_turbo` populated — and
+  `synthesis_turbo` consumed — coefficient columns an `mres > 1` transform has no
+  storage for. The disagreement with `analysis`/`synthesis` was O(1), not
+  roundoff. Both now share the core's cached `m` ordering. The turbo pair also
+  nested `@threads :static`, which is illegal inside an outer threaded region;
+  they now fall back to a serial loop there, using the same predicate the core
+  orchestrators use.
+- **Rotation pullbacks read coefficients the primal had already overwritten.**
+  `SH_Yrotate`, `shtns_rotation_apply_cplx` and `shtns_rotation_apply_real`
+  captured their primal *input* and read it lazily, so an in-place call
+  (`Rlm === Qlm`) — or any caller reusing the buffer before the pullback ran —
+  silently corrupted the angle gradient. Each rule now snapshots what it needs at
+  primal time. Applies to both the ChainRules and Zygote adjoints.
+- **`rrule`s for `analysis`/`synthesis`/`analysis_sphtor`/`synthesis_sphtor`
+  declared fewer keyword arguments than their primals.** Passing `use_rfft` or
+  `fft_scratch` — even at its default — made ChainRules skip the rule entirely and
+  fall through to source tracing. The keywords select a different FFT
+  implementation of the same linear operator, so the adjoints are unchanged.
+- **`synthesis_qst` and `analysis_qst` had no `rrule` at all**, so
+  differentiating a QST pipeline fell through to Zygote's source tracing and
+  crashed inside FFTW. Each adjoint is the existing scalar and sphtor adjoints
+  side by side.
+- **`shtns_rotation_apply_real` reported a bare size mismatch for `mres > 1`
+  configurations**, leaving the caller to reverse-engineer why. The message now
+  names `mres` and states the restriction, matching `dist_SH_Yrotate`.
 - **Silent precision loss in batch QST/sphtor transforms.** `analysis_qst_batch`,
   `_synthesis_qst_batch` and the sphtor batch pair derived their output element
   type from one input array instead of promoting across all of them, truncating
@@ -95,6 +178,43 @@ compensated for the old scale downstream must drop that compensation.
 
 ### Internal
 
+- **Removed the dead parallel FFT-plan cache.** The `SHTNSKIT_CACHE_PENCILFFTS`
+  environment variable and the `fft_plan_cache_enabled` / `set_fft_plan_cache!` /
+  `enable_fft_plan_cache!` / `disable_fft_plan_cache!` controls forwarded to a
+  cache in the parallel extension whose only reader, `_get_or_plan`, had no call
+  sites — the "plans" it stored were `NamedTuple` placeholders the FFT wrappers
+  ignored, so every knob was a no-op. The cache the transforms actually use is now
+  in `src/fftutils.jl`, shared by the serial and distributed paths, and the same
+  four controls address it without requiring the extension to be loaded.
+  `SHTNSKIT_FFT_PLAN_CACHE` is the current spelling; the old name still works.
+- **`DistributedSpectralPlan2D` no longer attaches a finalizer.** `close` frees
+  the plan's `l_comm` / `m_comm` sub-communicators, and `MPI_Comm_free` is
+  collective; a finalizer runs at whatever point that rank's garbage collector
+  fires, which is rank-local and nondeterministic. Cleanup is explicit only — call
+  `close(plan)` collectively. Leaking two communicators until `MPI_Finalize` is
+  strictly better than a nondeterministic collective.
+- **`spatial_view(cfg, A)` is exported**, the missing bridge in the padding API:
+  `allocate_padded_spatial` returns an array with `nlat_padded ≥ nlat` rows while
+  every transform requires exactly `nlat`, so the padded buffer could not be passed
+  to `analysis` at all. The view keeps the padded column stride, so it preserves
+  what the padding is for.
+- **`set_batch_size!` is documented as advisory.** `howmany` / `spec_dist` mirror
+  the SHTns C batch descriptors and are stored for interoperability, but the Julia
+  batch entry points take the field count from `size(fields, 3)`; the old docstring
+  claimed otherwise.
+- **Regular-grid quadrature exactness is documented.** Fejér and Clenshaw–Curtis
+  rules with `nlat` nodes are exact only through degree `nlat - 1`, and analysis
+  integrates degree `2*lmax`, so the equiangular grids need `nlat ≥ 2*lmax + 1` —
+  where Gauss–Legendre needs `nlat = lmax + 1`. Below that threshold nothing warns
+  and `analysis ∘ synthesis` is not an identity (7.2e-2 relative error at
+  `lmax = 8, nlat = 10`). Both `create_regular_config` and `docs/src/grids.md` now
+  say so with measured numbers.
+- New regression coverage: planned-vs-`cfg` conformance across the tables / rfft /
+  Robert-form matrix; the `mmax < lmax` rotation guard; angle-gradient survival
+  under an in-place primal; evaluator `phi_scale` agreement; turbo `mres`; and a
+  4-rank `test_mpi_2d_alignment.jl` for the 2D spectral-plan alignment
+  preconditions, wired into CI.
+
 - `pack_lm!`/`pack_lm`/`unpack_lm!`/`unpack_lm` in `src/layout.jl` replace six
   open-coded copies of the packed↔dense `(l,m)` mapping. The `m % mres` guard had
   to be fixed three separate times across those copies.

docs/Distributed_SHTnsKit_Guide.md

@@ -148,7 +148,7 @@ SHTnsKit.dist_synthesis!(spln, fθφ_out, PencilArray(Alm))
 
 Enable plan caching across calls (optional)
 ```bash
-export SHTNSKIT_CACHE_PENCILFFTS=1
+export SHTNSKIT_FFT_PLAN_CACHE=1   # legacy alias: SHTNSKIT_CACHE_PENCILFFTS
 ```
 
 ---
@@ -242,7 +242,10 @@ Both use the distributed transform paths internally and return gradients in the
   ```
 - Robert form: for vector transforms, set `robert_form=true` in your config to stabilize polar behavior.
 - Normalization/phase: match `cfg.norm` and `cfg.cs_phase` to your data; conversions are handled internally on input/output.
-- FFT plan caching: `ENV["SHTNSKIT_CACHE_PENCILFFTS"] = "1"` to reuse PencilFFTs plans.
+- FFT plan caching: `ENV["SHTNSKIT_FFT_PLAN_CACHE"] = "1"` (default) to reuse the
+  per-rank φ-FFT plans. Equivalently `SHTnsKit.enable_fft_plan_cache!()` /
+  `disable_fft_plan_cache!()`. The cache is shared with the serial transforms;
+  cap it with `SHTnsKit.fft_plan_cache_max!(n)`.
 - Y-rotation strategy: truncated gather typically reduces bandwidth; switch to allgather for high-m–dominated spectra.
 
 ---

docs/phi_scaling.md

@@ -20,25 +20,32 @@ The `phi_scale` field in `SHTConfig` controls how the longitude (φ) dimension i
   - Regular grids with poles (simple trapezoidal rule)
 - **Rationale**: Adjusts for the φ integration measure `dφ` where ∫₀²ᵖ f dφ ≈ (2π/nlon) Σ f_j
 
-### `:auto` - Automatic Selection
-- **Behavior**: Chooses based on `grid_type`
-  - `:gauss` → `:dft`
-  - `:driscoll_healy` → `:dft`
-  - `:regular`, `:regular_poles` → `:quad`
+### `:auto` - Unset
+- **Behavior**: treated as `:dft`.
+- Historically this keyed on `grid_type` and handed every non-Gauss grid `:quad`,
+  so a regular grid built through the exported `SHTConfig(; ...)` keyword
+  constructor (which defaulted to `:auto`) disagreed by 2π with the identical
+  grid from `create_regular_config`, which sets `:dft` explicitly. Both
+  constructors emit `:dft`, so an unset value now means `:dft` as well.
 
 ## Configuration
 
 ### In Code
 ```julia
-# Explicit control
-cfg = create_gauss_config(lmax, nlat; phi_scale=:dft)
-cfg = create_regular_config(lmax, nlat; phi_scale=:quad)
+# Every constructor emits :dft; there is no phi_scale keyword on them.
+cfg = create_gauss_config(lmax, nlat)     # :dft
+cfg = create_regular_config(lmax, nlat)   # :dft
 
-# Automatic (recommended)
-cfg = create_gauss_config(lmax, nlat)  # Uses :dft
-cfg = create_regular_config(lmax, nlat)  # Uses :quad
+# To opt into the quadrature convention, set it on the config:
+cfg.phi_scale = :quad
 ```
 
+!!! note
+    An earlier version of this page showed `create_regular_config(lmax, nlat;
+    phi_scale=:quad)` and claimed regular grids default to `:quad`. Neither was
+    true: those constructors take no `phi_scale` keyword (that call raises) and
+    both set `:dft`.
+
 ### Via Environment Variable
 ```bash
 # Override for all grids
@@ -60,18 +67,24 @@ function phi_inv_scale(cfg::SHTConfig)
         return cfg.nlon
     end
 
-    # 2. Use config-specified mode
+    # 2. Use config-specified mode. `:dft` and anything unset (`:auto`) both
+    #    mean the DFT convention — every constructor emits `:dft`, so an unset
+    #    value is not a signal to guess from the grid type.
     if cfg.phi_scale === :quad
         return cfg.nlon / (2π)
-    elseif cfg.phi_scale === :dft
-        return cfg.nlon
+    else
+        return Float64(cfg.nlon)
     end
-
-    # 3. Fall back to grid-type heuristic
-    return cfg.grid_type == :gauss ? cfg.nlon : cfg.nlon / (2π)
 end
 ```
 
+!!! note "The old grid-type fallback is gone"
+    An unset `phi_scale` used to fall back to `grid_type == :gauss ? nlon : nlon/2π`.
+    That made a regular grid built through the exported keyword `SHTConfig(...)`
+    constructor (whose `phi_scale` defaults to `:auto`) disagree by 2π with the
+    identical grid from `create_regular_config`, which sets `:dft` explicitly.
+    Unset now resolves to `:dft`, matching every constructor.
+
 ## Why This Matters
 
 Incorrect φ-scaling leads to round-trip errors:
@@ -90,3 +103,22 @@ The φ-scaling factor must match the quadrature weight convention to ensure:
 - Commit fc1d114: Regular grids changed to quadrature scaling (`nlon/(2π)`)
 - Commit 2441db0: Formalized with auto-detection
 - Current: Explicit `phi_scale` field for clarity and control
+
+
+## Invariant
+
+Whichever mode is selected, `analysis` and `synthesis` are mutual inverses:
+
+```julia
+analysis(cfg, synthesis(cfg, alm)) ≈ alm    # exact under :dft and :quad
+```
+
+`synthesis` scales its Fourier bins by `phi_inv_scale(cfg)` and the inverse FFT
+divides by `nlon`, a net spatial factor `σ`; `analysis` carries `cphi/σ` so the
+two cancel. `:quad` therefore changes the scale of the *spatial* field (by 1/2π)
+without changing what a round trip returns. Point and latitude evaluators
+(`synthesis_point`, `SH_to_lat`, `SHqst_to_lat`, …) apply the same `σ`, so they
+always agree with the grid `synthesis` produces.
+
+Before this was fixed, `analysis` ignored `phi_scale` entirely, so under `:quad`
+a round trip returned `alm/2π` and every evaluator was 2π off from the grid.

docs/src/grids.md

@@ -37,6 +37,26 @@ Gauss–Legendre quadrature with `nlat` points integrates polynomials in
 quadrature is exact at the intended band limit when `nlat = 2*(lmax + 1)`;
 `nlat` must be even.
 
+!!! warning "Equiangular grids need `nlat ≥ 2*lmax + 1`"
+    Fejér and Clenshaw–Curtis rules with `nlat` nodes are exact only through
+    degree ``n_\mathrm{lat} - 1``, and analysis integrates a product of two
+    degree-`lmax` Legendre functions — degree ``2\,l_\mathrm{max}``. So the
+    two regular grids reproduce `alm` from `analysis(synthesis(alm))` only from
+    `nlat = 2*lmax + 1` upward, where Gauss–Legendre needs just `nlat = lmax + 1`.
+
+    Below the threshold nothing warns; the answer is simply wrong. Measured
+    relative round-trip error at `lmax = 8`:
+
+    | `nlat` | `:gauss` | regular midpoint | regular with poles |
+    |---:|---:|---:|---:|
+    | 10 | 8e-16 | 7.2e-2 | 3.6e-1 |
+    | 14 | 8e-16 | 1.3e-2 | 7.1e-3 |
+    | 16 | 8e-16 | 4.2e-3 | 1.3e-3 |
+    | 17 (`2*lmax+1`) | 8e-16 | 7e-16 | 7e-16 |
+
+    Size the grid deliberately, or use Gauss–Legendre (or Driscoll–Healy, exact
+    at `nlat = 2*(lmax+1)`) when round-trip accuracy matters.
+
 !!! tip "Which grid should I choose?"
     Start with Gauss–Legendre unless you need to exchange data with a particular
     equiangular layout. Use regular midpoint for cell-centred data, regular with

ext/ParallelLocal.jl

@@ -80,6 +80,8 @@ end
 
 `Q_p`/`S_p`/`T_p` use the coefficient convention configured by `cfg` (see
 [`dist_SH_to_lat`](@ref)).
+When `cfg.robert_form` is enabled, the tangential components include the
+`sin(θ)` factor used by full-grid vector synthesis; the radial component is unchanged.
 """
 function SHTnsKit.dist_SHqst_to_point(cfg::SHTnsKit.SHTConfig, Q_p::PencilArray, S_p::PencilArray, T_p::PencilArray, cost::Real, phi::Real)
     return SHTnsKit.SHqst_to_point(cfg, Q_p, S_p, T_p, cost, phi)
@@ -91,6 +93,8 @@ end
 
 `Q_p`/`S_p`/`T_p` use the coefficient convention configured by `cfg` (see
 [`dist_SH_to_lat`](@ref)).
+When `cfg.robert_form` is enabled, the tangential components include the
+`sin(θ)` factor used by full-grid vector synthesis; the radial component is unchanged.
 """
 function SHTnsKit.dist_SHqst_to_lat(cfg::SHTnsKit.SHTConfig, Q_p::PencilArray, S_p::PencilArray, T_p::PencilArray, cost::Real;
                                     nphi::Int=cfg.nlon, ltr::Int=cfg.lmax, mtr::Int=cfg.mmax)
@@ -313,6 +317,12 @@ function _pencil_local_qst(cfg, Q::PencilArray, S::PencilArray,
         Vt .*= sinth
         Vp .*= sinth
     end
+    # Same φ convention factor the serial evaluators apply (1 under :dft;
+    # 1/2π under :quad), so a PencilArray input agrees with the dense one.
+    sphi_scale = RT(SHTnsKit._evaluator_phi_scale(cfg))
+    if sphi_scale != 1
+        Vr .*= sphi_scale; Vt .*= sphi_scale; Vp .*= sphi_scale
+    end
     combined = vcat(Vr, Vt, Vp)
     _record_local_payload!(length(combined))
     MPI.Allreduce!(combined, +, comm)
@@ -1135,7 +1145,7 @@ function _analysis_mode_pencil(cfg::SHTnsKit.SHTConfig, im::Int,
     RT = typeof(real(zero(CT)))
     P = Vector{Float64}(undef, ltr + 1)
     rank = MPI.Comm_rank(comm)
-    phi_scale = axisymmetric ? cfg.cphi * cfg.nlon : cfg.cphi
+    phi_scale = axisymmetric ? SHTnsKit._analysis_phi_scale(cfg) * cfg.nlon : SHTnsKit._analysis_phi_scale(cfg)
     for root in 0:(MPI.Comm_size(comm) - 1)
         send = zeros(CT, counts[root + 1])
         @inbounds for (i, θindex) in pairs(θglobals)