Equilibrium - BUGFIX! - Build the root-area weight with the inverse Fourier transform

regression-harness/cases/diiid_n1.toml

@@ -411,3 +411,14 @@ extract = "value"
 label = "Runtime (s)"
 noise_threshold = 0.0
 order = 999
+
+# Root-area-weighted (coordinate-invariant) forcing field b̃ on the control surface: the
+# applied coil spectrum pushed through the √(J|∇ψ|) weight operator. Tracks the field-space
+# conversion that every ResponseMatrices/ field-space output depends on.
+[quantities.forcing_b_rootarea_norm]
+h5path = "PerturbedEquilibrium/forcing_b_root_area"
+type = "complex_vector"
+extract = "norm"
+label = "||forcing b~|| (root-area-weighted)"
+noise_threshold = 1e-10
+order = 1000

src/Equilibrium/CoordinateInvariant.jl

@@ -54,13 +54,14 @@ w(θ) = √(J·|∇ψ|).
 
 Operationally, sqrtamat is the mode-space √weight operator: for a field b with
 Fourier coefficients b_fft, it satisfies the identity
-`‖sqrtamat·b_fft‖² = N² · ∫ |b|² · J|∇ψ| dθ`
-which is Jacobian-invariant on a given flux surface (see
-`scripts/test_power_norm_invariance.jl`).
-
-The backward step uses exp(−imθ)/(1/N) normalization paired with the Julia
-forward FT exp(+imθ) so that round-trip = identity and the convolution
-structure is correct.
+`‖sqrtamat·b_fft‖² = ∫ |b|² · J|∇ψ| dθ`  (θ normalized to [0, 1))
+which is Jacobian-invariant on a given flux surface. Its diagonal is the θ-average of
+√(J|∇ψ|), so `sqrtamat/√jarea` has a diagonal close to (and never above) one; this is the
+matrix Fortran GPEC writes as `J_surf_2` after that division.
+
+Each column is the unit mode e^{+i m_k θ} taken to θ-space with `FourierTransforms.inverse`,
+weighted pointwise, and brought back with the forward transform `ft(...)`, the same pair whose
+round trip is the identity, so the convolution structure is correct.
 """
 function compute_sqrtamat(
     equil::PlasmaEquilibrium,
@@ -78,8 +79,8 @@ function compute_sqrtamat(
         e_k .= 0.0
         e_k[k] = 1.0 + 0.0im
 
-        # Standard backward FT: f(θ_j) = (1/N) Σ_m c_m exp(-imθ_j) = (transpose(basis) * c) / N
-        theta_vec = (transpose(ft.basis) * e_k) ./ mtheta
+        # Inverse FT of the unit mode, f(θ_j) = exp(+i m_k θ_j), through the library's own inverse
+        theta_vec = Utilities.FourierTransforms.inverse(ft, e_k)
 
         # Multiply pointwise by √(J·|∇ψ|) in theta-space
         theta_vec .*= sqrt_jdp

test/runtests_coordinate_invariant.jl

@@ -43,41 +43,41 @@ end
     fm = PE.field_space_response_matrices(Λ, L, P, ϱ, S, jarea)
 
     @testset "Flux recovery contract (round-trip via R = S·A)" begin
-        @test R * fm.permeability / R ≈ P                rtol = 1e-10
-        @test R * fm.surface_inductance * R' ≈ L         rtol = 1e-10
-        @test R * fm.plasma_inductance * R' ≈ Λ          rtol = 1e-10
-        @test (R') \ fm.reluctance / R ≈ ϱ               rtol = 1e-10
+        @test R * fm.permeability / R ≈ P rtol = 1e-10
+        @test R * fm.surface_inductance * R' ≈ L rtol = 1e-10
+        @test R * fm.plasma_inductance * R' ≈ Λ rtol = 1e-10
+        @test (R') \ fm.reluctance / R ≈ ϱ rtol = 1e-10
     end
 
     @testset "Area-weighted (b̄) recovery via S (= flux/A²)" begin
         # b̄-space inductance L_b̄ = S·L̃·S† = L/A² since S·R⁻¹ = A⁻¹·I.
-        @test S * fm.surface_inductance * S' ≈ L ./ jarea^2   rtol = 1e-10
-        @test S * fm.plasma_inductance * S' ≈ Λ ./ jarea^2    rtol = 1e-10
-        @test S * fm.permeability / S ≈ P                     rtol = 1e-10   # similarity: A cancels
+        @test S * fm.surface_inductance * S' ≈ L ./ jarea^2 rtol = 1e-10
+        @test S * fm.plasma_inductance * S' ≈ Λ ./ jarea^2 rtol = 1e-10
+        @test S * fm.permeability / S ≈ P rtol = 1e-10   # similarity: A cancels
     end
 
     @testset "Internal consistency of the b̃ transform rules" begin
-        @test fm.permeability ≈ fm.plasma_inductance / fm.surface_inductance      rtol = 1e-10
+        @test fm.permeability ≈ fm.plasma_inductance / fm.surface_inductance rtol = 1e-10
         L̃inv = inv(fm.surface_inductance)
-        @test fm.reluctance ≈ L̃inv * (fm.plasma_inductance - fm.surface_inductance) * L̃inv  rtol = 1e-10
+        @test fm.reluctance ≈ L̃inv * (fm.plasma_inductance - fm.surface_inductance) * L̃inv rtol = 1e-10
     end
 
     @testset "Energy-scalar invariance (flux ↔ b̃)" begin
         Φ = ComplexF64[cis(0.3k) / k for k in 1:n]
         b̃ = R \ Φ                          # root-area-weighted field
-        @test dot(Φ, L \ Φ) ≈ dot(b̃, fm.surface_inductance \ b̃)   rtol = 1e-10
-        @test dot(Φ, Λ \ Φ) ≈ dot(b̃, fm.plasma_inductance \ b̃)    rtol = 1e-10
-        @test dot(Φ, ϱ * Φ) ≈ dot(b̃, fm.reluctance * b̃)           rtol = 1e-10
+        @test dot(Φ, L \ Φ) ≈ dot(b̃, fm.surface_inductance \ b̃) rtol = 1e-10
+        @test dot(Φ, Λ \ Φ) ≈ dot(b̃, fm.plasma_inductance \ b̃) rtol = 1e-10
+        @test dot(Φ, ϱ * Φ) ≈ dot(b̃, fm.reluctance * b̃) rtol = 1e-10
     end
 
     @testset "Three-field vector relations (b, b̃, b̄; flux = A·b̄)" begin
         b̃ = ComplexF64[cis(0.21k) / (1 + k) for k in 1:n]
         b̄ = S * b̃                          # area-weighted field
         b = (S .* sqrt(jarea)) \ b̃          # bare field b = Σ⁻¹·b̃
         Φ = R * b̃                           # poloidal flux
-        @test Φ ≈ jarea .* b̄               rtol = 1e-12   # Φ = A·b̄
-        @test b̄ ≈ Φ ./ jarea               rtol = 1e-12
-        @test (S .* sqrt(jarea)) * b ≈ b̃    rtol = 1e-12   # Σ·b = b̃
+        @test Φ ≈ jarea .* b̄ rtol = 1e-12   # Φ = A·b̄
+        @test b̄ ≈ Φ ./ jarea rtol = 1e-12
+        @test (S .* sqrt(jarea)) * b ≈ b̃ rtol = 1e-12   # Σ·b = b̃
     end
 end
 
@@ -97,3 +97,40 @@ end
     e_old = sort(real.(eigvals(Mold' * W * Mold)))
     @test e_new ≈ e_old rtol = 1e-12
 end
+
+# The √weight operator on a real flux surface: its diagonal is the θ-average of √(J|∇ψ|), it is
+# Hermitian, it carries the area-weighted field energy exactly (Parseval with the weight), and
+# it inverts against area_to_rootarea_weight. A wrong transform convention (sign of the
+# exponent or a stray 1/N) fails every one of these.
+@testset "sqrtamat on the Solovev surface" begin
+    using TOML
+    EQ = GeneralizedPerturbedEquilibrium.Equilibrium
+    equil_dir = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example")
+    inputs = TOML.parsefile(joinpath(equil_dir, "gpec.toml"))
+    eq_config = EQ.EquilibriumConfig(inputs["Equilibrium"], equil_dir)
+    equil = EQ.setup_equilibrium(eq_config, EQ.SolovevConfig(inputs["SOL_INPUT"]))
+    psi = equil.rzphi_xs[end]
+    mtheta = length(equil.rzphi_ys)
+    mpert, mlow = 21, -8
+    ft = GeneralizedPerturbedEquilibrium.Utilities.FourierTransforms.FourierTransform(mtheta, mpert, mlow)
+    w = EQ.compute_sqrt_jac_delpsi(equil, psi, mtheta)
+    jarea = EQ.flux_surface_area(equil, psi, mtheta)
+    Σ = EQ.compute_sqrtamat(equil, psi, ft)
+    S = EQ.rootarea_to_area_weight(equil, psi, ft)
+    @test size(Σ) == (mpert, mpert)
+    @test all(isapprox.(diag(Σ), sum(w) / mtheta; rtol=1e-12))
+    @test norm(Σ - Σ') / norm(Σ) < 1e-12
+    @test all(0 .< real.(diag(S)) .<= 1)          # Cauchy–Schwarz: ⟨√(J|∇ψ|)⟩ ≤ √⟨J|∇ψ|⟩
+    @test S ≈ Σ ./ sqrt(jarea)
+    # Each column is the forward transform of the weighted unit mode (the definition).
+    b = ComplexF64[cis(0.7k) * (1 + 0.1k) for k in 1:mpert]
+    @test Σ * b ≈ ft(w .* (adjoint(ft.basis) * b)) rtol = 1e-12
+    # Weighted Parseval on the complete basis (mpert = mtheta, no truncation):
+    # ‖Σ·b‖² = (1/N) Σ_j J|∇ψ|_j |f_j|², f the θ-space field of b (θ ∈ [0,1)).
+    ft_full = GeneralizedPerturbedEquilibrium.Utilities.FourierTransforms.FourierTransform(mtheta, mtheta, -(mtheta ÷ 2))
+    Σ_full = EQ.compute_sqrtamat(equil, psi, ft_full)
+    b_full = ComplexF64[cis(0.3k) / (1 + abs(k - mtheta ÷ 2)) for k in 1:mtheta]
+    f_full = adjoint(ft_full.basis) * b_full
+    @test sum(abs2, Σ_full * b_full) ≈ sum(w .^ 2 .* abs2.(f_full)) / mtheta rtol = 1e-10
+    @test S * EQ.area_to_rootarea_weight(equil, psi, ft) ≈ I atol = 1e-10
+end