diff --git a/docs/src/conventions.md b/docs/src/conventions.md index cc1afd539..ffaadfd29 100644 --- a/docs/src/conventions.md +++ b/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 diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index e3bb3fdd8..32509e77a 100644 --- a/src/ForceFreeStates/Free.jl +++ b/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 diff --git a/src/Utilities/FourierTransforms.jl b/src/Utilities/FourierTransforms.jl index f1037990c..e7a977cbb 100644 --- a/src/Utilities/FourierTransforms.jl +++ b/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 diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index bab0c50f9..a817e9a27 100644 --- a/src/Vacuum/DataTypes.jl +++ b/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 diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index 128616e44..3c3dd6940 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/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 diff --git a/src/Vacuum/Utilities.jl b/src/Vacuum/Utilities.jl index f6a2c17d4..d0d7e04b5 100644 --- a/src/Vacuum/Utilities.jl +++ b/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) diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index b3942b161..b0e008575 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -74,7 +74,8 @@ Green's functions are internal scratch only. # Loop over all decoupled toroidal modes for (idx_n, n) in enumerate(inputs.n_modes) - ft = FourierTransform(inputs.mtheta, mpert, mlow; n=n, ν=plasma_surf.ν) + # On the ϕ = 0 plane, each point sits at ζ = ϕ - ν = -ν + ft = FourierTransform(inputs.mtheta, mpert, mlow; n=n, ζ=(-plasma_surf.ν)) # Diagonal block of wv (and I_v when requested) block_idx = ((idx_n-1)*mpert+1):(idx_n*mpert) @@ -91,7 +92,7 @@ Green's functions are internal scratch only. # Plasma–Plasma block compute_2D_kernel_matrices!(grad_green, green_temp, plasma_surf, plasma_surf, n) - # Project plasma observer onto source basis exp(i*(mθ - nν)) + # Project plasma observer onto source basis exp(i*(mθ - nζ)) mul!(view(grre, 1:num_points_surf, :), green_temp, ft.basis') if !wall.nowall @@ -101,7 +102,7 @@ Green's functions are internal scratch only. compute_2D_kernel_matrices!(grad_green, green_temp, wall, wall, n) # Wall–Plasma block compute_2D_kernel_matrices!(grad_green, green_temp, wall, plasma_surf, n) - # Project wall observer onto source basis exp(i*(mθ - nν)) + # Project wall observer onto source basis exp(i*(mθ - nζ)) mul!(view(grre, (num_points_surf+1):num_points_total, :), green_temp, ft.basis') end @@ -112,32 +113,28 @@ Green's functions are internal scratch only. grad_green_interior = similar!(pool, grad_green) grad_green_interior .= grad_green - # Exterior operator D_ext = 2I + 𝒦 (Chance 1997 eq. 89); the solve gives - # grre = -(2π)²χ^(vo), the vacuum-outside potential. Overwrites the block to save memory. + # Exterior operator D_ext = 2I + 𝒦 (Chance 1997 eq. 89), normals out of the vacuum so δB·n = -b. Green's identity has + # -𝒢·b on the right and the solve uses +𝒢·b, so grre = -(2π)²χ^(vo); (2π)² from GPEC's Jacobian in turns. Overwrites the block. ldiv!(lu!(grad_green), grre) - # Interior operator D_int = D_ext - 2I: the double-layer jump between the two one-sided - # boundary limits is 2I here, giving the vacuum-inside potential grri = χ^(vi). + # Interior: same 𝒦 and right-hand side, diagonal lowered by the double-layer jump 2I, so grri = -(2π)²χ^(vi). for i in 1:num_points_total grad_green_interior[i, i] -= 2.0 end ldiv!(lu!(grad_green_interior), grri) - # Surface-current matrix, Park 2007 eq. 21b: μ₀I^v = χ^(vi) - χ^(vo) = grri - grre - # They are flipped because VACUUM builds the operators in its CW-θ frame while GPEC - # uses CCW-θ, flipping the outward-normal sign. + # Surface-current matrix, Park 2007 eq. 21b: μ₀I^v = χ^(vi) - χ^(vo), so grre - grri = (2π)²μ₀I^v g_diff = @view grri[1:num_points_surf, :] g_diff .= @view(grre[1:num_points_surf, :]) .- g_diff I_v_block = @view vac_data.I_v[block_idx, block_idx] mul!(I_v_block, ft.basis, g_diff) - conj!(I_v_block) # Flip θ_VAC → -θ_VAC to get I^v in GPEC's CCW-θ frame. I_v_block ./= num_points_surf else # Only need exterior system for wv ldiv!(lu!(grad_green), grre) end - # Project exterior kernel onto observer basis exp(-i*(mθ - nν)) and scale to get the response matrix + # Project exterior kernel onto observer basis exp(-i*(mθ - nζ)) and scale to get the response matrix mul!(wv_block, ft.basis, @view(grre[1:num_points_surf, :])) wv_block .*= 4π^2 / num_points_surf end @@ -361,7 +358,7 @@ class operator real and splitting a self-conjugate class into two half-size bloc _solve_projected_rhs!(grre_k, grri_k, Ẽ, op, real_scratch) if compute_Iv - # μ₀Iᵛ = χ^(vi) - χ^(vo) (Park 2007 eq. 21b), accumulated over the blocks + # μ₀Iᵛ = χ^(vi) - χ^(vo) (Park 2007 eq. 21b), so grre - grri = (2π)²μ₀Iᵛ as in the 2D path; accumulated over the blocks g_diff = @view grri_k[1:sz, :] g_diff .= @view(grre_k[1:sz, :]) .- g_diff mul!(Iv_block, Ẽ, g_diff, 1, 1) @@ -372,13 +369,12 @@ class operator real and splitting a self-conjugate class into two half-size bloc end wv_block .*= 4π^2 / num_points_per_fp - # Partner solved the conjugated system: conjugate wv back + # Partner solved the conjugated system: conjugate wv and I_v back partner && conj!(wv_block) @views vac_data.wv[mode_cols, mode_cols] .= wv_block if compute_Iv - # The θ_VAC → -θ_VAC flip on I_v cancels that conjugation for a partner, so only the representative conj!s I_v. - partner || conj!(Iv_block) + partner && conj!(Iv_block) Iv_block ./= num_points_per_fp @views vac_data.I_v[mode_cols, mode_cols] .= Iv_block end diff --git a/test/runtests_coordinate_invariant.jl b/test/runtests_coordinate_invariant.jl index 157c76c54..64300b0dd 100644 --- a/test/runtests_coordinate_invariant.jl +++ b/test/runtests_coordinate_invariant.jl @@ -134,3 +134,24 @@ end @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 + +# N must be the quadratic form ξ†Nξ = ⟨J|ξ|²⟩/(dV/dψ) with N[m, m'] = J_{m−m'}. An up-down asymmetric J makes the +# band complex, so the transposed (conjugate-frame) N misses by O(0.1). +@testset "power_norm_matrix! is the ∮J|ξ|² quadratic form" begin + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + mpert, mlow, ntheta, dV_dpsi = 5, -2, 64, 1.7 + θ = range(; start=0, length=ntheta, step=2π/ntheta) + J = 1 .+ 0.3 .* cos.(θ) .+ 0.2 .* sin.(θ) .+ 0.1 .* sin.(2θ) + # Fourfit storage: jmat[mpert − d] = J_d = ⟨J e^{−idθ}⟩ + jmat = ComplexF64[sum(J .* cis.(-d .* θ)) / ntheta for d in (mpert-1):-1:(-(mpert-1))] + m_modes = mlow:(mlow+mpert-1) + ξ = ComplexF64[cis(0.7m) * (1 + 0.2m) for m in m_modes] + ξθ = [sum(ξ[l] * cis(m * t) for (l, m) in enumerate(m_modes)) for t in θ] + expected = sum(J .* abs2.(ξθ)) / ntheta / dV_dpsi + + N = zeros(ComplexF64, mpert, mpert) + FFS.power_norm_matrix!(N, jmat, mpert, 1, dV_dpsi) + @test N ≈ N' + @test real(ξ' * N * ξ) ≈ expected rtol = 1e-12 + @test abs(ξ' * transpose(N) * ξ - expected) > 0.1 * expected +end diff --git a/test/runtests_fouriertransforms.jl b/test/runtests_fouriertransforms.jl index 005d66d62..b1013d9b5 100644 --- a/test/runtests_fouriertransforms.jl +++ b/test/runtests_fouriertransforms.jl @@ -80,17 +80,17 @@ end @testset "compute_fourier_coefficients" begin atol = ATOL_TIGHT - @testset "2D basis is exp(-i(mθ - nν))" begin + @testset "2D basis is exp(-i(mθ - nζ))" begin N, mlow, mpert = 32, -3, 7 m_modes = mlow:(mlow+mpert-1) n = 2 - ν = collect(range(; start=0.0, length=N, step=0.05)) - basis = compute_fourier_coefficients(N, m_modes, n, ν) + ζ = collect(range(; start=0.0, length=N, step=0.05)) + basis = compute_fourier_coefficients(N, m_modes, n, ζ) @test size(basis) == (mpert, N) θ = collect(range(; start=0, length=N, step=2π/N)) for (l, m) in enumerate(m_modes) - expected = exp.(-im .* (m .* θ .- n .* ν)) + expected = exp.(-im .* (m .* θ .- n .* ζ)) @test all(isapprox.(basis[l, :], expected; atol)) end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 2a5c56883..68b9728d4 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -50,7 +50,7 @@ mtheta_in=5, nzeta_in=1, x=[1.0, 1.1, 1.2, 1.1, 1.0], - z=[0.0, 0.1, 0.0, -0.1, 0.0], + z=[0.0, -0.1, 0.0, 0.1, 0.0], ν=zeros(5), mtheta=5, m_modes=[1], @@ -70,7 +70,7 @@ # Periodic spline requires at least 4 points inputs = VacuumInput(; x=[1.0, 1.1, 1.2, 1.1, 1.0], - z=[0.0, 0.1, 0.0, -0.1, 0.0], + z=[0.0, -0.1, 0.0, 0.1, 0.0], mtheta_in=5, nzeta_in=1, ν=zeros(5), @@ -840,11 +840,10 @@ @testset "compute_vacuum_response 3D I_v matches the 2D path" begin mtheta = 48 θ = range(; start=0, length=mtheta, step=2π / mtheta) - # Up-down asymmetric so Iᵛ is genuinely complex and the θ_VAC → -θ_VAC conjugation is observable + # Up-down asymmetric so Iᵛ is genuinely complex and a frame mismatch between the paths is observable R = 1.7 .+ 0.3 .* cos.(θ) Z = 0.3 .* sin.(θ) .+ 0.08 .* sin.(2θ) .+ 0.08 .* cos.(θ) - # Arrays are reversed for VACUUM's CW θ, as the equilibrium-based constructor does - make(nzeta) = VacuumInput(; x=collect(reverse(R)), z=collect(reverse(Z)), ν=zeros(mtheta), + make(nzeta) = VacuumInput(; x=collect(R), z=collect(Z), ν=zeros(mtheta), mtheta_in=mtheta, nzeta_in=1, m_modes=[-2, -1, 0, 1, 2], n_modes=[1], mtheta=mtheta, nzeta=nzeta) nowall = WallShapeSettings(; shape="nowall") @@ -853,10 +852,69 @@ @test norm(r3d.wv - r2d.wv) / norm(r2d.wv) < 1e-3 @test norm(r3d.I_v - r2d.I_v) / norm(r2d.I_v) < 5e-3 - # The imaginary parts must agree in sign, not be opposed — this is what pins the conjugation + # The imaginary parts must agree in sign, not be opposed — this is what pins the two paths to one frame @test norm(imag.(r3d.I_v) - imag.(r2d.I_v)) < norm(imag.(r3d.I_v) + imag.(r2d.I_v)) end + # ν enters the 3D path only through the geometry (ϕ = ζ + ν) and the 2D path only through the basis (ζ = -ν at ϕ = 0), + # so the two agree only if both use the same relation. Negating ν in 2D alone misses by O(0.1). + @testset "2D and 3D paths agree for a nonzero ν" begin + mtheta = 48 + θ = range(; start=0, length=mtheta, step=2π/mtheta) + R = 1.7 .+ 0.3 .* cos.(θ) + Z = 0.3 .* sin.(θ) .+ 0.08 .* sin.(2θ) .+ 0.08 .* cos.(θ) + ν = 0.25 .* sin.(θ) .+ 0.1 .* cos.(2θ) .+ 0.05 + make(nzeta, ν) = VacuumInput(x=collect(R), z=collect(Z), ν=collect(ν), + mtheta_in=mtheta, nzeta_in=1, m_modes=collect(-3:3), n_modes=[1], mtheta=mtheta, nzeta=nzeta) + nowall = WallShapeSettings(shape="nowall") + + r3d = compute_vacuum_response(make(48, ν), nowall) + @test norm(compute_vacuum_response(make(1, ν), nowall).wv - r3d.wv) / norm(r3d.wv) < 1e-2 + @test norm(compute_vacuum_response(make(1, -ν), nowall).wv - r3d.wv) / norm(r3d.wv) > 1e-1 + end + + # Moving the θ origin by a multiplies every harmonic by exp(i·m·a), so both operators must follow D·A·D† with + # D = diag(exp(i·m·a)). An operator left in the conjugate frame follows D†·A·D instead and fails this by O(0.1). + @testset "wv and I_v follow a θ-origin shift in the same Fourier frame" begin + mtheta, k = 64, 4 + θ = range(; start=0, length=mtheta, step=2π/mtheta) + R = 1.7 .+ 0.5 .* cos.(θ .+ 0.3 .* sin.(θ)) + Z = 0.8 .* sin.(θ) .+ 0.08 .* cos.(θ) .+ 0.1 .* cos.(2θ) + m_modes = collect(-4:6) + D = Diagonal(cis.(m_modes .* (2π * k / mtheta))) + nowall = WallShapeSettings(shape="nowall") + # 2D only: "3D I_v matches the 2D path" ties the 3D frame to this one + make(x, z) = VacuumInput(; x, z, ν=zeros(mtheta), mtheta_in=mtheta, nzeta_in=1, m_modes, n_modes=[1], mtheta, nzeta=1) + ref = compute_vacuum_response(make(collect(R), collect(Z)), nowall; compute_Iv=true) + rolled = compute_vacuum_response(make(circshift(R, -k), circshift(Z, -k)), nowall; compute_Iv=true) + for (A, B) in ((ref.wv, rolled.wv), (ref.I_v, rolled.I_v)) + @test norm(B - D * A * D') / norm(B) < 1e-3 + @test norm(B - D' * A * D) / norm(B) > 5e-2 + end + end + + @testset "clockwise contours" begin + mtheta = 32 + θ = range(; start=0, length=mtheta, step=2π/mtheta) + R, Z = 1.7 .+ 0.3 .* cos.(θ), 0.4 .* sin.(θ) + make(x, z) = VacuumInput(; x, z, ν=zeros(mtheta), mtheta_in=mtheta, nzeta_in=1, m_modes=collect(-2:2), n_modes=[1], mtheta, nzeta=1) + @test_throws ArgumentError compute_vacuum_response(make(reverse(R), reverse(Z)), WallShapeSettings(shape="nowall")) + + # A wall file written clockwise is turned round and gives the same operator as its counter-clockwise twin + Rw, Zw = 1.7 .+ 0.6 .* cos.(θ), 0.8 .* sin.(θ) + wall_response(x, z) = mktempdir() do dir + path = joinpath(dir, "wall.dat") + open(path, "w") do io + println(io, mtheta, "\n", 1.7, "\nindex R Z") + foreach(i -> println(io, i, " ", x[i], " ", z[i]), 1:mtheta) + end + compute_vacuum_response(make(collect(R), collect(Z)), WallShapeSettings(shape=path, equal_arc_wall=false)) + end + ccw = wall_response(Rw, Zw) + cw = @test_logs (:warn, r"runs clockwise") match_mode = :any wall_response(circshift(reverse(Rw), 1), circshift(reverse(Zw), 1)) + @test cw.wv ≈ ccw.wv + end + @testset "compute_vacuum_response 3D conformal wall on a non-axisymmetric boundary" begin # The "3D conformal wall" testset above uses _make_3d_inputs (nzeta_in = 1), so it exercises the # revolved branch. This is the only response-level coverage of the offset branch. @@ -984,7 +1042,7 @@ R = R0 + a * cos(θi) + b * cos(θi - nfp * ζ) + odd * sin(θi - nfp * ζ) push!(X, R * cos(ζ)) push!(Y, R * sin(ζ)) - push!(Z, -a * sin(θi) + b * sin(θi - nfp * ζ) + odd * cos(2θi - nfp * ζ)) + push!(Z, a * sin(θi) - b * sin(θi - nfp * ζ) - odd * cos(2θi - nfp * ζ)) end end return X, Y, Z