diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index 68de3c34..d9e8b663 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -48,9 +48,7 @@ ktw = 50.0 # Parameter for ktanh_flag (default: 50.0) ion_flag = true # Include ion dW_k when kin_flag is true electron_flag = false # Include electron dW_k when kin_flag is true -tol_nr = 1e-6 # Relative tolerance of dynamic integration steps away from rationals -tol_r = 1e-7 # Relative tolerance of dynamic integration steps near rationals -crossover = 1e-2 # Fractional distance from rational q at which tol switches +eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Maximum fraction of solutions allowed before re-normalized force_wv_symmetry = true # Forces vacuum energy matrix symmetry diff --git a/src/Utilities/FourierTransforms.jl b/src/Utilities/FourierTransforms.jl index ef018d42..36a53e0b 100644 --- a/src/Utilities/FourierTransforms.jl +++ b/src/Utilities/FourierTransforms.jl @@ -88,41 +88,41 @@ function compute_fourier_coefficients( npert::Int, nlow::Int; n_2D::Union{Nothing, Int}=nothing, - ν::Vector{Float64}=zeros(Float64, mtheta) + ν::Union{Nothing, Vector{Float64}}=nothing ) - # Validate inputs - @assert length(ν) == mtheta "ν must have length mtheta" - @assert !(nzeta > 1 && n_2D !== nothing) "If nzeta > 1, don't set n_2D since we use the full nlow - nhigh range for the toroidal modes" - @assert !(n_2D === nothing && nzeta == 1) "If nzeta == 1 (2D), set n_2D to the toroidal mode number" # Uniform theta grid: [0, 2π) θ_grid = range(; start=0, length=mtheta, step=2π/mtheta) - if !isnothing(n_2D) + if nzeta == 1 + @assert n_2D !== nothing "n_2D must be set for 2D" + @assert ν !== nothing "ν must be set for 2D" + @assert length(ν) == mtheta "ν must have length mtheta" + # In 2D, we only use one toroidal mode at a time # Compute sin(mθ - nν) and cos(mθ - nν) sin_mn_basis = sin.((mlow .+ (0:(mpert-1))') .* θ_grid .- n_2D .* ν) cos_mn_basis = cos.((mlow .+ (0:(mpert-1))') .* θ_grid .- n_2D .* ν) - else # nzeta > 1 + else # 3D + @assert (n_2D === nothing && ν === nothing) "n_2D and ν should be nothing for 3D" + # In 3D, we need to compute the basis for all modes and grid points - # Compute sin(mθ - nν - nϕ) and cos(mθ - nν - nϕ) - ϕ_grid = range(; start=0, length=nzeta, step=2π/nzeta) + # Compute sin(mθ - nζ) and cos(mθ - nζ) + ζ_grid = range(; start=0, length=nzeta, step=2π/nzeta) sin_mn_basis = zeros(mtheta * nzeta, mpert * npert) cos_mn_basis = zeros(mtheta * nzeta, mpert * npert) for idx_n in 1:npert n = nlow + idx_n - 1 n_col_offset = (idx_n - 1) * mpert - # Precompute n*(ν + ϕ) for all (theta, phi) pairs once per n — (mtheta, nzeta) - nν_nϕ = n .* (ν .+ ϕ_grid') for idx_m in 1:mpert m = mlow + idx_m - 1 col = idx_m + n_col_offset - # Broadcast m*θ over all phi; sincos halves trig evaluations - arg = m .* θ_grid .- nν_nϕ # (mtheta, nzeta) - for k in eachindex(arg) - s, c = sincos(arg[k]) - cos_mn_basis[k, col] = c - sin_mn_basis[k, col] = s + for (j, ζ) in enumerate(ζ_grid), (i, θ) in enumerate(θ_grid) + idx = i + (j-1)*mtheta + arg = m * θ - n * ζ + s, c = sincos(arg) + cos_mn_basis[idx, col] = c + sin_mn_basis[idx, col] = s end end end diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index 4d29ba24..63d029b4 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -2,28 +2,38 @@ VacuumInput Struct holding plasma boundary and mode data as provided from ForceFreeStates namelist and computed quantities. +For an axisymmetric boundary, nzeta_in = 1 and only the x 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. 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. # Fields - - `r::Vector{Float64}`: Plasma boundary R-coordinate as a function of poloidal angle - - `z::Vector{Float64}`: Plasma boundary Z-coordinate as a function of poloidal angle - - `ν::Vector{Float64}`: Free parameter in specifying toroidal angle, ϕ = 2πζ + ν(ψ, θ), on theta grid + - `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) + - `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) - `mlow::Int`: Lower poloidal mode number - `mpert::Int`: Number of poloidal modes - `nlow::Int`: Lower toroidal mode number - `npert::Int`: Number of toroidal modes - `mtheta::Int`: Number of vacuum calculation poloidal grid points - - `nzeta::Int`: Number of vacuum calculation toroidal grid points + - `nzeta::Int`: Number of vacuum calculation toroidal grid points (1 for 2D vacuum calculation, > 1 for 3D vacuum calculation) - `force_wv_symmetry::Bool`: Boolean flag to enforce symmetry in the vacuum response matrix """ @kwdef struct VacuumInput - r::Vector{Float64} = Float64[] + x::Vector{Float64} = Float64[] + y::Vector{Float64} = Float64[] z::Vector{Float64} = Float64[] ν::Vector{Float64} = Float64[] - mlow::Int = 0 - mpert::Int = 0 - nlow::Int = 0 - npert::Int = 0 + mtheta_in::Int = 0 + nzeta_in::Int = 1 + mlow::Int = 1 + mpert::Int = 1 + nlow::Int = 1 + npert::Int = 1 mtheta::Int = 1 nzeta::Int = 1 force_wv_symmetry::Bool = true @@ -72,9 +82,10 @@ function VacuumInput( r, z, ν = extract_plasma_surface_at_psi(equil, ψ) return VacuumInput(; - r=reverse(r), + x=reverse(r), z=reverse(z), ν=reverse(ν), + mtheta_in=length(r), mlow=mlow, mpert=mpert, nlow=nlow, @@ -102,11 +113,19 @@ Struct containing input settings for vacuum wall geometry. + `"filepath"`: Custom wall shape from the file you specify - `a::Float64`: Distance of wall from plasma in units of major radius (conformal), or shape parameter (others) + - `aw::Float64`: Half-thickness parameter for Dee-shaped walls + - `bw::Float64`: Elongation parameter for wall shapes + - `cw::Float64`: Offset of the center of the wall from the major radius + - `dw::Float64`: Triangularity parameter for wall shapes + - `tw::Float64`: Sharpness of the corners of the wall (try 0.05 as initial value) + + # Core shape selection + - `equal_arc_wall::Bool`: Flag to enforce equal arc length distribution of nodes on the wall (recommended unless wall is very close to plasma) """ @@ -199,12 +218,15 @@ We interpolate the input plasma boundary arrays from the inputs struct onto the """ function PlasmaGeometry(inputs::VacuumInput) + @assert !isempty(inputs.ν) "ν must be specified for 2D calculations" + @assert length(inputs.x) == length(inputs.z) == length(inputs.ν) "x, z, and ν must have the same length" + # Interpolate arrays from input onto mtheta grid - θ_in = range(0.0, 2π; length=length(inputs.r)) # VacuumInput uses [0, 2π] grid + θ_in = range(0.0, 2π; length=inputs.mtheta_in) # VacuumInput uses [0, 2π] grid θ_out = range(; start=0, length=inputs.mtheta, step=2π/inputs.mtheta) # VACUUM uses [0, 2π) grid # Use one-shot API with PeriodicBC - x = cubic_interp(θ_in, inputs.r, θ_out; bc=PeriodicBC()) # no endpoint handling needed! + x = cubic_interp(θ_in, inputs.x, θ_out; bc=PeriodicBC()) # no endpoint handling needed! z = cubic_interp(θ_in, inputs.z, θ_out; bc=PeriodicBC()) ν = cubic_interp(θ_in, inputs.ν, θ_out; bc=PeriodicBC()) @@ -225,7 +247,6 @@ that the gradient/area elements are scaled by dθ and dζ. - `mtheta::Int`: Number of poloidal grid points - `nzeta::Int`: Number of toroidal grid points - `r::Matrix{Float64}`: Surface points in Cartesian (X,Y,Z), shape (num_points, 3) - - `ν::Vector{Float64}`: Magnetic toroidal angle offset from geometric toroidal angle (same as 2D PlasmaGeometry) - `dr_dθ::Matrix{Float64}`: Poloidal tangent vector ∂r/∂θ × dθ, shape (num_gridpoints, 3) - `dr_dζ::Matrix{Float64}`: Toroidal tangent vector ∂r/∂ζ × dζ, shape (num_points, 3) - `normal::Matrix{Float64}`: Oriented normal vectors, shape (num_points, 3) @@ -235,7 +256,6 @@ struct PlasmaGeometry3D mtheta::Int nzeta::Int r::Matrix{Float64} - ν::Vector{Float64} dr_dθ::Matrix{Float64} dr_dζ::Matrix{Float64} normal::Matrix{Float64} @@ -243,26 +263,53 @@ struct PlasmaGeometry3D end """ - PlasmaGeometry3D(plasma_2d::PlasmaGeometry, nzeta::Int) + PlasmaGeometry3D(inputs::VacuumInput) -Construct a 3D axisymmetric toroidal surface from a 2D poloidal contour. +Construct a 3D toroidal plasma surface from vacuum input data. -# Algorithm +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`). - 0. Interpolate 2D arrays onto mtheta grid - 1. Map 2D (R, Z, ν) to 3D Cartesian: X = R cos(ϕ+ν), Y = R sin(ϕ+ν), Z = Z - 2. Fit periodic bicubic splines to (X, Y, Z) on (θ, ϕ) grid - 3. Compute tangent vectors via spline gradients - 4. Compute normals via cross product: n = ∂r/∂θ × ∂r/∂ζ +## Axisymmetric input (inputs.nzeta_in == 1) -# Arguments + 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). - - `plasma_2d`: 2D poloidal plasma geometry - - `nzeta`: Number of toroidal grid points +## Fully 3D input (inputs.nzeta_in > 1) -# Returns + 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. - - `PlasmaGeometry3D`: Complete 3D surface description +## 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 + + - `inputs::VacuumInput`: Vacuum calculation inputs defining the boundary + geometry and the desired `mtheta, nzeta` resolution. + +## Returns + + - `PlasmaGeometry3D`: Complete 3D surface description on the + `mtheta × nzeta` grid, including points, tangents, normals, and + orientation. """ function PlasmaGeometry3D(inputs::VacuumInput) @@ -272,7 +319,7 @@ function PlasmaGeometry3D(inputs::VacuumInput) dθ = 2π / mtheta dζ = 2π / nzeta θ_grid = range(; start=0, length=mtheta, step=dθ) - ϕ_grid = range(; start=0, length=nzeta, step=dζ) + ζ_grid = range(; start=0, length=nzeta, step=dζ) # Allocate output arrays r = zeros(num_points, 3) @@ -280,18 +327,32 @@ function PlasmaGeometry3D(inputs::VacuumInput) dr_dζ = zeros(num_points, 3) normal = zeros(num_points, 3) - # Interpolate arrays from input onto mtheta grid (same as 2D) - surf_2D = PlasmaGeometry(inputs) - - # Build 3D surface point-by-point from 2D contour - for i in 1:mtheta, (j, ϕ) in enumerate(ϕ_grid) - r[i+mtheta*(j-1), 1] = surf_2D.x[i] * cos(ϕ) - r[i+mtheta*(j-1), 2] = surf_2D.x[i] * sin(ϕ) - r[i+mtheta*(j-1), 3] = surf_2D.z[i] + if inputs.nzeta_in == 1 + # 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]) + r[i+mtheta*(j-1), 3] = surf_2D.z[i] + end + else + # Interpolate 3D inputs onto vacuum grid + @assert !isempty(inputs.y) "y must be provided for 3D (nzeta_in > 1) inputs" + @assert length(inputs.x) == length(inputs.y) == length(inputs.z) == inputs.mtheta_in * inputs.nzeta_in "x, y, z must have length mtheta_in * nzeta_in" + θ_in = range(0.0, 2π; length=inputs.mtheta_in) + ζ_in = range(0.0, 2π; length=inputs.nzeta_in) + θ_flat = repeat(collect(θ_grid), nzeta) + ζ_flat = repeat(collect(ζ_grid); inner=mtheta) + grid_points = (θ_flat, ζ_flat) + for (k, data) in enumerate((inputs.x, inputs.y, inputs.z)) + itp = cubic_interp((θ_in, ζ_in), reshape(data, inputs.mtheta_in, inputs.nzeta_in); bc=(PeriodicBC(), PeriodicBC())) + r[:, k] = itp(grid_points) + end end # Compute tangent vectors and normal vectors via periodic bicubic splines - itps = [cubic_interp((θ_grid, ϕ_grid), reshape(r[:, k], mtheta, nzeta); bc=(PeriodicBC(; endpoint=:exclusive), PeriodicBC(; endpoint=:exclusive))) for k in 1:3] + itps = [cubic_interp((θ_grid, ζ_grid), reshape(r[:, k], mtheta, nzeta); bc=(PeriodicBC(; endpoint=:exclusive), PeriodicBC(; endpoint=:exclusive))) for k in 1:3] for i in 1:mtheta, j in 1:nzeta idx = i + (j - 1) * mtheta for k in 1:3 @@ -318,7 +379,6 @@ function PlasmaGeometry3D(inputs::VacuumInput) mtheta, nzeta, r, - surf_2D.ν, dr_dθ, dr_dζ, normal, @@ -555,6 +615,8 @@ function WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings) ) end + inputs.nzeta_in > 1 && error("3D wall geometry not yet implemented for non-axisymmetric inputs") + # Plasma surface coordinates (2D) surf_2D = PlasmaGeometry(inputs) wall_2D = WallGeometry(inputs, surf_2D, wall_settings) diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 4d7887dc..2f6335d6 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -57,7 +57,9 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC + 3D: `num_modes × num_modes` (full coupled) - `grri`: Interior Green's function matrix. + - `grre`: Exterior Green's function matrix. + - `xzpts`: Coordinate array (mtheta×4 for 2D, mtheta*nzeta×4 for 3D) [R_plasma, Z_plasma, R_wall, Z_wall]. """ @with_pool pool function _compute_vacuum_response_single!( @@ -77,7 +79,8 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC wall = inputs.nzeta > 1 ? WallGeometry3D(inputs, wall_settings) : WallGeometry(inputs, plasma_surf, wall_settings) # Compute Fourier basis coefficients - cos_mn_basis, sin_mn_basis = compute_fourier_coefficients(inputs.mtheta, inputs.mpert, inputs.mlow, inputs.nzeta, inputs.npert, inputs.nlow; n_2D=n_override, ν=plasma_surf.ν) + ν = hasproperty(plasma_surf, :ν) ? plasma_surf.ν : nothing + cos_mn_basis, sin_mn_basis = compute_fourier_coefficients(inputs.mtheta, inputs.mpert, inputs.mlow, inputs.nzeta, inputs.npert, inputs.nlow; n_2D=n_override, ν=ν) num_points_surf, num_modes = size(cos_mn_basis) # Create kernel parameters structs used to dispatch to the correct kernel diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index a8136a1f..debb8299 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -6,13 +6,16 @@ @testset "VacuumInput" begin @testset "default constructor" begin vac = VacuumInput() - @test vac.r == Float64[] + @test vac.x == Float64[] + @test vac.y == Float64[] @test vac.z == Float64[] @test vac.ν == Float64[] - @test vac.mlow == 0 - @test vac.mpert == 0 - @test vac.nlow == 0 - @test vac.npert == 0 + @test vac.mtheta_in == 0 + @test vac.nzeta_in == 1 + @test vac.mlow == 1 + @test vac.mpert == 1 + @test vac.nlow == 1 + @test vac.npert == 1 @test vac.mtheta == 1 @test vac.nzeta == 1 @test vac.force_wv_symmetry == true @@ -50,7 +53,9 @@ @testset "PlasmaGeometry" begin @testset "from VacuumInput" begin inputs = VacuumInput( - r=[1.0, 1.1, 1.2, 1.1, 1.0], + 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], ν=zeros(5), mtheta=5, @@ -72,8 +77,10 @@ @testset "edge: mtheta larger than input length" begin # Periodic spline requires at least 4 points inputs = VacuumInput( - r=[1.0, 1.1, 1.2, 1.1, 1.0], + x=[1.0, 1.1, 1.2, 1.1, 1.0], z=[0.0, 0.1, 0.0, -0.1, 0.0], + mtheta_in=5, + nzeta_in=1, ν=zeros(5), mtheta=8, mpert=1, @@ -92,7 +99,9 @@ # ------------------------------------------------------------------------- @testset "WallGeometry" begin _circle_inputs(mtheta) = VacuumInput( - r=1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta)), + mtheta_in=mtheta, + nzeta_in=1, + x=1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta)), z=0.3 .* sin.(range(0, 2π, length=mtheta)), ν=zeros(mtheta), mtheta=mtheta, @@ -149,7 +158,9 @@ inputs = _circle_inputs(16) plasma_surf_near_zero = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry( VacuumInput( - r=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)), + mtheta_in=16, + nzeta_in=1, + x=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)), z=0.03 .* sin.(range(0, 2π, length=16)), ν=zeros(16), mtheta=16, @@ -170,7 +181,9 @@ inputs = _circle_inputs(16) plasma_surf_near_zero = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry( VacuumInput( - r=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)), + mtheta_in=16, + nzeta_in=1, + x=0.05 .+ 0.03 .* cos.(range(0, 2π, length=16)), z=0.03 .* sin.(range(0, 2π, length=16)), ν=zeros(16), mtheta=16, @@ -360,7 +373,9 @@ # ------------------------------------------------------------------------- @testset "compute_vacuum_response" begin _make_inputs(; mtheta=128, mtheta_eq=17, mpert=2, nlow=1, npert=1) = VacuumInput( - r=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))), + mtheta_in=mtheta_eq, + nzeta_in=1, + x=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))), z=collect(0.3 .* sin.(range(0, 2π, length=mtheta_eq))), ν=zeros(mtheta_eq), mlow=1, @@ -437,7 +452,9 @@ # Kernel requires mtheta, nzeta >= PATCH_DIM (23 for default KernelParams3D(11, 20, 5)) @testset "Vacuum.jl (3D)" begin _make_3d_inputs(; mtheta=32, mtheta_eq=17, mpert=2, nlow=0, npert=2, nzeta=32) = VacuumInput( - r=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))), + mtheta_in=mtheta_eq, + nzeta_in=1, + x=collect(1.7 .+ 0.3 .* cos.(range(0, 2π, length=mtheta_eq))), z=collect(0.3 .* sin.(range(0, 2π, length=mtheta_eq))), ν=zeros(mtheta_eq), mlow=1, @@ -448,6 +465,43 @@ mtheta=mtheta ) + # Helper: simple nonaxisymmetric (3D) plasma boundary built from SFL-style (θ, ζ) coordinates. + _make_3d_nonaxis_inputs(; mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) = begin + θ_in = range(0, 2π, length=mtheta_in) + ζ_in = range(0, 2π, length=nzeta_in) + + X = zeros(mtheta_in, nzeta_in) + Y = zeros(mtheta_in, nzeta_in) + Z = zeros(mtheta_in, nzeta_in) + + R0 = 1.7 + a = 0.3 + ε = 0.05 + + for (i, θ) in enumerate(θ_in), (j, ζ) in enumerate(ζ_in) + # Base circular cross‑section with a small toroidal corrugation + R = R0 + a * cos(θ) + ε * cos(2ζ) * cos(θ) + Z_ij = 0.3 * sin(θ) + ε * sin(2ζ) * sin(θ) + X[i, j] = R * cos(ζ) + Y[i, j] = R * sin(ζ) + Z[i, j] = Z_ij + end + + VacuumInput( + x=vec(X), + y=vec(Y), + z=vec(Z), + mtheta_in=mtheta_in, + nzeta_in=nzeta_in, + mlow=1, + mpert=mpert, + nlow=nlow, + npert=npert, + mtheta=mtheta, + nzeta=nzeta + ) + end + @testset "VacuumInput nzeta > 1" begin vac = VacuumInput(mtheta=32, nzeta=24, mpert=2, npert=2) @test vac.nzeta == 24 @@ -471,8 +525,6 @@ @test size(surf.dr_dθ) == (num_points, 3) @test size(surf.dr_dζ) == (num_points, 3) @test size(surf.normal) == (num_points, 3) - # ν is from 2D contour: one value per poloidal point (length mtheta) - @test length(surf.ν) == inputs.mtheta @test surf.normal_orient in (1, -1) @test all(isfinite, surf.r) @test all(isfinite, surf.normal) @@ -481,6 +533,22 @@ @test !isapprox(surf.r[1, 1], surf.r[1+32, 1]) || !isapprox(surf.r[1, 2], surf.r[1+32, 2]) end + @testset "PlasmaGeometry3D nonaxisymmetric input" begin + inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12) + surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs) + num_points = inputs.mtheta * inputs.nzeta + + @test surf.mtheta == 24 + @test surf.nzeta == 24 + @test size(surf.r) == (num_points, 3) + @test size(surf.normal) == (num_points, 3) + @test surf.normal_orient in (1, -1) + @test all(isfinite, surf.r) + @test all(isfinite, surf.normal) + # Nonaxisymmetric boundary should have genuine 3D structure (non‑trivial Y variation) + @test maximum(abs, surf.r[:, 2]) > 0 + end + @testset "WallGeometry3D nowall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, WallShapeSettings(shape="nowall")) @@ -526,6 +594,27 @@ @test isapprox(wv, wv', rtol=1e-12) end + @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin + inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) + wall_settings = WallShapeSettings(shape="nowall") + wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + + numpoints = inputs.mtheta * inputs.nzeta + num_modes = inputs.mpert * inputs.npert + + @test size(wv) == (num_modes, num_modes) + @test eltype(wv) == ComplexF64 + @test all(isfinite, wv) + @test size(grri) == (2 * numpoints, 2 * num_modes) + @test size(grre) == (2 * numpoints, 2 * num_modes) + @test all(isfinite, grri) + @test all(isfinite, grre) + @test size(plasma_pts) == (numpoints, 3) + @test all(isfinite, plasma_pts) + @test size(wall_pts) == (numpoints, 3) + @test isapprox(wv, wv', rtol=1e-12) + end + @testset "compute_vacuum_response 3D conformal wall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="conformal", a=0.3)