Vacuum/ForceFreeStates - BUGFIX! - Put the vacuum and norm matrices in the plasma's Fourier frame

docs/src/conventions.md

@@ -57,9 +57,9 @@ kernel ``\exp(-im\theta)``; the inverse transform reconstructs with ``\exp(+im\t
 
 ### Toroidal Coordinate ``\zeta`` and ``\phi``
 
-- The magnetic coordinate toroidal angle is ``\zeta = \phi/(2\pi) + \nu(\psi,\theta)``, where ``\nu`` is
-  a single-valued straight-field-line offset that depends on the working coordinate. PEST
-  coordinates have ``\nu = 0``.
+- The magnetic coordinate toroidal angle is ``\zeta = (\phi - \nu(\psi,\theta))/(2\pi)``, where ``\phi`` is GPEC's internal
+  toroidal angle (the physical angle is ``-\mathrm{helicity}\,\phi``, below) and ``\nu`` is a single-valued
+  straight-field-line offset that depends on the working coordinate. PEST coordinates have ``\nu = 0``.
 - The physical toroidal angle is reconstructed as
   ``\phi = -\,\mathrm{helicity}\,(2\pi\zeta + \nu)`` (`sample_boundary_grid` in
   `src/ForcingTerms/CoilFourier.jl`). Thus ``\phi`` is effectively **counter-clockwise** (viewed
@@ -196,8 +196,11 @@ conjugate is taken.
 
 ### Interfacing with Vacuum
 
-The Vacuum code uses CCW ``\phi`` and downward-outboard ``\theta``. GPEC uses the complex conjugate
-of RH configurations when interfacing with Vacuum.
+Unlike the Fortran VACUUM code, the GPEC Vacuum module works in the same angles as the rest of the code: ``\theta`` upward-outboard (counter-clockwise in the
+``(R, Z)`` plane) and the Fourier basis ``\exp(-i(m\theta - 2\pi n\zeta))``. A plasma contour passed in clockwise is rejected;
+a clockwise wall file is reversed with a warning. The ``\nu`` passed to Vacuum is GPEC's ``\nu = \phi - 2\pi\zeta``.
+Chance 1997 eq. 97 writes ``\zeta = \phi + \nu`` for the same ``\nu`` because the Fortran VACUUM code's ``\theta`` and ``\phi`` both
+run opposite to GPEC's.
 
 ## Field Amplitudes and Units
 

src/ForceFreeStates/Free.jl

@@ -51,9 +51,9 @@ function power_norm_matrix!(Nmat::AbstractMatrix{ComplexF64}, jmat::AbstractVect
     fill!(Nmat, 0.0 + 0.0im)
     for ipert_n in 1:npert
         off = (ipert_n - 1) * mpert
-        # Toeplitz band index: harmonic difference (m'−m) ∈ [−(mpert−1), mpert−1] maps to 1…2·mpert−1, midpoint mpert = the m=0 Jacobian coefficient
+        # N[m, m'] = J_{m−m'}; jmat stores J_d at index mpert − d, so J_0 sits at mpert
         for ipert_m in 1:mpert, jpert_m in 1:mpert
-            Nmat[off+jpert_m, off+ipert_m] = jmat[jpert_m-ipert_m+mpert] / dV_dpsi
+            Nmat[off+jpert_m, off+ipert_m] = jmat[ipert_m-jpert_m+mpert] / dV_dpsi
         end
     end
     return Nmat

src/Utilities/FourierTransforms.jl

@@ -3,7 +3,7 @@
 
 Pre-computed complex Fourier basis and functor interface for θ ↔ mode transforms.
 
-`basis[ℓ, i] = exp(-i(m_ℓ θ_i - n ν_i))` with shape `(mpert, mtheta)`. Forward: `basis * data / mtheta`;
+`basis[ℓ, i] = exp(-i(m_ℓ θ_i - n ζ_i))` with shape `(mpert, mtheta)`. Forward: `basis * data / mtheta`;
 inverse: `adjoint(basis) * modes` (Fortran `iscdftf`/`iscdftb`; see `docs/src/conventions.md`).
 """
 module FourierTransforms
@@ -14,27 +14,27 @@ export FourierTransform, inverse, inverse_transform!, transform!
 export compute_fourier_coefficients
 
 """
-    compute_fourier_coefficients(mtheta, m_modes, n, ν)
+    compute_fourier_coefficients(mtheta, m_modes, n, ζ)
 
-Build complex basis ``\\exp(-i(m\\theta - n\\nu))`` on the uniform poloidal grid.
+Build complex basis ``\\exp(-i(m\\theta - n\\zeta))`` at the points ``(\\theta_i, \\zeta_i)`` of the uniform poloidal grid.
 
 ## Arguments
 
   - `mtheta`: number of poloidal grid points
   - `m_modes`: poloidal mode numbers (one row per mode)
   - `n`: toroidal mode number
-  - `ν`: toroidal angle offset on the poloidal grid, length `mtheta`
+  - `ζ`: toroidal angle of each grid point in radians, length `mtheta`
 
 ## Returns
 
   - Basis matrix, size `(length(m_modes), mtheta)`
 """
-function compute_fourier_coefficients(mtheta::Int, m_modes::AbstractVector{<:Integer}, n::Integer, ν::Vector{Float64})
+function compute_fourier_coefficients(mtheta::Int, m_modes::AbstractVector{<:Integer}, n::Integer, ζ::Vector{Float64})
 
-    @assert length(ν) == mtheta "ν must have length mtheta"
+    @assert length(ζ) == mtheta "ζ must have length mtheta"
 
     θ_grid = range(; start=0, length=mtheta, step=2π/mtheta)
-    arg = m_modes' .* θ_grid .- n .* ν
+    arg = m_modes' .* θ_grid .- n .* ζ
     return transpose(exp.(-im .* arg))
 end
 
@@ -97,7 +97,7 @@ Struct with precomputed complex Fourier basis for repeated θ ↔ mode transform
   - `mtheta`: poloidal grid size
   - `mpert`: number of poloidal modes
   - `mlow`: lowest poloidal mode number
-  - `basis`: ``\\exp(-i(m\\theta - n\\nu))``, size `(mpert, mtheta)`
+  - `basis`: ``\\exp(-i(m\\theta - n\\zeta))``, size `(mpert, mtheta)`
 """
 struct FourierTransform
     mtheta::Int
@@ -107,7 +107,7 @@ struct FourierTransform
 end
 
 """
-    FourierTransform(mtheta, mpert, mlow; n=0, ν=zeros(mtheta))
+    FourierTransform(mtheta, mpert, mlow; n=0, ζ=zeros(mtheta))
 
 Construct a transform with precomputed basis for contiguous modes `mlow:(mlow+mpert-1)`.
 
@@ -120,7 +120,7 @@ Construct a transform with precomputed basis for contiguous modes `mlow:(mlow+mp
 ## Keyword Arguments
 
   - `n`: toroidal mode number (default 0)
-  - `ν`: toroidal angle offset on the poloidal grid, length `mtheta`
+  - `ζ`: toroidal angle of each grid point in radians, length `mtheta`
 
 ## Returns
 
@@ -131,9 +131,9 @@ function FourierTransform(
     mpert::Int,
     mlow::Int;
     n::Int=0,
-    ν::Vector{Float64}=zeros(Float64, mtheta)
+    ζ::Vector{Float64}=zeros(Float64, mtheta)
 )
-    basis = compute_fourier_coefficients(mtheta, mlow:(mlow+mpert-1), n, ν)
+    basis = compute_fourier_coefficients(mtheta, mlow:(mlow+mpert-1), n, ζ)
     return FourierTransform(mtheta, mpert, mlow, basis)
 end
 

src/Vacuum/DataTypes.jl

@@ -6,13 +6,15 @@ For an axisymmetric boundary, nzeta_in = 1 and only the x and z arrays need to b
 be run with nzeta = 1 for 2D vacuum calculation or nzeta > 1 for 3D vacuum calculation. For a non-axisymmetric boundary,
 nzeta_in > 1 and the x, y, and z arrays need to be provided - the code can then be run with nzeta = 1 for 2D vacuum calculation or
 nzeta > 1 for 3D vacuum calculation. The arrays should be for a single field period only, with excluded endpoints.
+A 2D contour (nzeta_in = 1), or the ζ = 0 cross-section of a 3D boundary, must run counter-clockwise in the (R, Z) plane, GPEC's θ direction; the returned operators are in the
+Fourier frame of the grid passed in.
 
 # Fields
 
   - `x::Vector{Float64}`: Plasma boundary X-coordinate (length mtheta_in * nzeta_in)
   - `y::Vector{Float64}`: Plasma boundary Y-coordinate (length mtheta_in * nzeta_in)
   - `z::Vector{Float64}`: Plasma boundary Z-coordinate (length mtheta_in * nzeta_in)
-  - `ν::Vector{Float64}`: Free parameter in specifying toroidal angle, ζ = ϕ + ν(θ), on input theta grid (axisymmetric only, length mtheta_in)
+  - `ν::Vector{Float64}`: Toroidal offset ν = ϕ - ζ in radians (GPEC's `rzphi_nu`; ζ here is in radians, not turns) on input theta grid (axisymmetric only, length mtheta_in)
   - `mtheta_in::Int`: Number of input poloidal grid points
   - `nzeta_in::Int`: Number of input toroidal grid points (1 for axisymmetric, > 1 for non-axisymmetric)
   - `m_modes::Vector{Int}`: Vector of poloidal mode numbers. E.g. `collect(mlow:mhigh)` for a contiguous range.
@@ -72,13 +74,11 @@ function VacuumInput(
     # Extract plasma surface geometry at this psi
     r, z, ν = extract_plasma_surface_at_psi(equil, ψ)
 
-    # Remove the last point to go from the [0, 2π] grid to VACUUM's [0, 2π) grid
-    # and reverse the arrays for VACUUM's CW θ direction (θ_VAC = -θ_GPEC). This handedness is why
-    # operators returned to GPEC (e.g. the surface-inductance current matrix) are conjugated.
+    # Remove the last point to go from the [0, 2π] grid to the vacuum [0, 2π) grid
     return VacuumInput(;
-        x=reverse(r)[1:(end-1)],
-        z=reverse(z)[1:(end-1)],
-        ν=reverse(ν)[1:(end-1)],
+        x=r[1:(end-1)],
+        z=z[1:(end-1)],
+        ν=ν[1:(end-1)],
         mtheta_in=length(r)-1,
         m_modes=collect(Int, m_modes),
         n_modes=collect(Int, n_modes),
@@ -249,9 +249,9 @@ of length `mtheta`, where `mtheta` is the number of poloidal grid points and θ
 
 # Fields
 
-  - `x::Vector{Float64}`: Plasma surface R-coordinate on VACUUM theta grid
-  - `z::Vector{Float64}`: Plasma surface Z-coordinate on VACUUM theta grid
-  - `ν::Vector{Float64}`: Magnetic toroidal angle offset from geometric toroidal angle
+  - `x::Vector{Float64}`: Plasma surface R-coordinate on the vacuum (counter-clockwise) theta grid
+  - `z::Vector{Float64}`: Plasma surface Z-coordinate on the vacuum (counter-clockwise) theta grid
+  - `ν::Vector{Float64}`: Toroidal offset ν = ϕ - ζ in radians on the vacuum theta grid
 """
 struct PlasmaGeometry
     x::Vector{Float64}
@@ -305,6 +305,7 @@ function PlasmaGeometry(inputs::VacuumInput)
     z = cubic_interp(θ_in, inputs.z, θ_out; bc=PeriodicBC(; endpoint=:exclusive))
     ν = cubic_interp(θ_in, inputs.ν, θ_out; bc=PeriodicBC(; endpoint=:exclusive))
 
+    assert_counterclockwise(x, z, "Plasma boundary")
     return PlasmaGeometry(x, z, ν)
 end
 
@@ -337,53 +338,19 @@ struct PlasmaGeometry3D
 end
 
 """
-    PlasmaGeometry3D(inputs::VacuumInput)
+    PlasmaGeometry3D(inputs::VacuumInput) -> PlasmaGeometry3D
 
-Construct a 3D toroidal plasma surface from vacuum input data.
+Build the 3D plasma surface on the `mtheta × nzeta` grid in (θ, ζ). An axisymmetric boundary (`nzeta_in == 1`) is
+revolved from its 2D contour, placing each point at ϕ = ζ + ν(θ); a 3D boundary is interpolated from its input grid.
+Tangents come from periodic bicubic splines and normals point into the plasma.
 
-This constructor builds a `PlasmaGeometry3D` directly from the `VacuumInput`
-struct, handling both axisymmetric (2D boundary, `nzeta_in == 1`) and fully
-3D input boundaries (`nzeta_in > 1`).
-
-## Axisymmetric input (inputs.nzeta_in == 1)
-
- 1. Build a 2D poloidal contour on the vacuum `mtheta` grid using
-    `PlasmaGeometry(inputs)` to obtain R(theta), Z(theta), and nu(theta).
- 2. Toroidally extrude this contour onto a uniform `nzeta` grid using the
-    SFL angle zeta = phi - nu(theta) and map to Cartesian coordinates:
-    X = R(theta) * cos(zeta - nu(theta)),
-    Y = R(theta) * sin(zeta - nu(theta)),
-    Z = Z(theta).
-
-## Fully 3D input (inputs.nzeta_in > 1)
-
- 1. Interpolate the input (x, y, z) arrays from the original
-    mtheta_in × nzeta_in grid onto the vacuum mtheta × nzeta grid using
-    periodic bicubic interpolation in both angles. The inputs are assumed
-    to already be equally spaced on the SFL angle grid.
-
-## Steps
-
- 1. Fit periodic bicubic splines to each Cartesian component on the
-    (theta, zeta) grid.
- 2. Compute tangent vectors dr/dtheta and dr/dzeta from spline derivatives,
-    scaled by the grid spacings.
- 3. Form oriented normals via the cross product
-    n = (dr/dtheta) × (dr/dzeta) and enforce a consistent orientation
-    (inward for the plasma surface).
- 4. Compute average poloidal/toroidal grid spacings and report the
-    aspect ratio for diagnostics.
-
-## Arguments
+# Arguments
 
-  - `inputs::VacuumInput`: Vacuum calculation inputs defining the boundary
-    geometry and the desired `mtheta, nzeta` resolution.
+  - `inputs::VacuumInput`: Struct containing plasma boundary data and the vacuum grid size
 
-## Returns
+# Returns
 
-  - `PlasmaGeometry3D`: Complete 3D surface description on the
-    `mtheta × nzeta` grid, including points, tangents, normals, and
-    orientation.
+  - `PlasmaGeometry3D`: Struct containing surface points, tangents, and oriented normals
 """
 function PlasmaGeometry3D(inputs::VacuumInput)
 
@@ -405,9 +372,9 @@ function PlasmaGeometry3D(inputs::VacuumInput)
         # Build 3D surface point-by-point from 2D contour
         surf_2D = PlasmaGeometry(inputs)
         for i in 1:mtheta, (j, ζ) in enumerate(ζ_grid)
-            # Our 3D grids are the SFL angle ζ = ϕ - ν
-            r[i+mtheta*(j-1), 1] = surf_2D.x[i] * cos(ζ - surf_2D.ν[i])
-            r[i+mtheta*(j-1), 2] = surf_2D.x[i] * sin(ζ - surf_2D.ν[i])
+            # The grid is in ζ, so each point sits at the geometric angle ϕ = ζ + ν
+            r[i+mtheta*(j-1), 1] = surf_2D.x[i] * cos(ζ + surf_2D.ν[i])
+            r[i+mtheta*(j-1), 2] = surf_2D.x[i] * sin(ζ + surf_2D.ν[i])
             r[i+mtheta*(j-1), 3] = surf_2D.z[i]
         end
     else
@@ -423,6 +390,7 @@ function PlasmaGeometry3D(inputs::VacuumInput)
             itp = cubic_interp((θ_in, ζ_in), reshape(data, inputs.mtheta_in, inputs.nzeta_in); bc=(PeriodicBC(; endpoint=:exclusive), PeriodicBC(; endpoint=:exclusive)))
             r[:, k] = itp(grid_points)
         end
+        assert_counterclockwise(hypot.(r[1:mtheta, 1], r[1:mtheta, 2]), r[1:mtheta, 3], "Plasma boundary ζ = 0 cross-section")
     end
 
     # Compute tangent vectors and normal vectors via periodic bicubic splines
@@ -525,7 +493,7 @@ function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_set
             j = mod1(i - 1, mtheta)
             k = mod1(i + 1, mtheta)
             # Normal vector calculation
-            alph = atan(x_plasma[k] - x_plasma[j], z_plasma[j] - z_plasma[k])
+            alph = atan(x_plasma[j] - x_plasma[k], z_plasma[k] - z_plasma[j])
             x_wall[i] = max(centerstack_min, x_plasma[i] + a * r_minor * cos(alph))
             z_wall[i] = z_plasma[i] + a * r_minor * sin(alph)
         end
@@ -546,25 +514,25 @@ function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_set
         for i in 1:mtheta
             the = (i - 1) * (2π / mtheta)
             x_wall[i] = r_major + a * cos(the)
-            z_wall[i] = -bw_eff * a * sin(the)
+            z_wall[i] = bw_eff * a * sin(the)
         end
 
     elseif wall_settings.shape == "dee"
         wcentr = r_major + cw * r_minor
-        @info "Calculating dee-shaped wall with R = $((@sprintf "%.2e" wcentr)) + $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = -$((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) - $((@sprintf "%.2e" aw)) * $((@sprintf "%.2e" r_minor)) * sin(2θ)."
+        @info "Calculating dee-shaped wall with R = $((@sprintf "%.2e" wcentr)) + $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = $((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) + $((@sprintf "%.2e" aw)) * $((@sprintf "%.2e" r_minor)) * sin(2θ)."
         for i in 1:mtheta
             the = (i - 1) * (2π / mtheta)
             x_wall[i] = wcentr + r_minor * (1.0 + a - cw) * cos(the + dw * sin(the))
-            z_wall[i] = -bw * r_minor * (1.0 + a - cw) * sin(the + tw * sin(2.0*the)) - aw * r_minor * sin(2.0*the)
+            z_wall[i] = bw * r_minor * (1.0 + a - cw) * sin(the + tw * sin(2.0*the)) + aw * r_minor * sin(2.0*the)
         end
 
     elseif wall_settings.shape == "mod_dee"
-        @info "Calculating modified dee-shaped wall with R = $((@sprintf "%.2e" cw)) + $((@sprintf "%.2e" a)) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = -$((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" a)) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) - $((@sprintf "%.2e" aw)) * sin(2θ)."
+        @info "Calculating modified dee-shaped wall with R = $((@sprintf "%.2e" cw)) + $((@sprintf "%.2e" a)) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = $((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" a)) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) + $((@sprintf "%.2e" aw)) * sin(2θ)."
         wcentr = cw
         for i in 1:mtheta
             the = (i - 1) * (2π / mtheta)
             x_wall[i] = cw + a * cos(the + dw * sin(the))
-            z_wall[i] = -bw * a * sin(the + tw * sin(2.0*the)) - aw * sin(2.0*the)
+            z_wall[i] = bw * a * sin(the + tw * sin(2.0*the)) + aw * sin(2.0*the)
         end
 
     else
@@ -592,6 +560,12 @@ function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_set
                 z_wall[i] = parse(Float64, line[3])
             end
         end
+        # Wall files written for the old clockwise convention are turned round, keeping their first point
+        if signed_area(x_wall, z_wall) < 0
+            @warn "Wall file $filepath runs clockwise; reversing its point order to GPEC's counter-clockwise θ."
+            reverse!(@view x_wall[2:end])
+            reverse!(@view z_wall[2:end])
+        end
     end
 
     # Optional: Re-parameterization for equal arc length spacing of wall points
@@ -602,6 +576,7 @@ function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_set
 
     # To add support for x<0 walls, be sure to carefully replicate Chance's fortran code x<0 handling in the kernel function to account for the additional singularities associated with this
     any(x_wall .<= 0.0) && error("Wall R-coordinates contain non-physical values (R <= 0). Check wall geometry.")
+    assert_counterclockwise(x_wall, z_wall, "Wall")
 
     return WallGeometry(false, x_wall, z_wall)
 end

src/Vacuum/Kernel2D.jl

@@ -215,9 +215,9 @@ grad_greenfunction is not zeroed since it fills a different block of the
         end
     end
 
-    # Normals need to point outward from vacuum region. In VACUUM clockwise θ convention, normal points
-    # out of vacuum for wall but inward for plasma, so we multiply by -1 for plasma sources
-    if source isa PlasmaGeometry
+    # Normals need to point outward from the vacuum region. With counter-clockwise θ the normal used in `green`
+    # points into the enclosed region: out of the vacuum for the plasma but into it for the wall, so wall sources flip.
+    if source isa WallGeometry
         grad_greenfunction_block .*= -1
     end
 

src/Vacuum/Utilities.jl

@@ -63,6 +63,27 @@ function extract_plasma_surface_at_psi(equil::Equilibrium.PlasmaEquilibrium, ψ:
     return r, z, ν
 end
 
+"""
+    signed_area(x, z)
+
+Shoelace area of the closed (x, z) contour: positive when it runs counter-clockwise.
+"""
+function signed_area(x::AbstractVector{<:Real}, z::AbstractVector{<:Real})
+    n = length(x)
+    return sum(x[i] * z[mod1(i + 1, n)] - x[mod1(i + 1, n)] * z[i] for i in 1:n) / 2
+end
+
+"""
+    assert_counterclockwise(x, z, name)
+
+Throw unless the closed (x, z) contour runs counter-clockwise. The 2D kernel takes its normals from the
+direction of travel, so a clockwise contour would give a wrong operator silently.
+"""
+function assert_counterclockwise(x::AbstractVector{<:Real}, z::AbstractVector{<:Real}, name::String)
+    signed_area(x, z) > 0 || throw(ArgumentError("$name contour runs clockwise in the (R, Z) plane. The vacuum solve needs it " *
+                        "counter-clockwise, GPEC's θ direction: reverse the point order."))
+end
+
 """
     distribute_to_equal_arc_grid(xin, zin)