diff --git a/benchmarks/benchmark_fourier_transforms.jl b/benchmarks/benchmark_fourier_transforms.jl deleted file mode 100644 index 8cb150c8..00000000 --- a/benchmarks/benchmark_fourier_transforms.jl +++ /dev/null @@ -1,237 +0,0 @@ -""" -Benchmark script comparing FourierTransform implementations to FFTW. - -This script compares: -1. FourierTransform allocating API: ft(data) -2. FourierTransform in-place API: transform!(output, ft, data) -3. Low-level matrix API: fourier_transform!(gil, gij, cs, m00, l00) -4. FFTW: fft() and rfft() - -The custom FourierTransform is designed for truncated Fourier series with -arbitrary mode ranges (mlow:mhigh), while FFTW computes full DFT (modes 0:N-1). -""" - -using GeneralizedPerturbedEquilibrium -using GeneralizedPerturbedEquilibrium.Utilities.FourierTransforms -using FFTW -using BenchmarkTools -using Printf - -println("="^80) -println("Fourier Transform Benchmark Comparison") -println("="^80) - -# Helper function to extract specific modes from full FFT -function extract_modes(fft_result, mlow, mhigh, mtheta) - """Extract modes mlow:mhigh from full FFT result.""" - modes = zeros(ComplexF64, mhigh - mlow + 1) - for (i, m) in enumerate(mlow:mhigh) - if m >= 0 - # Positive frequencies - modes[i] = fft_result[m+1] / mtheta # FFT normalization - else - # Negative frequencies (wrap around) - modes[i] = fft_result[mtheta+m+1] / mtheta - end - end - return modes -end - -# Test configurations -test_cases = [ - (name="Small (mtheta=128, mpert=10)", mtheta=128, mpert=10, mlow=-5), - (name="Medium (mtheta=256, mpert=20)", mtheta=256, mpert=20, mlow=-10), - (name="Large (mtheta=480, mpert=40)", mtheta=480, mpert=40, mlow=-20), - (name="Very Large (mtheta=1024, mpert=80)", mtheta=1024, mpert=80, mlow=-40) -] - -for test in test_cases - println("\n" * "="^80) - println(test.name) - println("="^80) - - mtheta = test.mtheta - mpert = test.mpert - mlow = test.mlow - mhigh = mlow + mpert - 1 - - # Create test data - theta = range(0, 2π; length=mtheta+1)[1:(end-1)] - data = sin.(3 .* theta) .+ 0.5 .* cos.(7 .* theta) .+ 0.2 .* sin.(11 .* theta) - - # Initialize FourierTransform - ft = FourierTransform(mtheta, mpert, mlow) - - # Pre-allocate buffers for in-place operations - modes_buffer = zeros(ComplexF64, mpert) - theta_buffer = zeros(ComplexF64, mtheta) - - # Pre-allocate for low-level API - cslth, snlth = compute_fourier_coefficients(mtheta, mpert, mlow) - gij = reshape(data, mtheta, 1) # Matrix form - gil = zeros(Float64, mtheta, mpert) - - # Pre-compute FFTW plan for fair comparison - fft_plan = plan_fft(data) - - println("\n--- Forward Transform ---") - - # 1. Our allocating API - print("FourierTransform (allocating): ") - t1 = @benchmark $ft($data) - display(t1) - modes_alloc = ft(data) - - # 2. Our in-place API - print("\nFourierTransform (in-place): ") - t2 = @benchmark transform!($modes_buffer, $ft, $data) - display(t2) - transform!(modes_buffer, ft, data) - - # 3. FFTW - print("\nFFTW (full DFT): ") - t3 = @benchmark $fft_plan * $data - display(t3) - fft_result = fft_plan * data - fft_modes = extract_modes(fft_result, mlow, mhigh, mtheta) - - # Verify results match (for overlapping modes) - # Note: Our transform uses a different normalization and basis - println("\n--- Accuracy Check ---") - println("FourierTransform allocating vs in-place: ", - @sprintf("%.2e", maximum(abs.(modes_alloc .- modes_buffer)))) - - # Compare magnitudes of modes (since basis might differ) - println("Mode magnitudes comparison (FourierTransform vs FFTW):") - println(" FourierTransform peak: ", @sprintf("%.4f", maximum(abs.(modes_alloc)))) - println(" FFTW peak: ", @sprintf("%.4f", maximum(abs.(fft_modes)))) - - println("\n--- Inverse Transform ---") - - # Use modes from our forward transform - modes_test = modes_alloc - - # 1. Our allocating inverse - print("inverse(ft, modes) [allocating]: ") - t4 = @benchmark inverse($ft, $modes_test) - display(t4) - theta_alloc = inverse(ft, modes_test) - - # 2. Our in-place inverse - print("\ninverse_transform! [in-place]: ") - t5 = @benchmark inverse_transform!($theta_buffer, $ft, $modes_test) - display(t5) - inverse_transform!(theta_buffer, ft, modes_test) - - # 3. IFFT - # Note: IFFT requires full spectrum, not just truncated modes - print("\nIFFT (full DFT, not comparable): ") - full_modes = zeros(ComplexF64, mtheta) - for (i, m) in enumerate(mlow:mhigh) - if m >= 0 - full_modes[m+1] = modes_test[i] - else - full_modes[mtheta+m+1] = modes_test[i] - end - end - t6 = @benchmark ifft($full_modes) - display(t6) - - # Accuracy check - println("\n--- Inverse Accuracy Check ---") - println("inverse() allocating vs in-place: ", - @sprintf("%.2e", maximum(abs.(theta_alloc .- theta_buffer)))) - println("Round-trip error (real part): ", - @sprintf("%.2e", maximum(abs.(real.(theta_alloc) .- data)))) - - # Performance summary - println("\n--- Performance Summary ---") - println(@sprintf("Forward transform speedup (in-place vs allocating): %.2fx", - median(t1).time / median(t2).time)) - println(@sprintf("Allocations eliminated: %d → %d", - t1.allocs, t2.allocs)) - - # Compare to FFTW - println(@sprintf("\nFourier vs FFTW (forward): %.2fx %s", - abs(median(t2).time / median(t3).time), - median(t2).time < median(t3).time ? "faster" : "slower")) - println("Note: FFTW computes full DFT (all N modes), we compute truncated series ($mpert modes)") -end - -# Additional test: Matrix transforms (batch processing) -println("\n" * "="^80) -println("Batch Processing Test (Multiple Functions at Once)") -println("="^80) - -mtheta = 256 -mpert = 20 -mlow = -10 -nbatch = 10 # Transform 10 functions simultaneously - -ft = FourierTransform(mtheta, mpert, mlow) -theta = range(0, 2π; length=mtheta+1)[1:(end-1)] - -# Create batch data -data_matrix = zeros(Float64, mtheta, nbatch) -for i in 1:nbatch - data_matrix[:, i] = sin.(i .* theta) -end - -modes_matrix = zeros(ComplexF64, mpert, nbatch) - -println("\nTransforming $nbatch functions of length $mtheta:") - -print("Allocating (loop): ") -@btime for i in 1:($nbatch) - modes = $ft($data_matrix[:, i]) -end - -print("Allocating (matrix):") -@btime $ft($data_matrix) - -print("In-place (loop): ") -modes_buffer = zeros(ComplexF64, mpert) -@btime for i in 1:($nbatch) - transform!($modes_buffer, $ft, $data_matrix[:, i]) -end - -print("In-place (matrix): ") -@btime transform!($modes_matrix, $ft, $data_matrix) - -# Low-level matrix API test -println("\n" * "="^80) -println("Low-Level Matrix API Test (for Vacuum module)") -println("="^80) - -mtheta = 128 -mpert = 10 -mlow = 1 - -# Setup for low-level API -cslth, snlth = compute_fourier_coefficients(mtheta, mpert, mlow) -gij = randn(mtheta, mtheta) # Green's function matrix -gil = zeros(Float64, mtheta, mpert) - -println("\nLow-level fourier_transform! with offset support:") -print("With offsets (m00=0, l00=0): ") -@btime fourier_transform!($gil, $gij, $cslth, 0, 0) - -print("With offsets (m00=64, l00=5): ") -gil_large = zeros(Float64, 256, 20) -@btime fourier_transform!($gil_large, $gij, $cslth, 64, 5) - -println("\n" * "="^80) -println("Benchmark Complete") -println("="^80) - -println("\nKey Takeaways:") -println("1. In-place API eliminates allocations with same performance") -println("2. Custom transforms optimized for truncated series (arbitrary mode ranges)") -println("3. FFTW is faster but computes full DFT (all modes, fixed range 0:N-1)") -println("4. Use FourierTransform when you need:") -println(" - Arbitrary mode ranges (e.g., m = -20:20)") -println(" - Truncated series (fewer modes than grid points)") -println(" - Phase-shifted basis functions (n*qa*delta terms)") -println("5. Use FFTW when you need:") -println(" - Full DFT (all modes 0:N-1)") -println(" - Maximum performance for dense spectral data") diff --git a/src/Equilibrium/CoordinateInvariant.jl b/src/Equilibrium/CoordinateInvariant.jl index fb5c2492..497c365f 100644 --- a/src/Equilibrium/CoordinateInvariant.jl +++ b/src/Equilibrium/CoordinateInvariant.jl @@ -54,7 +54,7 @@ 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θ` +`‖sqrtamat·b_fft‖² = N² · ∫ |b|² · J|∇ψ| dθ` which is Jacobian-invariant on a given flux surface (see `scripts/test_power_norm_invariance.jl`). @@ -78,13 +78,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) - # exp(-imθ) = cos(mθ) - i·sin(mθ), so: - # Re(f) = (1/N)(cslth·Re(c) + snlth·Im(c)) - # Im(f) = (1/N)(cslth·Im(c) - snlth·Re(c)) - real_part = (ft.cslth * real.(e_k) .+ ft.snlth * imag.(e_k)) ./ mtheta - imag_part = (ft.cslth * imag.(e_k) .- ft.snlth * real.(e_k)) ./ mtheta - theta_vec = complex.(real_part, imag_part) + # 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 # Multiply pointwise by √(J·|∇ψ|) in theta-space theta_vec .*= sqrt_jdp diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index c8702e00..f03a7bbb 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -11,8 +11,8 @@ A mutable struct holding data related to the singular surfaces in the equilibriu - `n::Vector{Int}` - Toroidal mode number(s) - `q::Float64` - Safety factor (= m/n) - `q1::Float64` - Derivative of safety factor with respect to ψ - - `grri::Array{Float64,2}` - Interior Green's function at this surface [2*mthvac, 2*mpert] - - `grre::Array{Float64,2}` - Exterior Green's function at this surface [2*mthvac, 2*mpert] + - `grri::Array{ComplexF64,2}` - Interior Green's function at this surface [mthvac, mpert] + - `grre::Array{ComplexF64,2}` - Exterior Green's function at this surface [mthvac, mpert] - `delta_prime::Vector{ComplexF64}` - **STUB (not physically valid)**. Per-surface ca-based Δ' estimate retained for future work / debugging only. The physically valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`, computed via the STRIDE global BVP (Glasser 2018 PoP 25, 032501). Do not use this field for tearing-stability analysis; do not expect agreement with `delta_prime_matrix`. - `delta_prime_col::Matrix{ComplexF64}` - **STUB (not physically valid)**. Per-surface ca-based Δ' column retained for future work / debugging only. Shape (numpert_total × n_res_modes); `delta_prime_col[j, i] = (ca_r[j,ipert_res_i,2] - ca_l[j,ipert_res_i,2]) / (4π²·psio)`. The diagonal element matches the (also stubbed) `delta_prime[i]`. Only populated for the Riccati/parallel FM paths. The physically valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`; this field exists for future development on intra-surface coupling diagnostics, not for production use. """ @@ -23,8 +23,8 @@ A mutable struct holding data related to the singular surfaces in the equilibriu n::Vector{Int} = Int[] q::Float64 = 0.0 q1::Float64 = 0.0 - grri::Array{Float64,2} = Array{Float64}(undef, 0, 0) - grre::Array{Float64,2} = Array{Float64}(undef, 0, 0) + grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 0, 0) + grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 0, 0) delta_prime::Vector{ComplexF64} = ComplexF64[] delta_prime_col::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) ua_left::Array{ComplexF64,3} = Array{ComplexF64}(undef, 0, 0, 0) # asymptotic basis at left inner-layer boundary @@ -413,8 +413,8 @@ Populated in `Free.jl`. - `et::Vector{ComplexF64}` - Total eigenvalues of plasma + vacuum - `n_tor_idx::Vector{Int}` - 0-based toroidal mode number index of each sorted eigenvalue (numpert_total). Needed in `write_imas` - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the vacuum matrix wv, clamped to zero - - `grri::Array{Float64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × 2 * numpert_total) - - `grre::Array{Float64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × 2 * numpert_total) + - `grri::Array{ComplexF64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × numpert_total) + - `grre::Array{ComplexF64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × numpert_total) - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z) - `wall_pts::Array{Float64, 3}` - Cartesian coordinates of wall points, shape (mthvac * nzvac) × 3 for (x, y, z) """ @@ -444,8 +444,8 @@ Populated in `Free.jl`. rootA_wv::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) rootA_wt::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) - grri::Array{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total) - grre::Array{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total) + grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) + grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) plasma_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) wall_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) end diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index a5184af4..51860c30 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -27,7 +27,7 @@ and data dumping. end # Compute vacuum response matrix in-place (handles 2D single-n, 2D multi-n block-diagonal, and 3D) - vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mthvac, ctrl.nzvac, mpert, mlow, npert, nlow; force_wv_symmetry=ctrl.force_wv_symmetry) + vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mthvac, ctrl.nzvac, mlow:mhigh, nlow:nhigh; force_wv_symmetry=ctrl.force_wv_symmetry) Vacuum.compute_vacuum_response!(vac_data, vac_inputs, wall_settings) # Scale by (m - n*q)(m' - n'*q) [Chance Phys. Plasmas 1997 2161 eq. 126] @@ -159,7 +159,7 @@ q-window minimum. ) # Compute raw vacuum matrix at the actual scan psi (singfac NOT applied; free_compute_total applies it analytically) - vac_inputs = Vacuum.VacuumInput(equil, psi_array[i], ctrl.mthvac, ctrl.nzvac, intr.mpert, intr.mlow, intr.npert, intr.nlow; force_wv_symmetry=ctrl.force_wv_symmetry) + vac_inputs = Vacuum.VacuumInput(equil, psi_array[i], ctrl.mthvac, ctrl.nzvac, intr.mlow:intr.mhigh, intr.nlow:intr.nhigh; force_wv_symmetry=ctrl.force_wv_symmetry) wv, _, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) @views wv_array[i, :, :] .= wv end diff --git a/src/ForcingTerms/CoilFourier.jl b/src/ForcingTerms/CoilFourier.jl index ba856f0d..19366a9d 100644 --- a/src/ForcingTerms/CoilFourier.jl +++ b/src/ForcingTerms/CoilFourier.jl @@ -23,16 +23,17 @@ const NZETA_POINTS_PER_PERIOD = 32 Pre-computed plasma boundary grid for evaluating coil fields. ## Fields -- `mtheta`, `nzeta`: grid dimensions -- `R`, `Z`: cylindrical coordinates `[mtheta]` (same for all ζ, axisymmetric) -- `phi_grid`: base toroidal angle grid `[nzeta]` in radians = `-helicity × 2π × j/nzeta` -- `phi_offset`: per-θ toroidal angle correction `[mtheta]` from SFL coordinates: - `ν(ψ, θ)` scaled by `-helicity`, so the physical toroidal angle at `(i, j)` is - `phi_grid[j] + phi_offset[i]`. Matches Fortran's `phi = -helicity*(2π*ζ + dphi(ψ,θ))`. - Zero for axisymmetric equilibria on-axis; non-zero off-axis due to SFL coordinate - transform (Hamada/SFL θ ≠ geometric θ introduces a toroidal offset). -- `dR_dtheta`, `dZ_dtheta`: poloidal derivatives w.r.t. unit-norm angle θ_norm ∈ [0,1] - `dR_dtheta[i] = dR/dθ_norm = 2π × dR/dθ_phys` + + - `mtheta`, `nzeta`: grid dimensions + - `R`, `Z`: cylindrical coordinates `[mtheta]` (same for all ζ, axisymmetric) + - `phi_grid`: base toroidal angle grid `[nzeta]` in radians = `-helicity × 2π × j/nzeta` + - `phi_offset`: per-θ toroidal angle correction `[mtheta]` from SFL coordinates: + `ν(ψ, θ)` scaled by `-helicity`, so the physical toroidal angle at `(i, j)` is + `phi_grid[j] + phi_offset[i]`. Matches Fortran's `phi = -helicity*(2π*ζ + dphi(ψ,θ))`. + Zero for axisymmetric equilibria on-axis; non-zero off-axis due to SFL coordinate + transform (Hamada/SFL θ ≠ geometric θ introduces a toroidal offset). + - `dR_dtheta`, `dZ_dtheta`: poloidal derivatives w.r.t. unit-norm angle θ_norm ∈ [0,1] + `dR_dtheta[i] = dR/dθ_norm = 2π × dR/dθ_phys` """ struct BoundaryGrid mtheta::Int @@ -57,13 +58,13 @@ Poloidal derivatives dR/dθ_norm and dZ/dθ_norm (unit-norm θ_norm ∈ [0,1]) a via periodic cubic splines on the resulting R(θ_norm), Z(θ_norm) data. The toroidal grid direction follows the Fortran GPEC convention: - `phi_j = -helicity × 2π × j/nzeta` where `helicity = sign(Bt) × sign(Ip)`. +`phi_j = -helicity × 2π × j/nzeta` where `helicity = sign(Bt) × sign(Ip)`. This is derived from `equil.params.bt_sign` and `equil.params.crnt`. For DIII-D (Bt < 0, Ip > 0 → helicity = -1): phi increases with j (standard direction). For positive-helicity machines (Bt > 0, Ip > 0 → helicity = +1): phi decreases with j. """ function sample_boundary_grid(equil::Equilibrium.PlasmaEquilibrium, mtheta::Int, nzeta::Int; - psi::Float64=1.0) + psi::Float64=1.0) # Build uniform theta grid (same convention as equil.rzphi_ys, but potentially finer) theta_grid = range(0; length=mtheta, step=1.0/mtheta) @@ -73,7 +74,7 @@ function sample_boundary_grid(equil::Equilibrium.PlasmaEquilibrium, mtheta::Int, for (i, θ_sfl) in enumerate(theta_grid) r_minor = sqrt(equil.rzphi_rsquared((psi, θ_sfl); hint=hint2d)) - θ_cyl = 2π * (θ_sfl + equil.rzphi_offset((psi, θ_sfl); hint=hint2d)) + θ_cyl = 2π * (θ_sfl + equil.rzphi_offset((psi, θ_sfl); hint=hint2d)) R_arr[i] = equil.ro + r_minor * cos(θ_cyl) Z_arr[i] = equil.zo + r_minor * sin(θ_cyl) end @@ -97,7 +98,7 @@ function sample_boundary_grid(equil::Equilibrium.PlasmaEquilibrium, mtheta::Int, bt_sign = !isnothing(equil.params.bt_sign) ? equil.params.bt_sign : 1 ip_sign = !isnothing(equil.params.crnt) ? Int(sign(equil.params.crnt)) : 1 helicity = bt_sign * ip_sign - phi_grid = collect(range(0; length=nzeta, step=-helicity * 2π/nzeta)) + phi_grid = collect(range(0; length=nzeta, step=(-helicity * 2π/nzeta))) # Toroidal angle offset ν(ψ, θ_SFL): in SFL coordinates the physical toroidal angle at # grid point (θ_SFL, ζ_SFL) is φ_phys = -helicity*(2π*ζ_SFL + ν(ψ,θ_SFL)). @@ -119,7 +120,7 @@ Project the cylindrical magnetic field (B_R, B_Z) onto the plasma boundary norma direction ∇ψ and store in `bn[mtheta, nzeta]`. Computes the unit-norm flux element Phi_x per (θ_norm, ζ_norm) cell [T·m²]: - bn(θ_norm, ζ_norm) = 2π × R(θ_norm) × (B_R × ∂Z/∂θ_norm - B_Z × ∂R/∂θ_norm) +bn(θ_norm, ζ_norm) = 2π × R(θ_norm) × (B_R × ∂Z/∂θ_norm - B_Z × ∂R/∂θ_norm) The `2π` factor comes from the toroidal Jacobian ∂r/∂ζ_norm = 2π·R·ê_φ in the cross-product ∂r/∂θ_norm × ∂r/∂ζ_norm. The derivatives dR/dθ_norm and dZ/dθ_norm @@ -129,9 +130,10 @@ are stored in `grid.dR_dtheta` and `grid.dZ_dtheta` (unit-norm convention from The output matches Fortran GPEC's `Phi_x` convention directly (no extra factor needed). ## Arguments -- `bn`: output array `[mtheta, nzeta]`; overwritten in-place -- `B_R`, `B_Z`: cylindrical field components, length `mtheta × nzeta` (flat, θ-major) -- `grid`: pre-computed boundary geometry from `sample_boundary_grid` + + - `bn`: output array `[mtheta, nzeta]`; overwritten in-place + - `B_R`, `B_Z`: cylindrical field components, length `mtheta × nzeta` (flat, θ-major) + - `grid`: pre-computed boundary geometry from `sample_boundary_grid` """ function project_normal_flux!( bn::Matrix{Float64}, @@ -140,7 +142,7 @@ function project_normal_flux!( grid::BoundaryGrid ) mtheta = grid.mtheta - nzeta = grid.nzeta + nzeta = grid.nzeta @assert size(bn) == (mtheta, nzeta) @assert length(B_R) == mtheta * nzeta @assert length(B_Z) == mtheta * nzeta @@ -160,7 +162,7 @@ end mode number `n` and poloidal range `m_low:m_high`. Coefficients are computed as: - bmn = (2 / (mtheta × nzeta)) × Σ_{i,j} bn[i,j] × exp(-i(m×θᵢ - n×ζⱼ)) +bmn = (2 / (mtheta × nzeta)) × Σ_{i,j} bn[i,j] × exp(-i(m×θᵢ - n×ζⱼ)) The factor 2 matches the GPEC/DCON convention for real signals where positive and negative m modes are related by conjugation. @@ -169,7 +171,7 @@ When called after `project_normal_flux!`, the returned amplitudes are in unit-no convention equal to Fortran `Phi_x` (T·m² per unit-norm cell). Uses `compute_fourier_coefficients` from `Utilities.FourierTransforms` with the -3D (mtheta×nzeta, mpert) basis matrix (npert=1, nlow=n). +3D (mpert, mtheta×nzeta) basis matrix (npert=1, nlow=n). """ function fourier_decompose_bn( bn::Matrix{Float64}, @@ -179,23 +181,20 @@ function fourier_decompose_bn( m_high::Int ) mtheta = grid.mtheta - nzeta = grid.nzeta - mpert = m_high - m_low + 1 + nzeta = grid.nzeta - # Build 2D basis: cos(m*θ - n*ζ) and sin(m*θ - n*ζ) - # Using 3D call with npert=1, nlow=n gives shape (mtheta*nzeta, mpert) - cos_basis, sin_basis = compute_fourier_coefficients(mtheta, mpert, m_low, nzeta, 1, n) + # Build Fourier basis: exp(-i*(m*θ - n*ζ)) + # Using 3D call with npert=1, nlow=n gives shape (mpert, mtheta*nzeta) + basis = compute_fourier_coefficients(mtheta, m_low:m_high, nzeta, [n]) bn_flat = vec(bn) # column-major: bn_flat[i + (j-1)*mtheta] = bn[i,j] ✓ scale = 2.0 / (mtheta * nzeta) - bmn_real = scale .* (cos_basis' * bn_flat) - bmn_imag = scale .* (-sin_basis' * bn_flat) + bmn = scale .* (basis * bn_flat) modes = ForcingMode[] for (idx, m) in enumerate(m_low:m_high) - push!(modes, ForcingMode(; n=n, m=m, - amplitude=complex(bmn_real[idx], bmn_imag[idx]))) + push!(modes, ForcingMode(; n=n, m=m, amplitude=bmn[idx])) end return modes end @@ -207,11 +206,12 @@ Top-level entry point: compute Fourier mode amplitudes of the normal magnetic flux from all coil sets on the plasma boundary. ## Pipeline -- Build boundary grid at (mtheta × nzeta) resolution, with helicity from `equil.params` -- Lay out observation points in (R, φ, Z) for all (θ, ζ) combinations -- Run threaded Biot-Savart summation over all coil sets (matches Fortran `field_bs_psi`) -- Project B field onto plasma boundary normal flux (`project_normal_flux!`) -- 2D Fourier decompose to get bmn amplitudes for mode range + + - Build boundary grid at (mtheta × nzeta) resolution, with helicity from `equil.params` + - Lay out observation points in (R, φ, Z) for all (θ, ζ) combinations + - Run threaded Biot-Savart summation over all coil sets (matches Fortran `field_bs_psi`) + - Project B field onto plasma boundary normal flux (`project_normal_flux!`) + - 2D Fourier decompose to get bmn amplitudes for mode range Output amplitudes are in unit-norm convention (= Fortran `Phi_x`). No normalization conversion is needed when using these modes with `compute_plasma_response!`. @@ -226,8 +226,8 @@ function compute_coil_forcing_modes!( n::Int, m_low::Int, m_high::Int; - psi::Float64 = 1.0, - verbose::Bool = false + psi::Float64=1.0, + verbose::Bool=false ) nzeta = cfg.nzeta_coil > 0 ? cfg.nzeta_coil : NZETA_POINTS_PER_PERIOD * max(1, abs(n)) mtheta = cfg.mtheta_coil @@ -237,23 +237,23 @@ function compute_coil_forcing_modes!( grid = sample_boundary_grid(equil, mtheta, nzeta; psi) # Lay out observation points: (theta_i, zeta_j) → cylindrical (R, phi, Z) - nobs = mtheta * nzeta - obs_R = zeros(nobs) + nobs = mtheta * nzeta + obs_R = zeros(nobs) obs_phi = zeros(nobs) - obs_Z = zeros(nobs) + obs_Z = zeros(nobs) for j in 1:nzeta for i in 1:mtheta idx = i + (j - 1) * mtheta - obs_R[idx] = grid.R[i] + obs_R[idx] = grid.R[i] obs_phi[idx] = grid.phi_grid[j] + grid.phi_offset[i] - obs_Z[idx] = grid.Z[i] + obs_Z[idx] = grid.Z[i] end end - B_R = zeros(nobs) + B_R = zeros(nobs) B_phi = zeros(nobs) - B_Z = zeros(nobs) + B_Z = zeros(nobs) compute_biot_savart_boundary!(B_R, B_phi, B_Z, obs_R, obs_phi, obs_Z, coil_sets) verbose && @info " Max |B_R| = $(maximum(abs, B_R)) T, Max |B_Z| = $(maximum(abs, B_Z)) T" @@ -282,6 +282,7 @@ SFL-flux convention. These are scaled by (2π)² to reach unit-norm. If `from_tag == "normal_field_T"`, the amplitudes represent Fourier modes of B·n̂ in Tesla (2π-angle convention). Conversion to unit-norm Phi_x: + - Inverse-Fourier reconstruct B·n̂(θ, ζ) from the input modes - Multiply pointwise by 2π × R(θ) × |dr/dθ_norm(θ)| - Re-Fourier transform to get unit-norm mode amplitudes @@ -290,13 +291,14 @@ The 2π factor comes from the toroidal Jacobian ∂r/∂ζ_norm = 2π·R·ê_φ. This conversion is a mode-mixing operation because R and |dr/dθ| vary poloidally. ## Arguments -- `modes`: `Vector{ForcingMode}` with amplitudes to convert (modified in place) -- `from_tag`: source normalization — `"normal_field_T"` or `"sfl_flux_Wb"` -- `equil`: `PlasmaEquilibrium` providing boundary geometry -- `n`: toroidal mode number -- `m_low`, `m_high`: poloidal mode range (must cover all modes in `modes`) -- `psi`: flux surface for boundary geometry (default 1.0) -- `mtheta`, `nzeta`: grid resolution for the conversion (defaults: 256, 64) + + - `modes`: `Vector{ForcingMode}` with amplitudes to convert (modified in place) + - `from_tag`: source normalization — `"normal_field_T"` or `"sfl_flux_Wb"` + - `equil`: `PlasmaEquilibrium` providing boundary geometry + - `n`: toroidal mode number + - `m_low`, `m_high`: poloidal mode range (must cover all modes in `modes`) + - `psi`: flux surface for boundary geometry (default 1.0) + - `mtheta`, `nzeta`: grid resolution for the conversion (defaults: 256, 64) """ function convert_forcing_normalization!( modes::Vector{ForcingMode}, @@ -305,9 +307,9 @@ function convert_forcing_normalization!( n::Int, m_low::Int, m_high::Int; - psi::Float64 = 1.0, - mtheta::Int = 256, - nzeta::Int = 64 + psi::Float64=1.0, + mtheta::Int=256, + nzeta::Int=64 ) if from_tag == "sfl_flux_Wb" # User provided 2π-angle SFL flux; scale by (2π)² to reach unit-norm (= Phi_x) @@ -325,20 +327,18 @@ function convert_forcing_normalization!( grid = sample_boundary_grid(equil, mtheta, nzeta; psi=psi) # Reconstruct real-space B·n̂(θ, ζ) from input Fourier modes - cos_basis, sin_basis = compute_fourier_coefficients(mtheta, mpert, m_low, nzeta, 1, n) + basis = compute_fourier_coefficients(mtheta, m_low:m_high, nzeta, [n]) - # Build amplitude vectors (real, imag) ordered m_low:m_high - amp_real = zeros(mpert) - amp_imag = zeros(mpert) + # Build amplitude vector ordered m_low:m_high + amp = zeros(ComplexF64, mpert) for mode in modes idx = mode.m - m_low + 1 1 <= idx <= mpert || continue - amp_real[idx] = real(mode.amplitude) - amp_imag[idx] = imag(mode.amplitude) + amp[idx] = mode.amplitude end # Inverse DFT: reconstruct B·n̂(θ, ζ) at grid points - bn_hat = (cos_basis * amp_real .- sin_basis * amp_imag) # length mtheta*nzeta + bn_hat = real.(adjoint(basis) * amp) # length mtheta*nzeta bn_field = reshape(bn_hat, mtheta, nzeta) # Multiply by 2π × R × |dr/dθ_norm| to convert B·n̂ → unit-norm Phi_x integrand. @@ -354,14 +354,13 @@ function convert_forcing_normalization!( # Re-Fourier transform bn_field → unit-norm (Phi_x) mode amplitudes bn_flat = vec(bn_field) scale = 2.0 / (mtheta * nzeta) - new_real = scale .* (cos_basis' * bn_flat) - new_imag = scale .* (-sin_basis' * bn_flat) + bmn = scale .* (basis * bn_flat) # Write back into modes vector (in place, same ordering) for mode in modes idx = mode.m - m_low + 1 1 <= idx <= mpert || continue - mode.amplitude = complex(new_real[idx], new_imag[idx]) + mode.amplitude = bmn[idx] end end diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl index 2272e723..1869053f 100644 --- a/src/ForcingTerms/ForcingTerms.jl +++ b/src/ForcingTerms/ForcingTerms.jl @@ -49,7 +49,9 @@ Data structure for a single forcing mode. ## Fields - `n::Int` - Toroidal mode number + - `m::Int` - Poloidal mode number + - `amplitude::ComplexF64` - Complex amplitude in unit-norm convention (= Fortran Phi_x, T·m² per unit-norm cell) after loading. File inputs are tagged by their input convention: @@ -151,7 +153,7 @@ function load_forcing_ascii!( end end - data = readdlm(filepath; comments=true, comment_char='#') + data = readdlm(filepath; comments=true, comment_char=('#')) nrows = size(data, 1) ncols = size(data, 2) diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index 8b61feac..ba282724 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -37,15 +37,12 @@ function compute_plasma_response!( # Surface inductance L from Green's functions at psilim. # Requires a 2D (nzvac=1) vacuum response so rows are theta points only, - # matching Fortran gpeq_surface which uses grri(2*mthsurf, 2*mpert). - # vac_data.grri has shape [2*mthvac*nzvac, 2*mpert] and cannot be used directly. nn = ffs_intr.nlow - vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, vac_data.mthvac, 1, - ffs_intr.mpert, ffs_intr.mlow, 1, nn) + vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, vac_data.mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn]) wall_nowall = Vacuum.WallShapeSettings(; shape="nowall") _, grri_2d_raw, grre_2d_raw, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall) - grri_2d = Matrix{Float64}(grri_2d_raw) - grre_2d = Matrix{Float64}(grre_2d_raw) + grri_2d = Matrix{ComplexF64}(grri_2d_raw) + grre_2d = Matrix{ComplexF64}(grre_2d_raw) ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac) permeability = calc_permeability(plasma_inductance, surface_inductance) diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index d929e9d3..a4b798a6 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -24,30 +24,36 @@ plasma surface from ForceFreeStates eigenmode solutions. ## What's extracted: -1. **Boundary displacement**: ξ_ψ from `u_store[:, :, 1, end]` - - This is the radial (normal) component of the eigenmode displacement - - At the last radial integration point (plasma edge) - - Dimensions: [numpert_total, numpert_total] + 1. **Boundary displacement**: ξ_ψ from `u_store[:, :, 1, end]` -2. **Flux surface spacing**: dΨ/dρ at boundary - - From equilibrium bicubic spline evaluation - - Needed to convert displacement to magnetic field + + This is the radial (normal) component of the eigenmode displacement + + At the last radial integration point (plasma edge) + + Dimensions: [numpert_total, numpert_total] -3. **Safety factor**: q at boundary - - Used to compute singular factors (m - n*q) - - Identifies resonant surfaces + 2. **Flux surface spacing**: dΨ/dρ at boundary + + + From equilibrium bicubic spline evaluation + + Needed to convert displacement to magnetic field + + 3. **Safety factor**: q at boundary + + + Used to compute singular factors (m - n*q) + + Identifies resonant surfaces ## Arguments -- `equil`: Equilibrium solution containing flux surfaces and q-profile -- `ForceFreeStates_results`: ODE integration results containing u_store with eigenmodes -- `intr`: ForceFreeStates internal state with boundary location (psilim) + + - `equil`: Equilibrium solution containing flux surfaces and q-profile + - `ForceFreeStates_results`: ODE integration results containing u_store with eigenmodes + - `intr`: ForceFreeStates internal state with boundary location (psilim) ## Returns + Named tuple with: -- `ξ_psi_boundary`: Boundary displacement [numpert_total, numpert_total] -- `dPsi_drho`: Flux surface spacing at boundary (scalar) -- `q_boundary`: Safety factor at boundary (scalar) -- `psi_boundary`: Normalized flux at boundary (scalar) + + - `ξ_psi_boundary`: Boundary displacement [numpert_total, numpert_total] + - `dPsi_drho`: Flux surface spacing at boundary (scalar) + - `q_boundary`: Safety factor at boundary (scalar) + - `psi_boundary`: Normalized flux at boundary (scalar) """ function extract_boundary_displacements( equil::Equilibrium.PlasmaEquilibrium, @@ -73,10 +79,10 @@ function extract_boundary_displacements( dPsi_drho = (2π)^2 * equil.psio return ( - ξ_psi_boundary = ξ_psi_boundary, - dPsi_drho = dPsi_drho, - q_boundary = q_boundary, - psi_boundary = psi_boundary + ξ_psi_boundary=ξ_psi_boundary, + dPsi_drho=dPsi_drho, + q_boundary=q_boundary, + psi_boundary=psi_boundary ) end @@ -96,22 +102,25 @@ at the plasma surface. From the ideal MHD constraint [Park Phys. Plasmas 2009 05 where ξ_ψ is the radial displacement eigenfunction. ## Physical Interpretation [Park Phys. Plasmas 2007 052110 Section II]: -- ξ_ψ[i,j]: Displacement of mode i due to eigenmode j -- singfac[i] = m[i] - n*q: Singular factor measuring distance from rational surface -- dΨ/dρ: Converts displacement to flux perturbation (poloidal flux gradient) -- Factor of i: Phase relationship for oscillating fields in complex representation + + - ξ_ψ[i,j]: Displacement of mode i due to eigenmode j + - singfac[i] = m[i] - n*q: Singular factor measuring distance from rational surface + - dΨ/dρ: Converts displacement to flux perturbation (poloidal flux gradient) + - Factor of i: Phase relationship for oscillating fields in complex representation ## Arguments -- `boundary_data`: Output from extract_boundary_displacements() - - ξ_psi_boundary: Boundary displacement [numpert_total, numpert_total] - - dPsi_drho: Flux surface spacing at boundary (scalar) - - q_boundary: Safety factor at boundary (scalar) - - psi_boundary: Normalized flux at boundary (scalar) -- `intr`: ForceFreeStates internal state with mode arrays (mlow, mhigh, nlow, etc.) + + - `boundary_data`: Output from extract_boundary_displacements() + + ξ_psi_boundary: Boundary displacement [numpert_total, numpert_total] + + dPsi_drho: Flux surface spacing at boundary (scalar) + + q_boundary: Safety factor at boundary (scalar) + + psi_boundary: Normalized flux at boundary (scalar) + - `intr`: ForceFreeStates internal state with mode arrays (mlow, mhigh, nlow, etc.) ## Returns -- `bwp_mn[numpert_total, numpert_total]`: Normal magnetic field matrix where bwp_mn[i,j] - is the normal field of Fourier mode i in response to eigenmode j + + - `bwp_mn[numpert_total, numpert_total]`: Normal magnetic field matrix where bwp_mn[i,j] + is the normal field of Fourier mode i in response to eigenmode j """ function compute_normal_magnetic_field( boundary_data::NamedTuple, @@ -163,19 +172,22 @@ In GPEC, this comes from `bwp_mn` (boundary normal field) computed from eigenmod displacements. The flux matrix relates eigenmode displacements to vacuum poloidal flux: -1. Extract eigenmode displacement at plasma boundary from u_store -2. Compute normal magnetic field: B_ψ = i×(dΨ/dρ)×(m - n×q)×ξ_ψ -3. Result is flux[mode_i, eigenmode_j] = bwp_mn[i,j] + + 1. Extract eigenmode displacement at plasma boundary from u_store + 2. Compute normal magnetic field: B_ψ = i×(dΨ/dρ)×(m - n×q)×ξ_ψ + 3. Result is flux[mode_i, eigenmode_j] = bwp_mn[i,j] ## Arguments -- `equil`: Equilibrium solution containing flux surfaces and q-profile -- `ForceFreeStates_results`: ForceFreeStates ODE integration results containing eigenmodes -- `vac_data`: Vacuum response data from free boundary calculation -- `intr`: ForceFreeStates internal state with mode information + + - `equil`: Equilibrium solution containing flux surfaces and q-profile + - `ForceFreeStates_results`: ForceFreeStates ODE integration results containing eigenmodes + - `vac_data`: Vacuum response data from free boundary calculation + - `intr`: ForceFreeStates internal state with mode information ## Returns -- `flxmats[numpert_total, numpert_total]`: Complex flux matrix where flxmats[i,j] is the - vacuum flux of mode i in response to eigenmode j + + - `flxmats[numpert_total, numpert_total]`: Complex flux matrix where flxmats[i,j] is the + vacuum flux of mode i in response to eigenmode j """ function build_flux_matrix( equil::Equilibrium.PlasmaEquilibrium, @@ -214,12 +226,14 @@ vacuum term via the scaling in `free_run!`. The t₁/t₂ factors divide by s_i correctly recovering the properly-normalized inductance. ## Arguments -- `vac_data`: Vacuum data containing wt0 (total energy matrix before eigenvector sorting) -- `ffs_intr`: ForceFreeStates internal state with mode info (mlow, mpert, nlow, qlim) -- `psio`: Total toroidal flux [Wb/rad] from equilibrium (equil.psio) + + - `vac_data`: Vacuum data containing wt0 (total energy matrix before eigenvector sorting) + - `ffs_intr`: ForceFreeStates internal state with mode info (mlow, mpert, nlow, qlim) + - `psio`: Total toroidal flux [Wb/rad] from equilibrium (equil.psio) ## Returns -- Plasma inductance matrix Lambda [numpert_total × numpert_total] + + - Plasma inductance matrix Lambda [numpert_total × numpert_total] """ function calc_plasma_inductance( vac_data::VacuumData, @@ -228,9 +242,9 @@ function calc_plasma_inductance( )::Matrix{ComplexF64} mpert = ffs_intr.numpert_total - chi1 = 2π * psio # = Fortran's chi1 = twopi*psio - n = ffs_intr.nlow - qlim = ffs_intr.qlim # q at psilim + chi1 = 2π * psio # = Fortran's chi1 = twopi*psio + n = ffs_intr.nlow + qlim = ffs_intr.qlim # q at psilim # Singular factors s_i = m_i - n*qlim (same as Fortran: mfac(i) - nn*qlim) s = [((i-1) % ffs_intr.mpert + ffs_intr.mlow) - n * qlim for i in 1:mpert] @@ -245,36 +259,12 @@ function calc_plasma_inductance( for i in 1:mpert, j in 1:mpert t1 = im / (chi1 * s[i] * 2π) t2 = -im / (chi1 * s[j] * 2π) - temp2[i,j] = 2 * t1 * wt0_norm[i,j] * t2 + temp2[i, j] = 2 * t1 * wt0_norm[i, j] * t2 end return inv(temp2) end -""" - pack_complex_grouped!(packed, modes) - -Pack complex mode coefficients into grouped real/imaginary format for Green's function application. - -Converts [a+bi, c+di, ...] to [a, c, ..., b, d, ...] (real block first, then imaginary block). - -This matches the column layout of Julia's Vacuum grri/grre matrices, which store -cos-response columns (1:mpert) followed by sin-response columns (mpert+1:2mpert). -[Vacuum.jl: fourier_transform!(grre, cos_mn_basis) then fourier_transform!(grre, sin_mn_basis; col_offset=mpert)] - -## Arguments -- `packed`: Output array [2*mpert] -- `modes`: Complex Fourier mode coefficients [mpert] -""" -function pack_complex_grouped!(packed::AbstractVector{Float64}, modes::AbstractVector{ComplexF64})::AbstractVector{Float64} - mpert = length(modes) - for i in 1:mpert - packed[i] = real(modes[i]) - packed[i+mpert] = imag(modes[i]) - end - return packed -end - """ calc_permeability( plasma_inductance::Matrix{ComplexF64}, @@ -284,11 +274,13 @@ end Calculate permeability matrix P = Λ·L⁻¹ (matches Fortran `gpresp_permeab`). ## Arguments -- `plasma_inductance`: Plasma inductance matrix Lambda -- `surface_inductance`: Surface inductance matrix L + + - `plasma_inductance`: Plasma inductance matrix Lambda + - `surface_inductance`: Surface inductance matrix L ## Returns -- Permeability matrix P = Lambda * L^{-1} [mpert, mpert] + + - Permeability matrix P = Lambda * L^{-1} [mpert, mpert] """ function calc_permeability( plasma_inductance::Matrix{ComplexF64}, @@ -319,7 +311,7 @@ coordinate-invariant (b̃) matrices in the area-weighted field or recover flux function build_control_surface_rootarea_to_area_weight( equil::Equilibrium.PlasmaEquilibrium, ffs_intr::ForceFreeStatesInternal -)::Tuple{Matrix{ComplexF64}, Float64} +)::Tuple{Matrix{ComplexF64},Float64} mpert = ffs_intr.mpert npert = ffs_intr.npert Npert = ffs_intr.numpert_total @@ -333,7 +325,7 @@ function build_control_surface_rootarea_to_area_weight( S_full = zeros(ComplexF64, Npert, Npert) for in in 1:npert - r = ((in - 1) * mpert + 1):(in * mpert) + r = ((in-1)*mpert+1):(in*mpert) S_full[r, r] .= S_block end return (S_full, jarea) @@ -351,6 +343,7 @@ and the scalar surface area `jarea`. The brief internal flux-conform operator is (`Φ = R·b̃`); poloidal flux never leaves this function. The matrices fall into two algebraic classes: + - **Operators** (map flux → flux): permeability `P` (Φ_tot = P·Φ_x) transforms by similarity `P̃ = R⁻¹·P·R`. Its singular values are coordinate-invariant. - **Quadratic generators** (energy = Φ†·G⁻¹·Φ): inductances `Λ`, `L` transform by congruence @@ -395,11 +388,13 @@ T·m² per unit-norm cell). Files loaded in `normal_field_T` or `sfl_flux_Wb` co are automatically converted to unit-norm on load (see `ForcingMode` docstring). ## Arguments -- `forcing_modes`: External forcing modes (amplitudes in unit-norm / Phi_x convention) -- `intr`: ForceFreeStates internal state with mode arrays + + - `forcing_modes`: External forcing modes (amplitudes in unit-norm / Phi_x convention) + - `intr`: ForceFreeStates internal state with mode arrays ## Returns -- Forcing vector in eigenmode basis [mpert] + + - Forcing vector in eigenmode basis [mpert] """ function map_forcing_to_eigenmodes( forcing_modes::Vector{ForcingMode}, @@ -441,11 +436,13 @@ Compute plasma response to external forcing. Response = Permeability * Forcing ## Arguments -- `permeability`: Permeability matrix -- `forcing_vector`: External forcing in eigenmode basis + + - `permeability`: Permeability matrix + - `forcing_vector`: External forcing in eigenmode basis ## Returns -- Plasma response vector in eigenmode basis + + - Plasma response vector in eigenmode basis """ function compute_plasma_response_vector( permeability::Matrix{ComplexF64}, diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 2a2bec04..eb41cd76 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -86,9 +86,7 @@ function compute_singular_coupling_metrics!( ) ctrl.verbose && @info "Computing singular coupling metrics (GPEC method)" - msing = ffs_intr.msing - mpert = ffs_intr.mpert - numpert_total = ffs_intr.numpert_total + (; msing, numpert_total, mlow, mhigh, nlow, nhigh) = ffs_intr if msing == 0 ctrl.verbose && @info "No singular surfaces found. Skipping singular coupling calculation." @@ -103,9 +101,6 @@ function compute_singular_coupling_metrics!( chi1 = 2π * equil.psio twopi = 2π mtheta = vac_data.mthvac - mlow = ffs_intr.mlow - nlow = ffs_intr.nlow - nhigh = ffs_intr.nhigh wall_settings = Vacuum.WallShapeSettings(; shape="nowall") # Phase 1: Collect all resonant (surface, n) pairs in psi order @@ -115,7 +110,7 @@ function compute_singular_coupling_metrics!( m_res_float = ffs_intr.sing[s].q * nn m_res = round(Int, m_res_float) abs(m_res_float - m_res) > 1e-6 && continue - (m_res < ffs_intr.mlow || m_res > ffs_intr.mhigh) && continue + (m_res < mlow || m_res > mhigh) && continue push!(resonant_pairs, (s, nn)) end end @@ -166,10 +161,10 @@ function compute_singular_coupling_metrics!( end # Compute Green's functions at this surface for this n (once per pair) - vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mpert, mlow, 1, nn) + vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn]) _, grri_raw, grre_raw, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings) - grri = Matrix{Float64}(grri_raw) - grre = Matrix{Float64}(grre_raw) + grri = Matrix{ComplexF64}(grri_raw) + grre = Matrix{ComplexF64}(grre_raw) ffs_intr.sing[s].grri = grri ffs_intr.sing[s].grre = grre @@ -425,8 +420,8 @@ end """ compute_surface_inductance_from_greens( - grri::Matrix{Float64}, - grre::Matrix{Float64}, + grri::Matrix{ComplexF64}, + grre::Matrix{ComplexF64}, ffs_intr::ForceFreeStatesInternal, nn::Int, ν::Vector{Float64} @@ -442,8 +437,8 @@ the DFT (matching Fortran `gpvacuum_flxsurf`'s `EXP(-ifac*nn*dphi)` phase correc ## Arguments - - `grri`: Interior Green's function [2*mtheta, 2*mpert] (GROUPED: cos cols 1:mpert, sin cols mpert+1:2*mpert) - - `grre`: Exterior Green's function [2*mtheta, 2*mpert] + - `grri`: Interior Green's function [mtheta, mpert] + - `grre`: Exterior Green's function [mtheta, mpert] - `ffs_intr`: ForceFreeStates internal state - `nn`: Toroidal mode number - `ν`: Toroidal angle offset on the vacuum theta grid [mtheta] @@ -453,14 +448,14 @@ the DFT (matching Fortran `gpvacuum_flxsurf`'s `EXP(-ifac*nn*dphi)` phase correc Surface inductance matrix [mpert × mpert] """ @with_pool pool function compute_surface_inductance_from_greens( - grri::Matrix{Float64}, - grre::Matrix{Float64}, + grri::Matrix{ComplexF64}, + grre::Matrix{ComplexF64}, ffs_intr::ForceFreeStatesInternal, nn::Int, ν::Vector{Float64} )::Matrix{ComplexF64} mpert = ffs_intr.mpert - mtheta = size(grri, 1) ÷ 2 + mtheta = length(ν) μ₀ = 4π * 1e-7 ft = FourierTransforms.FourierTransform(mtheta, mpert, ffs_intr.mlow) @@ -468,28 +463,21 @@ Surface inductance matrix [mpert × mpert] flux_matrix = zeros!(pool, ComplexF64, mpert, mpert) current_matrix = zeros!(pool, ComplexF64, mpert, mpert) + kax = zeros!(pool, ComplexF64, mtheta) grri_surf = @view grri[1:mtheta, :] grre_surf = @view grre[1:mtheta, :] - kax_re = zeros!(pool, Float64, mtheta) - kax_im = zeros!(pool, Float64, mtheta) - - # Toroidal phase correction: exp(-i*n*ν) matching Fortran gpvacuum_flxsurf - # EXP(-ifac*nn*dphi). Precompute since it's the same for all modes. - cos_nν = cos.(nn .* ν) - sin_nν = sin.(nn .* ν) + # Toroidal phase correction: exp(-i*n*ν) + phase = cis.(-nn .* ν) for i in 1:mpert flux_matrix[i, i] = 1.0 - for k in 1:mtheta - kax_re[k] = (grri_surf[k, i] + grre_surf[k, i]) / (μ₀ * (2π)^2) - kax_im[k] = (grri_surf[k, mpert+i] + grre_surf[k, mpert+i]) / (μ₀ * (2π)^2) - end + # Complex grri/e stores exp(i(mθ-nν)) projection, need conjugate for exp(-i(mθ-nν)) + kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ₀ * (2π)^2) - # Port of Fortran gpvacuum_flxsurf: apply the toroidal phase exp(-i*n*ν), - # reverse theta (VACUUM theta runs opposite to DCON theta), then forward-DFT. - g_phased = (kax_re .- im .* kax_im) .* (cos_nν .- im .* sin_nν) + # Port of Fortran gpvacuum_flxsurf: apply toroidal phase, reverse theta, forward-DFT. + g_phased = kax .* phase current_matrix[:, i] = ft(_reverse_theta(g_phased)) end diff --git a/src/Utilities/FourierTransforms.jl b/src/Utilities/FourierTransforms.jl index 47703882..30ebdb70 100644 --- a/src/Utilities/FourierTransforms.jl +++ b/src/Utilities/FourierTransforms.jl @@ -1,210 +1,130 @@ """ FourierTransforms -Utility module for efficient Fourier transforms using pre-computed basis functions. +Pre-computed complex Fourier basis and functor interface for θ ↔ mode transforms. -Provides a functor-based interface for Fourier transforms between theta-space and mode-space -representations. Supports both simple harmonic transforms (cos(m*θ), sin(m*θ)) and phase-shifted -transforms (cos(m*θ + n*qa*δ), sin(m*θ + n*qa*δ)) used in vacuum field calculations. - -# Features - -- Pre-computes trigonometric basis functions for efficiency -- Direct complex number support (no manual real/imaginary splitting) -- Type-stable functor pattern for high performance -- Supports different grids (mthvac, mtheta) with different mode ranges +`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 -# Example +using LinearAlgebra -```julia -using GeneralizedPerturbedEquilibrium.Utilities.FourierTransforms +export FourierTransform, inverse +export compute_fourier_coefficients -# For PerturbedEquilibrium (no phase shift) -ft = FourierTransform(mtheta, mpert, mlow) +""" + compute_fourier_coefficients(mtheta, m_modes, n, ν) -# For Vacuum (with n*qa*delta phase) -ft_vac = FourierTransform(mthvac, mpert, mlow; n=1, qa=2.5, delta=delta_array) +Build complex basis ``\\exp(-i(m\\theta - n\\nu))`` on the uniform poloidal grid. -# Forward transform: theta-space → Fourier modes -theta_data = randn(mtheta, 10) # 10 different theta-space functions -modes = ft(theta_data) # Complex modes [mpert, 10] +## Arguments -# Inverse transform: Fourier modes → theta-space -reconstructed = inverse(ft, modes) -``` + - `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` -# Sign convention +## Returns -The forward transform uses `exp(-imθ)` (Fortran `iscdftf`) for both real- and -complex-valued input; `inverse` uses `exp(+imθ)` (Fortran `iscdftb`). Hence -`inverse(ft, ft(data)) ≈ data` for any real- or complex-valued `data`. + - Basis matrix, size `(length(m_modes), mtheta)` """ -module FourierTransforms +function compute_fourier_coefficients(mtheta::Int, m_modes::AbstractVector{<:Integer}, n::Integer, ν::Vector{Float64}) -using LinearAlgebra + @assert length(ν) == mtheta "ν must have length mtheta" -export FourierTransform, inverse -export compute_fourier_coefficients -export transform!, inverse_transform! -export fourier_transform!, fourier_inverse_transform! + θ_grid = range(; start=0, length=mtheta, step=2π/mtheta) + arg = m_modes' .* θ_grid .- n .* ν + return transpose(exp.(-im .* arg)) +end """ - compute_fourier_coefficients(mtheta, mpert, mlow, nzeta, npert, nlow; n=nothing, ν=zeros(mtheta)) - -Compute Fourier basis function coefficients for transforms between physical-space and mode-space. -Supports both 2D and 3D geometries. In 2D, we only use one toroidal mode at a time, so we can -just use the n argument. In 3D, we need to compute the basis for all modes and grid points. - -# Arguments + compute_fourier_coefficients(mtheta, m_modes, nzeta, n_modes; nfp=1) -- `mtheta::Int`: Number of poloidal grid points (theta resolution) -- `mpert::Int`: Number of Fourier modes (spectral resolution) -- `mlow::Int`: Lowest mode number (mode numbering starts here) -- `nzeta::Int`: Number of toroidal grid points -- `npert::Int`: Number of toroidal modes -- `nlow::Int`: Lowest toroidal mode number +Build 3D basis for every `(m, n)` in `m_modes × n_modes` (rows m-fast, n-slow). -# Keyword Arguments +## Arguments -- `n_2D::Union{Nothing, Int}=nothing`: Toroidal mode number for 2D (default: nothing) -- `ν::Vector{Float64}=zeros(mtheta)`: Toroidal angle offset array (default: no offset, n*ν = 0) + - `mtheta`: poloidal grid points per toroidal plane + - `m_modes`: poloidal mode numbers + - `nzeta`: full-torus toroidal grid points (must be divisible by `nfp`) + - `n_modes`: toroidal mode numbers -# Returns +## Keyword Arguments -- 2D - - `cos_mn_basis::Matrix{Float64}`: Cosine coefficients `cos(m*θ - n*ν)` [mtheta, mpert] - - `sin_mn_basis::Matrix{Float64}`: Sine coefficients `sin(m*θ - n*ν)` [mtheta, mpert] -- 3D - - `cos_mn_basis::Matrix{Float64}`: Cosine coefficients `cos(m*θ - n*ν - n*ϕ)` [mtheta * nzeta, mpert * npert] - - `sin_mn_basis::Matrix{Float64}`: Sine coefficients `sin(m*θ - n*ν - n*ϕ)` [mtheta * nzeta, mpert * npert] + - `nfp`: number of field periods; when `> 1`, emits one period only (`nzeta ÷ nfp` planes) -# Notes +## Returns -The theta and phi grids are uniform: `θᵢ = 2π*i/mtheta` for `i = 0:mtheta-1` and `ϕⱼ = 2π*j/nzeta` for `j = 0:nzeta-1` - -When `n=0, ν=0` (default), this reduces to simple harmonic basis: -- `cos_mn_basis[i,l] = cos(m*θᵢ)` -- `sin_mn_basis[i,l] = sin(m*θᵢ)` + - Basis matrix, size `(length(m_modes) * length(n_modes), mtheta * (nzeta ÷ nfp))` """ function compute_fourier_coefficients( mtheta::Int, - mpert::Int, - mlow::Int, + m_modes::AbstractVector{<:Integer}, nzeta::Int, - npert::Int, - nlow::Int; - n_2D::Union{Nothing, Int}=nothing, - ν::Union{Nothing, Vector{Float64}}=nothing + n_modes::AbstractVector{<:Integer}; + nfp::Int=1 ) + @assert nzeta % nfp == 0 "nzeta ($nzeta) must be divisible by nfp ($nfp)" - # Uniform theta grid: [0, 2π) θ_grid = range(; start=0, length=mtheta, step=2π/mtheta) - - 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 # 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ζ) 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 - for idx_m in 1:mpert - m = mlow + idx_m - 1 - col = idx_m + n_col_offset - 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 + ζ_grid = range(; start=0, length=nzeta, step=2π/nzeta) + nzeta_out = nzeta ÷ nfp + + mpert = length(m_modes) + basis = zeros(ComplexF64, mpert * length(n_modes), mtheta * nzeta_out) + for (idx_n, n) in enumerate(n_modes) + n_row_offset = (idx_n - 1) * mpert + for (idx_m, m) in enumerate(m_modes) + row = idx_m + n_row_offset + for j in 1:nzeta_out, (i, θ) in enumerate(θ_grid) + col = i + (j - 1) * mtheta + basis[row, col] = cis(-(m * θ - n * ζ_grid[j])) end end end - return cos_mn_basis, sin_mn_basis + return basis end """ FourierTransform -Callable struct for efficient Fourier transforms with pre-computed basis functions. +Struct with precomputed complex Fourier basis for repeated θ ↔ mode transforms. # Fields -- `mtheta::Int`: Number of theta grid points -- `mpert::Int`: Number of Fourier modes -- `mlow::Int`: Lowest mode number -- `cslth::Matrix{Float64}`: Cosine basis functions [mtheta, mpert] -- `snlth::Matrix{Float64}`: Sine basis functions [mtheta, mpert] - -# Usage - -Create once with appropriate grid parameters, then use repeatedly: - -```julia -# Create transform (coefficients computed once) -ft = FourierTransform(mtheta, mpert, mlow) - -# Forward transform (callable): theta → modes -modes = ft(theta_data) - -# Inverse transform: modes → theta -theta_reconstructed = inverse(ft, modes) -``` - -See also: [`compute_fourier_coefficients`](@ref), [`inverse`](@ref) + - `mtheta`: poloidal grid size + - `mpert`: number of poloidal modes + - `mlow`: lowest poloidal mode number + - `basis`: ``\\exp(-i(m\\theta - n\\nu))``, size `(mpert, mtheta)` """ struct FourierTransform mtheta::Int mpert::Int mlow::Int - cslth::Matrix{Float64} - snlth::Matrix{Float64} + basis::Matrix{ComplexF64} end """ - FourierTransform(mtheta, mpert, mlow; n=0, qa=0.0, delta=zeros(mtheta)) + FourierTransform(mtheta, mpert, mlow; n=0, ν=zeros(mtheta)) -Construct a FourierTransform object with pre-computed basis functions. +Construct a transform with precomputed basis for contiguous modes `mlow:(mlow+mpert-1)`. -# Arguments +## Arguments -- `mtheta::Int`: Number of poloidal grid points -- `mpert::Int`: Number of Fourier modes -- `mlow::Int`: Lowest mode number + - `mtheta`: number of poloidal grid points + - `mpert`: number of poloidal modes + - `mlow`: lowest poloidal mode number -# Keyword Arguments +## Keyword Arguments -- `n::Int=0`: Toroidal mode number -- `qa::Float64=0.0`: Safety factor at boundary -- `delta::Vector{Float64}=zeros(mtheta)`: Toroidal phase shift array + - `n`: toroidal mode number (default 0) + - `ν`: toroidal angle offset on the poloidal grid, length `mtheta` -# Returns +## Returns -`FourierTransform` object ready for forward and inverse transforms. - -# Examples - -```julia -# Simple case: no phase shift (for PerturbedEquilibrium) -ft = FourierTransform(480, 40, -20) - -# With toroidal phase (for Vacuum calculations) -ft_vac = FourierTransform(480, 40, -20; n=1, ν=ν_array) -``` + - `FourierTransform` ready for forward and inverse transforms """ function FourierTransform( mtheta::Int, @@ -213,569 +133,55 @@ function FourierTransform( n::Int=0, ν::Vector{Float64}=zeros(Float64, mtheta) ) - cos_mn_basis, sin_mn_basis = compute_fourier_coefficients(mtheta, mpert, mlow, 1, 1, 1; n_2D=n, ν=ν) - return FourierTransform(mtheta, mpert, mlow, cos_mn_basis, sin_mn_basis) + basis = compute_fourier_coefficients(mtheta, mlow:(mlow+mpert-1), n, ν) + return FourierTransform(mtheta, mpert, mlow, basis) end """ - (ft::FourierTransform)(data::AbstractVecOrMat{<:Real}) - -Forward Fourier transform: theta-space → mode-space (callable functor). - -Transforms real-valued data at theta grid points into complex Fourier mode coefficients. - -# Arguments - -- `data::AbstractVecOrMat{Float64}`: Data at theta points - - If `Vector{Float64}` with length `mtheta`: single function to transform - - If `Matrix{Float64}` with size `(mtheta, n)`: n functions to transform simultaneously + (ft::FourierTransform)(data) -# Returns +Forward transform from θ-space to mode space. -- Complex Fourier coefficients - - If input is `Vector`: returns `Vector{ComplexF64}` of length `mpert` - - If input is `Matrix`: returns `Matrix{ComplexF64}` of size `(mpert, n)` +## Arguments -# Formula + - `data`: real or complex samples on the poloidal grid — `Vector{mtheta}` or `Matrix{mtheta, :}` -For each mode `l` (corresponding to mode number `m = mlow + l - 1`), using the -`exp(-imθ)` convention (Fortran `iscdftf`): +## Returns -``` -mode[l] = (1/N) Σᵢ data[i] * exp(-i*(m*θᵢ + phase)) - = (1/N) Σᵢ data[i] * (cos(m*θᵢ + phase) - im*sin(m*θᵢ + phase)) -``` - -Identical convention to the complex-input method, so `inverse(ft, ft(data)) ≈ data` -holds for both real- and complex-valued `data`. - -# Examples - -```julia -ft = FourierTransform(480, 40, -20) - -# Transform a single function -f_theta = sin.(theta_grid) -f_modes = ft(f_theta) # Vector{ComplexF64} of length 40 - -# Transform multiple functions at once -data = randn(480, 10) # 10 different functions -modes = ft(data) # Matrix{ComplexF64} of size (40, 10) -``` + - Mode coefficients — `Vector{mpert}` or `Matrix{mpert, :}` """ -function (ft::FourierTransform)(data::AbstractVecOrMat{<:Real}) - # Forward transform with 1/N normalization and exp(-imθ) convention (Fortran iscdftf), - # identical to the complex-input method below. - +function (ft::FourierTransform)(data::AbstractVecOrMat{<:Number}) if data isa AbstractVector - # For vector input: [mtheta] → [mpert] @assert length(data) == ft.mtheta "Input vector must have length mtheta=$(ft.mtheta)" - real_part = ft.cslth' * data - imag_part = ft.snlth' * data - return complex.(real_part, -imag_part) ./ ft.mtheta else - # For matrix input: [mtheta, n] → [mpert, n] @assert size(data, 1) == ft.mtheta "Input matrix first dimension must be mtheta=$(ft.mtheta)" - real_part = ft.cslth' * data - imag_part = ft.snlth' * data - return complex.(real_part, -imag_part) ./ ft.mtheta end -end - -""" - (ft::FourierTransform)(data::AbstractVecOrMat{<:Complex}) - -Forward Fourier transform for complex-valued theta-space data. - -Transforms complex data at theta grid points into complex Fourier mode coefficients. - -# Arguments - -- `data::AbstractVecOrMat{ComplexF64}`: Complex data at theta points - -# Returns - -- Complex Fourier coefficients with same shape as real input case - -# Formula - -For complex input `data = data_real + im*data_imag`, computes using exp(-imθ): - -``` -Re{mode[l]} = Σᵢ (data_real[i]*cos + data_imag[i]*sin) -Im{mode[l]} = Σᵢ (-data_real[i]*sin + data_imag[i]*cos) -``` -where cos and sin are the basis functions at theta point i for mode l. -This matches the Fortran `iscdftf` convention: `f_m = (1/N) Σ_j f(θ_j) exp(-i m θ_j)`. -""" -function (ft::FourierTransform)(data::AbstractVecOrMat{<:Complex}) - # Forward transform with 1/N normalization and exp(-imθ) convention - # (matches Fortran iscdftf: f_m = (1/N) Σ_j f_j * exp(-2πi*m*j/N)) - - if data isa AbstractVector - @assert length(data) == ft.mtheta "Input vector must have length mtheta=$(ft.mtheta)" - real_part = ft.cslth' * real.(data) .+ ft.snlth' * imag.(data) - imag_part = -(ft.snlth' * real.(data)) .+ ft.cslth' * imag.(data) - return complex.(real_part, imag_part) ./ ft.mtheta - else - @assert size(data, 1) == ft.mtheta "Input matrix first dimension must be mtheta=$(ft.mtheta)" - real_part = ft.cslth' * real.(data) .+ ft.snlth' * imag.(data) - imag_part = -(ft.snlth' * real.(data)) .+ ft.cslth' * imag.(data) - return complex.(real_part, imag_part) ./ ft.mtheta - end + return (ft.basis * data) ./ ft.mtheta end """ - inverse(ft::FourierTransform, modes::AbstractVecOrMat{<:Complex}) - -Inverse Fourier transform: mode-space → theta-space. - -Reconstructs theta-space data from complex Fourier mode coefficients. - -# Arguments - -- `modes::AbstractVecOrMat{ComplexF64}`: Fourier mode coefficients - - If `Vector{ComplexF64}` with length `mpert`: single mode expansion to reconstruct - - If `Matrix{ComplexF64}` with size `(mpert, n)`: n mode expansions to reconstruct - -# Returns - -- Reconstructed data at theta points - - If input is `Vector`: returns `Vector{ComplexF64}` of length `mtheta` - - If input is `Matrix`: returns `Matrix{ComplexF64}` of size `(mtheta, n)` + inverse(ft, modes) -# Formula +Inverse transform from mode space to θ-space. -For each theta point `i`: +## Arguments -``` -data[i] = (2π/mtheta) * Σₗ modes[l] * (cos(m*θᵢ + phase) + im*sin(m*θᵢ + phase)) -``` + - `ft`: precomputed transform + - `modes`: complex mode coefficients — `Vector{mpert}` or `Matrix{mpert, :}` -This is equivalent to: -``` -data[i] = (2π/mtheta) * Σₗ [Re{modes[l]}*cos - Im{modes[l]}*sin + - im*(Re{modes[l]}*sin + Im{modes[l]}*cos)] -``` +## Returns -# Notes - -The normalization factor `2π/mtheta` ensures proper Fourier series reconstruction. - -For real-valued theta data, the result will have negligible imaginary parts (machine precision). -Use `real.(inverse(ft, modes))` to extract the real part if needed. - -# Examples - -```julia -ft = FourierTransform(480, 40, -20) - -# Reconstruct from modes -modes = randn(ComplexF64, 40) -theta_data = inverse(ft, modes) # Vector{ComplexF64} of length 480 - -# For real reconstruction -theta_real = real.(inverse(ft, modes)) -``` + - θ-space data — `Vector{mtheta}` or `Matrix{mtheta, :}` """ function inverse(ft::FourierTransform, modes::AbstractVecOrMat{<:Complex}) - # Inverse transform without normalization (matches Fortran iscdftb convention) - # f(θ) = Σₗ cₗ * exp(i*m*θ) - # Round-trip: forward(inverse(x)) = x (1/N in forward cancels N summation terms) - - if modes isa AbstractVector - @assert length(modes) == ft.mpert "Input vector must have length mpert=$(ft.mpert)" - real_part = ft.cslth * real.(modes) .- ft.snlth * imag.(modes) - imag_part = ft.cslth * imag.(modes) .+ ft.snlth * real.(modes) - return complex.(real_part, imag_part) - else - @assert size(modes, 1) == ft.mpert "Input matrix first dimension must be mpert=$(ft.mpert)" - real_part = ft.cslth * real.(modes) .- ft.snlth * imag.(modes) - imag_part = ft.cslth * imag.(modes) .+ ft.snlth * real.(modes) - return complex.(real_part, imag_part) - end -end - -""" - inverse(ft::FourierTransform, modes::AbstractVecOrMat{<:Real}) - -Inverse Fourier transform for real mode coefficients (pure cosine expansion). - -Special case where all modes are real-valued, corresponding to a cosine-only Fourier series. - -# Arguments - -- `modes::AbstractVecOrMat{Float64}`: Real Fourier coefficients - -# Returns - -- Real-valued data at theta points with normalization factor `2π/mtheta` - -# Formula - -``` -data[i] = (2π/mtheta) * Σₗ modes[l] * cos(m*θᵢ + phase) -``` - -# Notes - -This is a special case and less commonly used. Most applications use complex modes. -""" -function inverse(ft::FourierTransform, modes::AbstractVecOrMat{<:Real}) - dth = 2π / ft.mtheta - - if modes isa AbstractVector - @assert length(modes) == ft.mpert "Input vector must have length mpert=$(ft.mpert)" - return (ft.cslth * modes) .* dth - else - @assert size(modes, 1) == ft.mpert "Input matrix first dimension must be mpert=$(ft.mpert)" - return (ft.cslth * modes) .* dth - end -end - -# ============================================================================== -# In-place transform functions for performance-critical applications -# ============================================================================== - -""" - transform!(output::AbstractVecOrMat{ComplexF64}, ft::FourierTransform, data::AbstractVecOrMat{<:Real}) - -In-place forward Fourier transform for real-valued theta-space data. - -More efficient than the allocating version `ft(data)` when called repeatedly, -as it reuses the pre-allocated output array. - -# Arguments - -- `output::AbstractVecOrMat{ComplexF64}`: Pre-allocated output array for Fourier modes - - For vector input: must have length `mpert` - - For matrix input: must have size `(mpert, n)` where `n = size(data, 2)` -- `ft::FourierTransform`: Fourier transform object with pre-computed coefficients -- `data::AbstractVecOrMat{Float64}`: Real-valued data at theta points - - For vector: length must be `mtheta` - - For matrix: first dimension must be `mtheta` - -# Returns - -- `output`: The modified output array (for chaining) - -# Performance - -This function uses in-place matrix multiplication (`mul!`) to avoid allocations, -making it suitable for tight loops and performance-critical code. - -# Example - -```julia -ft = FourierTransform(480, 40, -20) - -# Pre-allocate output buffer -modes = zeros(ComplexF64, 40) - -# Reuse buffer in loop -for i in 1:1000 - data = get_theta_data(i) - transform!(modes, ft, data) # No allocations - process_modes(modes) -end -``` -""" -function transform!(output::AbstractVecOrMat{ComplexF64}, ft::FourierTransform, data::AbstractVecOrMat{<:Real}) - # exp(-imθ) convention (Fortran iscdftf), matching the allocating `ft(data::Real)`. - if data isa AbstractVector - @assert length(data) == ft.mtheta "Input vector must have length mtheta=$(ft.mtheta)" - @assert length(output) == ft.mpert "Output vector must have length mpert=$(ft.mpert)" - - # Extract real and imaginary views - real_part = reinterpret(Float64, output) - real_view = @view real_part[1:2:end] # Real components - imag_view = @view real_part[2:2:end] # Imaginary components - - # In-place computation: real = cslth' * data, imag = -snlth' * data - mul!(real_view, ft.cslth', data) - mul!(imag_view, ft.snlth', data) - imag_view .*= -1 - output ./= ft.mtheta - else - @assert size(data, 1) == ft.mtheta "Input matrix first dimension must be mtheta=$(ft.mtheta)" - @assert size(output, 1) == ft.mpert "Output matrix first dimension must be mpert=$(ft.mpert)" - @assert size(output, 2) == size(data, 2) "Output and input must have same number of columns" - - # For matrices, we need temporary storage for real/imag parts - n_cols = size(data, 2) - real_part = similar(output, Float64, ft.mpert, n_cols) - imag_part = similar(output, Float64, ft.mpert, n_cols) - - # In-place computation - mul!(real_part, ft.cslth', data) - mul!(imag_part, ft.snlth', data) - - # Combine into complex output with 1/N normalization - output .= complex.(real_part, -imag_part) ./ ft.mtheta - end - - return output -end - -""" - transform!(output::AbstractVecOrMat{ComplexF64}, ft::FourierTransform, data::AbstractVecOrMat{<:Complex}) - -In-place forward Fourier transform for complex-valued theta-space data. - -# Arguments - -- `output::AbstractVecOrMat{ComplexF64}`: Pre-allocated output array for Fourier modes -- `ft::FourierTransform`: Fourier transform object -- `data::AbstractVecOrMat{ComplexF64}`: Complex-valued data at theta points - -# Returns - -- `output`: The modified output array - -# Formula - -For complex input `data = data_real + im*data_imag`, uses exp(-imθ): -``` -Re{mode[l]} = Σᵢ (data_real[i]*cos + data_imag[i]*sin) -Im{mode[l]} = Σᵢ (-data_real[i]*sin + data_imag[i]*cos) -``` -""" -function transform!(output::AbstractVecOrMat{ComplexF64}, ft::FourierTransform, data::AbstractVecOrMat{<:Complex}) - if data isa AbstractVector - @assert length(data) == ft.mtheta "Input vector must have length mtheta=$(ft.mtheta)" - @assert length(output) == ft.mpert "Output vector must have length mpert=$(ft.mpert)" - - # Temporary storage for intermediate results - real_part = similar(output, Float64) - imag_part = similar(output, Float64) - temp1 = similar(output, Float64) - temp2 = similar(output, Float64) - - # Re{mode} = cslth' * real(data) + snlth' * imag(data) - mul!(temp1, ft.cslth', real.(data)) - mul!(temp2, ft.snlth', imag.(data)) - real_part .= temp1 .+ temp2 - - # Im{mode} = -snlth' * real(data) + cslth' * imag(data) - mul!(temp1, ft.cslth', imag.(data)) - mul!(temp2, ft.snlth', real.(data)) - imag_part .= temp1 .- temp2 - - output .= complex.(real_part, imag_part) ./ ft.mtheta - else - @assert size(data, 1) == ft.mtheta "Input matrix first dimension must be mtheta=$(ft.mtheta)" - @assert size(output, 1) == ft.mpert "Output matrix first dimension must be mpert=$(ft.mpert)" - @assert size(output, 2) == size(data, 2) "Output and input must have same number of columns" - - n_cols = size(data, 2) - real_part = similar(output, Float64, ft.mpert, n_cols) - imag_part = similar(output, Float64, ft.mpert, n_cols) - temp1 = similar(real_part) - temp2 = similar(real_part) - - # Re{mode} = cslth' * real(data) + snlth' * imag(data) - mul!(temp1, ft.cslth', real.(data)) - mul!(temp2, ft.snlth', imag.(data)) - real_part .= temp1 .+ temp2 - - # Im{mode} = -snlth' * real(data) + cslth' * imag(data) - mul!(temp1, ft.cslth', imag.(data)) - mul!(temp2, ft.snlth', real.(data)) - imag_part .= temp1 .- temp2 - - output .= complex.(real_part, imag_part) ./ ft.mtheta - end - - return output -end - -""" - inverse_transform!(output::AbstractVecOrMat{ComplexF64}, ft::FourierTransform, modes::AbstractVecOrMat{<:Complex}) - -In-place inverse Fourier transform: mode-space → theta-space. - -More efficient than the allocating version `inverse(ft, modes)` when called repeatedly. - -# Arguments - -- `output::AbstractVecOrMat{ComplexF64}`: Pre-allocated output array for theta-space data - - For vector: must have length `mtheta` - - For matrix: must have size `(mtheta, n)` where `n = size(modes, 2)` -- `ft::FourierTransform`: Fourier transform object -- `modes::AbstractVecOrMat{ComplexF64}`: Complex Fourier mode coefficients - - For vector: length must be `mpert` - - For matrix: first dimension must be `mpert` - -# Returns - -- `output`: The modified output array (for chaining) - -# Formula - -``` -data[i] = (2π/mtheta) * Σₗ [Re{modes[l]}*cos - Im{modes[l]}*sin + - im*(Re{modes[l]}*sin + Im{modes[l]}*cos)] -``` - -# Example - -```julia -ft = FourierTransform(480, 40, -20) - -# Pre-allocate output buffer -theta_data = zeros(ComplexF64, 480) - -# Reuse buffer in loop -for i in 1:1000 - modes = get_modes(i) - inverse_transform!(theta_data, ft, modes) # No allocations - process_theta_data(theta_data) -end -``` -""" -function inverse_transform!(output::AbstractVecOrMat{ComplexF64}, ft::FourierTransform, modes::AbstractVecOrMat{<:Complex}) - # Inverse transform without normalization (matches Fortran iscdftb convention) - if modes isa AbstractVector @assert length(modes) == ft.mpert "Input vector must have length mpert=$(ft.mpert)" - @assert length(output) == ft.mtheta "Output vector must have length mtheta=$(ft.mtheta)" - - # Temporary storage - real_part = similar(output, Float64) - imag_part = similar(output, Float64) - temp1 = similar(output, Float64) - temp2 = similar(output, Float64) - - # Re{data} = cslth * real(modes) - snlth * imag(modes) - mul!(temp1, ft.cslth, real.(modes)) - mul!(temp2, ft.snlth, imag.(modes)) - real_part .= temp1 .- temp2 - - # Im{data} = cslth * imag(modes) + snlth * real(modes) - mul!(temp1, ft.cslth, imag.(modes)) - mul!(temp2, ft.snlth, real.(modes)) - imag_part .= temp1 .+ temp2 - - output .= complex.(real_part, imag_part) else @assert size(modes, 1) == ft.mpert "Input matrix first dimension must be mpert=$(ft.mpert)" - @assert size(output, 1) == ft.mtheta "Output matrix first dimension must be mtheta=$(ft.mtheta)" - @assert size(output, 2) == size(modes, 2) "Output and input must have same number of columns" - - n_cols = size(modes, 2) - real_part = similar(output, Float64, ft.mtheta, n_cols) - imag_part = similar(output, Float64, ft.mtheta, n_cols) - temp1 = similar(real_part) - temp2 = similar(real_part) - - # Re{data} = cslth * real(modes) - snlth * imag(modes) - mul!(temp1, ft.cslth, real.(modes)) - mul!(temp2, ft.snlth, imag.(modes)) - real_part .= temp1 .- temp2 - - # Im{data} = cslth * imag(modes) + snlth * real(modes) - mul!(temp1, ft.cslth, imag.(modes)) - mul!(temp2, ft.snlth, real.(modes)) - imag_part .= temp1 .+ temp2 - - output .= complex.(real_part, imag_part) end - return output -end - -# ============================================================================== -# Low-level matrix transforms with offsets (for Vacuum module compatibility) -# ============================================================================== - -""" - fourier_transform!(gil::AbstractMatrix{Float64}, gij::AbstractMatrix{Float64}, cs::Matrix{Float64}; row_offset::Int=0, col_offset::Int=0) - -Low-level Fourier transform with offset support for Vacuum module. - -Performs a truncated Fourier transform of `gij` onto `gil` using pre-computed -Fourier coefficients `cs`, with support for block offsets in the output matrix. - -This is used by the Vacuum module for transforming Green's function matrices that -have specific block structure (e.g., plasma and wall contributions packed together). - -# Arguments - -- `gil::AbstractMatrix{Float64}`: Output matrix, updated in-place at offset block -- `gij::AbstractMatrix{Float64}`: Input matrix (mtheta × mtheta) containing theta-space data -- `cs::Matrix{Float64}`: Fourier coefficient matrix (mtheta × mpert) -- `row_offset::Int`: Row offset in `gil` matrix -- `col_offset::Int`: Column offset in `gil` matrix - -# Operation - -Computes: `gil[row_offset+i, col_offset+l] = Σⱼ gij[i, j] * cs[j, l]` for i ∈ 1:size(cs, 1), l ∈ 1:size(cs, 2) - -# Notes - -- This function uses 0-based offset convention (add 1 for Julia indexing) -- Used for packing real and imaginary parts in separate blocks of `gil` - -# Example - -```julia -ft = FourierTransform(mtheta, mpert, mlow; n=n, qa=qa, delta=delta) -gil = zeros(2*mtheta, 2*mpert) # Packed real/imag structure -gij = ... # Some theta-space data - -# Transform using cosine coefficients, real part block (offset 0, 0) -fourier_transform!(gil, gij, ft.cslth; row_offset=0, col_offset=0) - -# Transform using sine coefficients, imaginary part block (offset 0, mpert) -fourier_transform!(gil, gij, ft.snlth; row_offset=0, col_offset=mpert) -``` -""" -function fourier_transform!(gil::AbstractMatrix{Float64}, gij::AbstractMatrix{Float64}, cs::Matrix{Float64}; row_offset::Int=0, col_offset::Int=0) - mul!(view(gil, row_offset+1:row_offset+size(cs, 1), col_offset+1:col_offset+size(cs, 2)), gij, cs) -end - -""" - fourier_inverse_transform!(gll::Matrix{Float64}, gil::Matrix{Float64}, cs::Matrix{Float64}, m00::Int, l00::Int) - -Low-level inverse Fourier transform with offset support for Vacuum module. - -Performs the inverse Fourier transform of `gil` onto `gll` using pre-computed -Fourier coefficients `cs`, with support for block offsets in the input matrix. - -This is used by the Vacuum module for inverse transforming Green's function matrices -back to mode space from theta space. - -# Arguments - -- `gll::Matrix{Float64}`: Output matrix (mpert × mpert), updated in-place -- `gil::Matrix{Float64}`: Input matrix containing Fourier-transformed data -- `cs::Matrix{Float64}`: Fourier coefficient matrix (mtheta × mpert), either `cslth` or `snlth` -- `m00::Int`: Row offset in `gil` matrix (0-based indexing convention) -- `l00::Int`: Column offset in `gil` matrix (0-based indexing convention) - -# Operation - -Computes: `gll[l2, l1] = dθdζ * Σᵢ cs[i, l2] * gil[m00+i, l00+l1]` - -# Notes - -- The normalization factor `1 / size(cs, 1)` matches the 1/N forward convention -- This function uses 0-based offset convention (add 1 for Julia indexing) -- The `cs` matrix should be either `cslth` or `snlth` from a `FourierTransform` object -# Example - -```julia -ft = FourierTransform(mtheta, mpert, mlow; n=n, qa=qa, delta=delta) -gil = ... # Transformed Green's function data -arr = zeros(mpert, mpert) - -# Inverse transform from real part block using cosine coefficients -fourier_inverse_transform!(arr, gil, ft.cslth) -``` -""" -function fourier_inverse_transform!(gll::AbstractMatrix{Float64}, gil::AbstractMatrix{Float64}, cs::Matrix{Float64}; row_offset::Int=0, col_offset::Int=0) - mul!(gll, cs', view(gil, row_offset+1:row_offset+size(cs, 1), col_offset+1:col_offset+size(cs, 2)), 1.0 / size(cs, 1), 0.0) + return adjoint(ft.basis) * modes end end # module FourierTransforms diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl index 7cd036f1..060a4857 100644 --- a/src/Utilities/Utilities.jl +++ b/src/Utilities/Utilities.jl @@ -19,8 +19,6 @@ include("GridUtilities.jl") using .FourierTransforms export FourierTransform, inverse, compute_fourier_coefficients -export transform!, inverse_transform! -export fourier_transform!, fourier_inverse_transform! export FourierCoefficients, empty_FourierCoefficients, get_complex_coeff, get_complex_coeffs! diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index 9a406199..4ba746ba 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -5,7 +5,7 @@ Struct holding plasma boundary and mode data as provided from ForceFreeStates na 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. +nzeta > 1 for 3D vacuum calculation. The arrays should be for a single field period only, with excluded endpoints. # Fields @@ -15,13 +15,12 @@ nzeta > 1 for 3D vacuum calculation. - `ν::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 + - `m_modes::Vector{Int}`: Vector of poloidal mode numbers. E.g. `collect(mlow:mhigh)` for a contiguous range. + - `n_modes::Vector{Int}`: Vector of toroidal mode numbers. E.g. `collect(nlow:nhigh)` for a contiguous range or `collect(nlow:n_stride:nhigh)` for a strided list for stellarator mode-family calculations. - `mtheta::Int`: Number of vacuum calculation poloidal 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 + - `nfp::Int`: Number of field periods """ @kwdef struct VacuumInput x::Vector{Float64} = Float64[] @@ -30,13 +29,12 @@ nzeta > 1 for 3D vacuum calculation. ν::Vector{Float64} = Float64[] mtheta_in::Int = 0 nzeta_in::Int = 1 - mlow::Int = 1 - mpert::Int = 1 - nlow::Int = 1 - npert::Int = 1 + m_modes::Vector{Int} = [1] + n_modes::Vector{Int} = [1] mtheta::Int = 1 nzeta::Int = 1 force_wv_symmetry::Bool = true + nfp::Int = 1 end """ @@ -44,9 +42,9 @@ end equil::Equilibrium.PlasmaEquilibrium, ψ::Float64, mtheta::Int, - mpert::Int, - mlow::Int, - n::Int, + nzeta::Int, + m_modes::AbstractVector{<:Integer}, + n_modes::AbstractVector{<:Integer}; force_wv_symmetry::Bool = true ) -> VacuumInput @@ -58,10 +56,10 @@ Extracts plasma geometry from equilibrium at the given flux surface and packages - `equil`: Equilibrium solution - `ψ`: Normalized flux coordinate - `mtheta`: Number of vacuum calculation poloidal points - - `mpert`: Number of perturbing poloidal modes - - `mlow`: Lowest poloidal mode number - - `n`: Toroidal mode number - - `force_wv_symmetry::Bool`: Boolean flag to enforce symmetry in the vacuum response matrix (default: true) + - `nzeta`: Number of vacuum calculation toroidal points (1 for 2D, >1 for 3D) + - `m_modes`: Poloidal mode numbers (e.g. `mlow:mhigh`) + - `n_modes`: Toroidal mode numbers (e.g. `[n]` for a single mode, or `nlow:nhigh`) + - `force_wv_symmetry::Bool`: Enforce Hermitian symmetry in the vacuum response matrix (default: true) ## Returns @@ -72,30 +70,86 @@ function VacuumInput( ψ::Float64, mtheta::Int, nzeta::Int, - mpert::Int, - mlow::Int, - npert::Int, - nlow::Int; + m_modes::AbstractVector{<:Integer}, + n_modes::AbstractVector{<:Integer}; force_wv_symmetry::Bool=true ) # 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 return VacuumInput(; - x=reverse(r), - z=reverse(z), - ν=reverse(ν), - mtheta_in=length(r), - mlow=mlow, - mpert=mpert, - nlow=nlow, - npert=npert, + x=reverse(r)[1:(end-1)], + z=reverse(z)[1:(end-1)], + ν=reverse(ν)[1:(end-1)], + mtheta_in=length(r)-1, + m_modes=collect(Int, m_modes), + n_modes=collect(Int, n_modes), mtheta=mtheta, nzeta=nzeta, force_wv_symmetry=force_wv_symmetry ) end +""" + expand_field_periods(inputs::VacuumInput) -> VacuumInput + +Reconstruct the full-torus boundary from a single field period using the device field +periodicity. When `inputs.nfp > 1` (and the boundary is 3D), the input `x, y, z` arrays +describe one field period spanning ζ ∈ [0, 2π/nfp) on an `mtheta_in × nzeta_in` grid; the +full torus is generated by `nfp` rigid rotations about the Z-axis by 2π·k/nfp. The toroidal +counts `nzeta_in` and `nzeta` are scaled to their full-torus totals and `nfp` is reset to 1, +so all downstream code operates on the full surface unchanged. + +This is the layer-1 interface: it shrinks the Python→Julia payload to one field period while +producing a full-torus geometry identical (to floating-point) to sending the whole torus. + +Returns `inputs` unchanged for axisymmetric (`nzeta_in == 1`) or single-period (`nfp <= 1`) cases. +""" +function expand_field_periods(inputs::VacuumInput) + (inputs.nfp <= 1 || inputs.nzeta_in <= 1) && return inputs + + nfp = inputs.nfp + mt = inputs.mtheta_in + nz = inputs.nzeta_in # points per field period + npts = mt * nz + @assert length(inputs.x) == length(inputs.y) == length(inputs.z) == npts "per-period x, y, z must have length mtheta_in * nzeta_in" + + nz_full = nz * nfp + xf = Vector{Float64}(undef, mt * nz_full) + yf = similar(xf) + zf = similar(xf) + + # Tile to full torus: period k is the single period rotated by α = 2π·k/nfp about Z. + for k in 0:(nfp-1) + α = 2π * k / nfp + ca, sa = cos(α), sin(α) + @inbounds for j in 1:nz, i in 1:mt + src = i + mt * (j - 1) + dst = i + mt * (j - 1 + k * nz) + xf[dst] = inputs.x[src] * ca - inputs.y[src] * sa + yf[dst] = inputs.x[src] * sa + inputs.y[src] * ca + zf[dst] = inputs.z[src] + end + end + + return VacuumInput(; + x=xf, + y=yf, + z=zf, + ν=inputs.ν, + mtheta_in=mt, + nzeta_in=nz_full, + m_modes=inputs.m_modes, + n_modes=inputs.n_modes, + mtheta=inputs.mtheta, + nzeta=inputs.nzeta * nfp, + force_wv_symmetry=inputs.force_wv_symmetry, + nfp=1 + ) +end + """ WallShapeSettings @@ -146,26 +200,6 @@ Struct containing input settings for vacuum wall geometry. equal_arc_wall::Bool = true end -""" - KernelParams2D(n::Int) - -Parameter struct for 2D vacuum kernel dispatch. Holds the toroidal mode number `n`. -""" -struct KernelParams2D - n::Int -end - -""" - KernelParams3D(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) - -Parameter struct for 3D vacuum kernel dispatch. Holds singular quadrature parameters. -""" -struct KernelParams3D - PATCH_RAD::Int - RAD_DIM::Int - INTERP_ORDER::Int -end - """ PlasmaGeometry @@ -222,13 +256,13 @@ function PlasmaGeometry(inputs::VacuumInput) @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=inputs.mtheta_in) # VacuumInput uses [0, 2π] grid - θ_out = range(; start=0, length=inputs.mtheta, step=2π/inputs.mtheta) # VACUUM uses [0, 2π) grid + θ_in = range(; start=0, length=inputs.mtheta_in, step=2π/inputs.mtheta_in) + θ_out = range(; start=0, length=inputs.mtheta, step=2π/inputs.mtheta) # Use one-shot API with PeriodicBC - 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()) + x = cubic_interp(θ_in, inputs.x, θ_out; bc=PeriodicBC(; endpoint=:exclusive)) + z = cubic_interp(θ_in, inputs.z, θ_out; bc=PeriodicBC(; endpoint=:exclusive)) + ν = cubic_interp(θ_in, inputs.ν, θ_out; bc=PeriodicBC(; endpoint=:exclusive)) return PlasmaGeometry(x, z, ν) end @@ -339,13 +373,13 @@ function PlasmaGeometry3D(inputs::VacuumInput) # 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) + θ_in = range(; start=0, length=inputs.mtheta_in, step=2π/inputs.mtheta_in) + ζ_in = range(; start=0, length=inputs.nzeta_in, step=2π/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())) + 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 end @@ -428,10 +462,6 @@ function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_set x_plasma = plasma_surf.x z_plasma = plasma_surf.z - # Output wall coordinate arrays - x_wall = zeros(Float64, mtheta) - z_wall = zeros(Float64, mtheta) - # Common geometric parameters xmin = minimum(x_plasma) xmax = maximum(x_plasma) @@ -615,6 +645,7 @@ function WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings) end inputs.nzeta_in > 1 && error("3D wall geometry not yet implemented for non-axisymmetric inputs") + inputs.nfp > 1 && error("3D wall geometry not yet implemented for nfp > 1") # Plasma surface coordinates (2D) surf_2D = PlasmaGeometry(inputs) diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index eac91b01..ea93f224 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/src/Vacuum/Kernel2D.jl @@ -17,7 +17,7 @@ end ) end -# Precomputed Gauss-Legendre rule used in the hot path (`kernel!`). +# Precomputed Gauss-Legendre rule const GL8 = gausslegendre_rule(Val(8)) """ @@ -59,7 +59,7 @@ const GL8_LAGRANGE_STENCILS = precompute_lagrange_stencils(GL8.x) # and per-n sinh/cosh cache are defined in PnQuadCache.jl. """ - kernel!(grad_greenfunction, greenfunction, observer, source, n) + compute_2D_kernel_matrices!(grad_greenfunction, greenfunction, observer, source, n) Compute kernels of integral equation for Laplace's equation in a torus. **WARNING: This kernel only supports closed toroidal walls currently. @@ -76,8 +76,8 @@ The residue calculation needs to be updated for open walls.** # Returns Modifies `grad_greenfunction` and `greenfunction` in place. -Note that greenfunction is zeroed each time this function is called, -but grad_greenfunction is not since it fills a different block of the +Note that greenfunction is zeroed only when the source is plasma; +grad_greenfunction is not zeroed since it fills a different block of the (2 * mtheta, 2 * mtheta) depending on the source/observer. # Notes @@ -107,10 +107,9 @@ but grad_greenfunction is not since it fills a different block of the ((col_index-1)*mtheta+1):(col_index*mtheta) ) - # Zero out greenfunction at start of each kernel call - fill!(greenfunction, 0.0) # 𝒢ⁿ only needed for plasma as source term (RHS of eqs. 26/27 in Chance 1997) populate_greenfunction = source isa PlasmaGeometry + populate_greenfunction && fill!(greenfunction, 0.0) # S₁ᵢ logarithmic correction factors [Chance Phys. Plasmas 1997 2161 eq. 78] log_correction_0=16.0*dtheta*(log(2*dtheta)-68.0/15.0)/15.0 @@ -233,17 +232,6 @@ but grad_greenfunction is not since it fills a different block of the end end -# Dispatch wrapper for unified 2D/3D vacuum: forwards to 5-arg compute_2D_kernel_matrices! with params.n -function kernel!( - grad_greenfunction::AbstractMatrix{Float64}, - greenfunction::AbstractMatrix{Float64}, - observer::Union{PlasmaGeometry,WallGeometry}, - source::Union{PlasmaGeometry,WallGeometry}, - params::KernelParams2D -) - return compute_2D_kernel_matrices!(grad_greenfunction, greenfunction, observer, source, params.n) -end - ############################################################# # Legendre function of the first kind eq.(47)~(50) , replacing aleg. (verified) ############################################################# diff --git a/src/Vacuum/Kernel3D.jl b/src/Vacuum/Kernel3D.jl index 55a4eb95..3d135295 100644 --- a/src/Vacuum/Kernel3D.jl +++ b/src/Vacuum/Kernel3D.jl @@ -334,30 +334,17 @@ grad_greenfunction is the double-layer kernel matrix, where each entry is ∇_{x_src} φ(x_obs, x_src) · n_src, and greenfunction is the single-layer kernel matrix, where each entry is φ(x_obs, x_src). +Takes advantage of field periodicity to evaluate the kernel only over a single field period. + # Arguments - `grad_greenfunction`: Double-layer kernel matrix (Nobs × Nsrc) filled in place - - - `greenfunction`: Single-layer kernel matrix (Nobs × Nsrc) filled in place - + - `greenfunction`: Single-layer kernel matrix (Nobs × Nsrc); filled only when `source` is plasma - `observer`: Observer geometry (PlasmaGeometry3D) - - `source`: Source geometry (PlasmaGeometry3D) - - `PATCH_RAD`: Number of points adjacent to source point to treat as singular - - + Total patch size in # of gridpoints = (2 * PATCH_RAD + 1) x (2 * PATCH_RAD + 1) - - `RAD_DIM`: Polar radial quadrature order. Angular order = 2 * RAD_DIM - - - `INTERP_ORDER`: Lagrange interpolation order - - + Must be ≤ (2 * PATCH_RAD + 1) - -# Threading - -This function automatically uses all available threads (`Threads.nthreads()`). -Start Julia with `julia -t auto` or set `JULIA_NUM_THREADS` to enable multi-threading. + - `INTERP_ORDER`: Lagrange interpolation order, must be ≤ (2 * PATCH_RAD + 1) """ function compute_3D_kernel_matrices!( grad_greenfunction::AbstractMatrix{Float64}, @@ -369,6 +356,7 @@ function compute_3D_kernel_matrices!( INTERP_ORDER::Int ) num_points = observer.mtheta * observer.nzeta + n_obs = size(greenfunction, 1) # num_points ÷ nfp dθdζ = 4π^2 / num_points # Get block of grad green function matrix @@ -376,14 +364,13 @@ function compute_3D_kernel_matrices!( row_index = (observer isa PlasmaGeometry3D ? 1 : 2) grad_greenfunction_block = view( grad_greenfunction, - ((row_index-1)*num_points+1):(row_index*num_points), + ((row_index-1)*n_obs+1):(row_index*n_obs), ((col_index-1)*num_points+1):(col_index*num_points) ) - # Zero out green function matrix - fill!(greenfunction, 0.0) # 𝒢ⁿ only needed for plasma as source term (RHS of eqs. 26/27 in Chance 1997) populate_greenfunction = source isa PlasmaGeometry3D + populate_greenfunction && fill!(greenfunction, 0.0) # This allows the code to run at lower resolution without erroring out, but will warn the user. if PATCH_RAD > (min(source.mtheta, source.nzeta) - 1) ÷ 2 @@ -400,7 +387,7 @@ function compute_3D_kernel_matrices!( workspaces = [KernelWorkspace(PATCH_DIM, RAD_DIM, ANG_DIM) for _ in 1:max_threadid] # Parallel loop through observer points - Threads.@threads for idx_obs in 1:num_points + Threads.@threads for idx_obs in 1:n_obs # Get thread-local workspace ws = workspaces[Threads.threadid()] (; r_patch, dr_dθ_patch, dr_dζ_patch, r_polar, dr_dθ_polar, dr_dζ_polar, @@ -488,35 +475,12 @@ function compute_3D_kernel_matrices!( # Use the same normalization as in the 2D kernel so we can just add I to the diagonal # This makes the grri logic identical to the 2D kernel. grad_greenfunction_block ./= 2π - greenfunction ./= 2π + populate_greenfunction && (greenfunction ./= 2π) # Add the term that comes from the volume integral of Green's identity if typeof(source) == typeof(observer) - for i in 1:num_points + for i in 1:n_obs grad_greenfunction_block[i, i] += 1.0 end end end - -""" - kernel!(grad_greenfunction, greenfunction, observer, source, params::KernelParams3D) - -Dispatch wrapper for 3D kernel that forwards to `compute_3D_kernel_matrices!` with params. -""" -function kernel!( - grad_greenfunction::AbstractMatrix{Float64}, - greenfunction::AbstractMatrix{Float64}, - observer::Union{PlasmaGeometry3D,WallGeometry3D}, - source::Union{PlasmaGeometry3D,WallGeometry3D}, - params::KernelParams3D -) - return compute_3D_kernel_matrices!( - grad_greenfunction, - greenfunction, - observer, - source, - params.PATCH_RAD, - params.RAD_DIM, - params.INTERP_ORDER - ) -end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index e42946fb..f485b647 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -9,7 +9,7 @@ using AdaptiveArrayPools # Import parent modules import ..Equilibrium -using ..Utilities.FourierTransforms: compute_fourier_coefficients, fourier_transform!, fourier_inverse_transform! +using ..Utilities.FourierTransforms: FourierTransform, compute_fourier_coefficients include("Utilities.jl") include("DataTypes.jl") @@ -24,175 +24,233 @@ export extract_plasma_surface_at_psi export PlasmaGeometry """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) -Compute the vacuum response matrix and both Green's functions using provided vacuum inputs. +2D (axisymmetric) vacuum response calculation. -Single entry point for vacuum calculations. +Each toroidal mode `n` decouples in 2D geometry, so the routine loops over `inputs.n_modes`, +building the double-/single-layer operators, solving the exterior and interior systems, and +filling the corresponding diagonal block of the response matrix and the matching column block +of the Green's functions. +""" +@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) - - For **3D** (`inputs.nzeta > 1`), computes the full coupled response across all (m,n) modes defined - by `inputs.(mlow, mpert, nlow, npert)`. + mpert = length(inputs.m_modes) + num_points_surf = inputs.mtheta - - For **2D geometry** (`inputs.nzeta == 1`), supports either: + vac_data.wv .= 0 - + **single-n** (`inputs.npert == 1`): computes (m,n) response for `n = inputs.nlow` - + **multi-n** (`inputs.npert > 1`): loops over `n = inputs.nlow:(inputs.nlow+inputs.npert-1)` and returns - **blocks** of the full response matrices with one block per toroidal mode number. + # Form the plasma and wall geometries + plasma_surf = PlasmaGeometry(inputs) + wall = WallGeometry(inputs, plasma_surf, wall_settings) -This is the pure Julia implementation that replaces the Fortran `mscvac` function. -It computes both interior (grri) and exterior (grre) Green's functions for GPEC response calculations. + # Loop over all decoupled toroidal modes + for (idx_n, n) in enumerate(inputs.n_modes) + ft = FourierTransform(inputs.mtheta, mpert, inputs.m_modes[1]; n=n, ν=plasma_surf.ν) -# Arguments + # Diagonal block of wv and matching column block of the Green's functions + block_idx = ((idx_n-1)*mpert+1):(idx_n*mpert) + wv_block = @view vac_data.wv[block_idx, block_idx] + grri_block = @view vac_data.grri[:, block_idx] + grre_block = @view vac_data.grre[:, block_idx] - - `inputs`: `VacuumInput` struct with mode numbers, grid resolution, and boundary info. - - `wall_settings::WallShapeSettings`: Wall geometry configuration. + # Active rows for computation (plasma only if no wall, plasma+wall if wall present) + num_points_total = wall.nowall ? num_points_surf : 2 * num_points_surf -# Returns + # Local work matrices + grad_green = zeros!(pool, num_points_total, num_points_total) + green_temp = zeros!(pool, num_points_surf, num_points_surf) + + # Views into output Green's function matrices for the active rows/columns + grre = @view grre_block[1:num_points_total, :] + grri = @view grri_block[1:num_points_total, :] - - `wv`: Complex vacuum response matrix. + # Plasma–Plasma block + compute_2D_kernel_matrices!(grad_green, green_temp, plasma_surf, plasma_surf, n) - + 2D single-n: `mpert × mpert` - + 2D multi-n: `(mpert*npert) × (mpert*npert)` (block diagonal) - + 3D: `num_modes × num_modes` (full coupled) + # Project plasma observer onto source basis exp(i*(mθ - nν)) + mul!(view(grre, 1:num_points_surf, :), green_temp, ft.basis') - - `grri`: Interior Green's function matrix. + if !wall.nowall + # Plasma–Wall block + compute_2D_kernel_matrices!(grad_green, green_temp, plasma_surf, wall, n) + # Wall–Wall block + 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ν)) + mul!(view(grre, (num_points_surf+1):num_points_total, :), green_temp, ft.basis') + end - - `grre`: Exterior Green's function matrix. + # Compute both Green's functions: exterior (kernelsign=+1) then interior (kernelsign=-1) + grri .= grre # start from same as exterior + grad_green_interior = similar!(pool, grad_green) + grad_green_interior .= grad_green - - `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!( - wv::AbstractMatrix{ComplexF64}, - grri_in::AbstractMatrix{Float64}, - grre_in::AbstractMatrix{Float64}, - plasma_pts::AbstractMatrix{Float64}, - wall_pts::AbstractMatrix{Float64}, - inputs::VacuumInput, - wall_settings::WallShapeSettings; - n_override::Union{Nothing,Int}=nothing -) - - # Initialize surface geometries - plasma_surf = inputs.nzeta > 1 ? PlasmaGeometry3D(inputs) : PlasmaGeometry(inputs) - wall = inputs.nzeta > 1 ? WallGeometry3D(inputs, wall_settings) : WallGeometry(inputs, plasma_surf, wall_settings) - - # Compute Fourier basis coefficients - ν = 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 - if inputs.nzeta > 1 - # Hardcode these values for now - can expose to the user in the future - kparams = KernelParams3D(11, 20, 5) - else - kparams = KernelParams2D(n_override) + # Solve exterior first, overwriting grad_green to save memory since we already have the interior kernel + ldiv!(lu!(grad_green), grre) + + # Interior flips the sign of the normal, but not the diagonal terms, so we multiply by -1 and add 2I to the diagonal + grad_green_interior .*= -1 + for i in 1:num_points_total + grad_green_interior[i, i] += 2.0 + end + ldiv!(lu!(grad_green_interior), grri) + + # 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 - # Active rows for computation (plasma only if no wall, plasma+wall if wall present) - num_points_total = wall.nowall ? num_points_surf : 2 * num_points_surf + inputs.force_wv_symmetry && hermitianpart!(vac_data.wv) + + # Populate coordinate arrays + @views begin + vac_data.plasma_pts[:, 1] .= plasma_surf.x + vac_data.plasma_pts[:, 2] .= 0.0 + vac_data.plasma_pts[:, 3] .= plasma_surf.z + vac_data.wall_pts[:, 1] .= wall.x + vac_data.wall_pts[:, 2] .= 0.0 + vac_data.wall_pts[:, 3] .= wall.z + end +end + +""" + _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + +3D (`inputs.nzeta > 1`) vacuum response via block-circulant field-period reduction. For `nfp == 1` +the block-circulant assembly and residue-class loop are skipped in favour of a more efficient +direct real solve on the built kernel matrices (equivalent to 2D). + +The `nfp`-periodic boundary makes the single-/double-layer operators `S`, `D` block-circulant in +the field-period index. Writing the response as `wv = (4π²/N)·Eᴴ·D⁻¹S·E` (`E` the complex Fourier +basis, `N = mtheta·nzeta_full`), the structure block-diagonalizes the problem by toroidal residue +class `k = mod(n, nfp)`: modes with different `k` do not couple, and within a class + + D̂ₖ = Σ_d D_d ω^{k d}, Ŝₖ = Σ_d S_d ω^{k d}, ω = exp(-2πi/nfp), - # Local work arrays - grad_green = zeros!(pool, num_points_total, num_points_total) - green_temp = zeros!(pool, num_points_surf, num_points_surf) +with `D_d`, `S_d` the blocks coupling observers in field period 0 to sources in period `d`. Each +class needs one solve `wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \\ Ŝₖ)|_plasma·E_local`. Only the first +block-row of the operators is built, so the kernel cost drops by `nfp` and the dense `O(N³)` +factorization is replaced by per-class `O(M³)` solves (`M = N/nfp`). - # Views into output Green's function matrices for the active rows/columns - grre = @view grre_in[1:num_points_total, :] - grri = @view grri_in[1:num_points_total, :] +Only `wv` is produced currently; `grri`/`grre` are returned zeroed and are not yet supported in the 3D path. +Extension point: per residue class, apply the per-period basis to `D̂ₖ⁻¹Ŝₖ` for the exterior columns and the +interior variant `-D + 2I` for the interior columns, then scatter back into the `[2N × 2·num_modes]` arrays. +""" +@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + + (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs + fill!(vac_data.wv, 0) + + # Full-torus geometry for source surface; observers are restricted to one field period + full = expand_field_periods(inputs) + plasma_surf = PlasmaGeometry3D(full) + wall = WallGeometry3D(full, wall_settings) + + num_points_per_fp = mtheta * nzeta # points per field period + num_points = num_points_per_fp * nfp # full-torus point count + mpert = length(m_modes) + nb = wall.nowall ? 1 : 2 # surface blocks: plasma, or [plasma; wall] + n_obs = nb * num_points_per_fp # number of observer points per field period + + # Complex Fourier basis exp(-i(mθ-nζ)) on a single field period + exp_mn_basis = compute_fourier_coefficients(mtheta, m_modes, nzeta * nfp, n_modes; nfp=nfp) - # Plasma–Plasma block - kernel!(grad_green, green_temp, plasma_surf, plasma_surf, kparams) + # Local work matrices + grad_green = zeros!(pool, n_obs, nb * num_points) # double-layer D (with +I on the period-0 diagonals) + green_temp = zeros!(pool, n_obs, num_points) # single-layer S (plasma sources only) - # Fourier transform obs=plasma, src=plasma block - fourier_transform!(grre, green_temp, cos_mn_basis) - fourier_transform!(grre, green_temp, sin_mn_basis; col_offset=num_modes) + # Kernel parameters, hardcoded for now + PATCH_RAD = 11 + RAD_DIM = 20 + INTERP_ORDER = 5 + + # Plasma–plasma block + compute_3D_kernel_matrices!(grad_green, view(green_temp, 1:num_points_per_fp, :), plasma_surf, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER) if !wall.nowall # Plasma–Wall block - kernel!(grad_green, green_temp, plasma_surf, wall, kparams) - # Wall–Wall block - kernel!(grad_green, green_temp, wall, wall, kparams) + compute_3D_kernel_matrices!(grad_green, view(green_temp, 1:num_points_per_fp, :), plasma_surf, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER) # Wall–Plasma block - kernel!(grad_green, green_temp, wall, plasma_surf, kparams) - # Fourier transform obs=wall, src=plasma block - fourier_transform!(grre, green_temp, cos_mn_basis; row_offset=num_points_surf) - fourier_transform!(grre, green_temp, sin_mn_basis; row_offset=num_points_surf, col_offset=num_modes) + compute_3D_kernel_matrices!(grad_green, view(green_temp, (num_points_per_fp+1):(2*num_points_per_fp), :), wall, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER) + # Wall–Wall block + compute_3D_kernel_matrices!(grad_green, view(green_temp, (num_points_per_fp+1):(2*num_points_per_fp), :), wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER) end - # Compute both Green's functions: exterior (kernelsign=+1) then interior (kernelsign=-1) - grri .= grre # start from same as exterior - grad_green_interior = similar!(pool, grad_green) - grad_green_interior .= grad_green - - # Solve exterior first, overwriting grad_green to save memory since we already have the interior kernel - F_ext = lu!(grad_green) - ldiv!(F_ext, grre) - - # Interior flips the sign of the normal, but not the diagonal terms, so we multiply by -1 and add 2I to the diagonal - grad_green_interior .*= -1 - for i in 1:num_points_total - grad_green_interior[i, i] += 2.0 - end - F_int = lu!(grad_green_interior) - ldiv!(F_int, grri) - - # Perform inverse Fourier transforms to get response matrix components [Chance Phys. Plasmas 2007 052506 eq. 115-118] - arr, aii, ari, air = ntuple(_ -> zeros(num_modes, num_modes), 4) - fourier_inverse_transform!(arr, grre, cos_mn_basis) - fourier_inverse_transform!(aii, grre, sin_mn_basis; col_offset=num_modes) - fourier_inverse_transform!(ari, grre, sin_mn_basis) - fourier_inverse_transform!(air, grre, cos_mn_basis; col_offset=num_modes) - - # fourier_inverse_transform! uses 1/N normalization; restore the 4π² physics factor - # from the vacuum Green's function integral [Chance 2007 eq. 114-118] - arr .*= 4π^2 - aii .*= 4π^2 - ari .*= 4π^2 - air .*= 4π^2 - - # Final form of vacuum response matrix [Chance Phys. Plasmas 2007 052506 eq. 114] - wv .= complex.(arr .+ aii, air .- ari) - inputs.force_wv_symmetry && hermitianpart!(wv) - - # Fill coordinate arrays - if inputs.nzeta > 1 # 3D - plasma_pts .= plasma_surf.r - wall_pts .= wall.r - else # 2D - @views begin - plasma_pts[:, 1] .= plasma_surf.x - plasma_pts[:, 2] .= 0.0 - plasma_pts[:, 3] .= plasma_surf.z - wall_pts[:, 1] .= wall.x - wall_pts[:, 2] .= 0.0 - wall_pts[:, 3] .= wall.z + if nfp == 1 + # nfp == 1 case: direct real solve on the built kernel matrices (equivalent to 2D) + ldiv!(lu!(grad_green), green_temp) + G_plasma = @view green_temp[1:num_points_per_fp, :] + vac_data.wv .= (4π^2 / num_points_per_fp) .* (exp_mn_basis * (G_plasma * exp_mn_basis')) + else + # nfp > 1 case: block-circulant solve on the built kernel matrices + # Solve one reduced (nb·M)×(nb·M) system per toroidal residue class k = mod(n, nfp) + D̂ = zeros!(pool, ComplexF64, n_obs, n_obs) + Ŝ = zeros!(pool, ComplexF64, n_obs, num_points_per_fp) + for k in unique(mod.(n_modes, nfp)) + # Reduced operators D̂ₖ = Σ_d D_d ω^{k d}, Ŝₖ = Σ_d S_d ω^{k d} (fused, allocation-free). + # D̂ keeps the nb-block source structure (plasma cols 1:M, wall cols M+1:2M); Ŝ is plasma- + # source-only over all nb·M observer rows. + fill!(D̂, 0) + fill!(Ŝ, 0) + for d in 0:(nfp-1) + phase = cis(-2π * (k * d) / nfp) + for s in 0:(nb-1) + src_cols = (s*num_points+d*num_points_per_fp+1):(s*num_points+(d+1)*num_points_per_fp) + dst_cols = (s*num_points_per_fp+1):((s+1)*num_points_per_fp) + @views @. D̂[:, dst_cols] += phase * grad_green[:, src_cols] + end + scols = (d*num_points_per_fp+1):((d+1)*num_points_per_fp) + @views @. Ŝ += phase * green_temp[:, scols] + end + + # Exterior response operator for this sector: solve D̂ₖ G = Ŝₖ in place (Ŝ ← G = D̂ₖ⁻¹Ŝₖ) + ldiv!(lu!(D̂), Ŝ) + mode_cols = [(idx_m + (idx_n-1)*mpert) for (idx_n, n) in enumerate(n_modes) if mod(n, nfp) == k for idx_m in 1:mpert] + E = @view exp_mn_basis[mode_cols, :] + G_plasma = @view Ŝ[1:num_points_per_fp, :] # plasma-observer rows of the combined exterior operator + vac_data.wv[mode_cols, mode_cols] .= (4π^2 / num_points_per_fp) .* (E * (G_plasma * E')) end end + + inputs.force_wv_symmetry && hermitianpart!(vac_data.wv) + + # Zero out the Green's function matrices (not tested in 3D yet) + fill!(vac_data.grri, 0) + fill!(vac_data.grre, 0) + + # Populate coordinate arrays + vac_data.plasma_pts .= plasma_surf.r + vac_data.wall_pts .= wall.r end """ compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) -Allocate and return the vacuum response matrix and Green's functions for the given -vacuum inputs. +Allocate and return the vacuum response matrix and Green's functions for the given vacuum +inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the +output arrays for the full torus and forwards to the same 2D/3D workers. For performance-critical +paths that already own preallocated storage (e.g. `ForceFreeStates.VacuumData`), prefer the +in-place method to avoid extra heap allocations. -This is a thin allocating wrapper around the in‑place [`compute_vacuum_response!`] -implementation. For performance‑critical paths that already own preallocated storage -(e.g. `ForceFreeStates.VacuumData`), prefer the in‑place method to avoid extra -heap allocations. +# Returns + + - `wv`: complex vacuum response matrix (`num_modes × num_modes`). + - `grri`, `grre`: interior/exterior Green's functions (zeroed on the 3D nowall path). + - `plasma_pts`, `wall_pts`: surface coordinate arrays. """ -@with_pool pool function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) - # Allocate storage for the vacuum response matrix and Green's functions - numpoints = inputs.mtheta * inputs.nzeta - num_modes = inputs.mpert * inputs.npert + num_points = inputs.mtheta * inputs.nzeta * inputs.nfp # mtheta for 2D + num_modes = length(inputs.m_modes) * length(inputs.n_modes) vac = ( wv=zeros(ComplexF64, num_modes, num_modes), - grri=zeros(2 * numpoints, 2 * num_modes), - grre=zeros(2 * numpoints, 2 * num_modes), - plasma_pts=zeros(numpoints, 3), - wall_pts=zeros(numpoints, 3) + grri=zeros(ComplexF64, 2 * num_points, num_modes), + grre=zeros(ComplexF64, 2 * num_points, num_modes), + plasma_pts=zeros(num_points, 3), + wall_pts=zeros(num_points, 3) ) compute_vacuum_response!(vac, inputs, wall_settings) @@ -202,63 +260,27 @@ end """ compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) -In-place variant that computes the vacuum response and directly populates the arrays -stored in `vac_data`. +In-place variant that computes the vacuum response and directly populates the arrays stored in +`vac_data`. Dispatches on dimensionality only: 2D (`inputs.nzeta == 1`) routes to +[`_compute_vacuum_response_2d!`], 3D to [`_compute_vacuum_response_3d!`]. -The `vac_data` argument is expected to provide the following writable fields with -compatible sizes: +The `vac_data` argument is expected to provide the following writable fields with compatible +sizes: - `wv::AbstractMatrix{ComplexF64}` – vacuum response matrix - - `grri::AbstractMatrix{Float64}` – interior Green's functions - - `grre::AbstractMatrix{Float64}` – exterior Green's functions + - `grri::AbstractMatrix{ComplexF64}` – interior Green's functions + - `grre::AbstractMatrix{ComplexF64}` – exterior Green's functions - `plasma_pts::AbstractMatrix{Float64}` – plasma surface coordinates - `wall_pts::AbstractMatrix{Float64}` – wall surface coordinates -This is designed to work with `ForceFreeStates.VacuumData` but does not depend on -its concrete type (duck-typed on field names only). +This is designed to work with `ForceFreeStates.VacuumData` but does not depend on its concrete +type (duck-typed on field names only). """ function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) - - mpert = inputs.mpert - npert = inputs.npert - numpert_total = mpert * npert - - # 3D vacuum: full coupled response across (m,n) from a single kernel call - if inputs.nzeta > 1 - _compute_vacuum_response_single!( - vac_data.wv, - vac_data.grri, - vac_data.grre, - vac_data.plasma_pts, - vac_data.wall_pts, - inputs, - wall_settings - ) + if inputs.nzeta == 1 + _compute_vacuum_response_2d!(vac_data, inputs, wall_settings) else - # 2D vacuum: fill diagonal blocks of the response matrix - vac_data.wv .= 0 - - # Each n is independent in 2D geometry → fill diagonal blocks inside preallocated arrays - ns = inputs.nlow:(inputs.nlow+inputs.npert-1) - for (idx_n, n) in enumerate(ns) - block_idx = ((idx_n-1)*mpert+1):(idx_n*mpert) - cols = vcat(block_idx, numpert_total .+ block_idx) - - wv_block = @view vac_data.wv[block_idx, block_idx] - grri_block = @view vac_data.grri[:, cols] - grre_block = @view vac_data.grre[:, cols] - - _compute_vacuum_response_single!( - wv_block, - grri_block, - grre_block, - vac_data.plasma_pts, - vac_data.wall_pts, - inputs, - wall_settings; - n_override=n - ) - end + _compute_vacuum_response_3d!(vac_data, inputs, wall_settings) end end diff --git a/test/runtests_fouriertransforms.jl b/test/runtests_fouriertransforms.jl index a6fd77f2..005d66d6 100644 --- a/test/runtests_fouriertransforms.jl +++ b/test/runtests_fouriertransforms.jl @@ -1,20 +1,13 @@ # Unit tests for the high-level FourierTransform functor interface in -# src/Utilities/FourierTransforms.jl. The low-level fourier_transform! / -# fourier_inverse_transform! matrix kernels are already covered in -# runtests_utilities.jl; this file targets the functor, inverse, in-place, -# and compute_fourier_coefficients paths. -# -# Expected mode amplitudes follow from grid orthogonality under the module's -# exp(-imθ), 1/N forward convention (Fortran iscdftf): on the uniform grid -# θ_i = 2π(i-1)/N, cos(m0θ) → 0.5 at ±m0, sin(m0θ) → ∓0.5im at ±m0. - -using GeneralizedPerturbedEquilibrium.Utilities: FourierTransform, inverse, transform!, inverse_transform!, compute_fourier_coefficients +# src/Utilities/FourierTransforms.jl. + +using GeneralizedPerturbedEquilibrium.Utilities: FourierTransform, inverse, compute_fourier_coefficients using GeneralizedPerturbedEquilibrium.Utilities: empty_FourierCoefficients const ATOL = 1e-10 # analytic-value comparisons const ATOL_TIGHT = 1e-12 # in-place vs allocating / formula-vs-formula checks -@testset "FourierTransform functor" begin +@testset "FourierTransform" begin atol = ATOL # Mode index l for mode number m (m = mlow + l - 1). @@ -82,92 +75,34 @@ const ATOL_TIGHT = 1e-12 # in-place vs allocating / formula-vs-formula checks @test isapprox(modes[midx(5, mlow), 2], -0.5im; atol) @test isapprox(modes[midx(-5, mlow), 2], 0.5im; atol) end - - @testset "real-mode inverse is a pure cosine expansion" begin - N, mlow, mpert = 32, -4, 9 - ft = FourierTransform(N, mpert, mlow) - θ = collect(range(; start=0, length=N, step=2π/N)) - - # A single real mode at m reconstructs (2π/N)·cos(mθ). - modes = zeros(Float64, mpert) - modes[midx(2, mlow)] = 1.0 - recon = inverse(ft, modes) - @test all(isapprox.(recon, (2π / N) .* cos.(2 .* θ); atol)) - end -end - -@testset "FourierTransform in-place == allocating" begin - atol = ATOL_TIGHT - N, mlow, mpert = 32, -6, 13 - ft = FourierTransform(N, mpert, mlow) - θ = collect(range(; start=0, length=N, step=2π/N)) - - @testset "transform! (vector)" begin - data_r = cos.(2 .* θ) .+ 0.3 .* sin.(4 .* θ) - out = zeros(ComplexF64, mpert) - transform!(out, ft, data_r) - @test all(isapprox.(out, ft(data_r); atol)) - - data_c = data_r .+ im .* sin.(3 .* θ) - transform!(out, ft, data_c) - @test all(isapprox.(out, ft(data_c); atol)) - end - - @testset "transform! (matrix)" begin - data = hcat(cos.(2 .* θ), sin.(3 .* θ), cos.(5 .* θ)) - out = zeros(ComplexF64, mpert, size(data, 2)) - transform!(out, ft, data) - @test all(isapprox.(out, ft(data); atol)) - - datac = ComplexF64.(data) .+ im - transform!(out, ft, datac) - @test all(isapprox.(out, ft(datac); atol)) - end - - @testset "inverse_transform!" begin - modes_v = ComplexF64[Float64(l) - im * l for l in 1:mpert] - out_v = zeros(ComplexF64, N) - inverse_transform!(out_v, ft, modes_v) - @test all(isapprox.(out_v, inverse(ft, modes_v); atol)) - - modes_m = reshape(modes_v, mpert, 1) .* [1.0 -2.0im] - out_m = zeros(ComplexF64, N, 2) - inverse_transform!(out_m, ft, modes_m) - @test all(isapprox.(out_m, inverse(ft, modes_m); atol)) - end end @testset "compute_fourier_coefficients" begin atol = ATOL_TIGHT - @testset "2D basis is cos(mθ - nν)/sin(mθ - nν)" begin - # Nonzero n and ν exercise the general phase-shifted form and lock both the - # grid convention (start=0, step=2π/N) and the -n·ν sign. + @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)) - cosb, sinb = compute_fourier_coefficients(N, mpert, mlow, 1, 1, 1; n_2D=n, ν=ν) - @test size(cosb) == (N, mpert) - @test size(sinb) == (N, mpert) + 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(mlow:(mlow+mpert-1)) - @test all(isapprox.(cosb[:, l], cos.(m .* θ .- n .* ν); atol)) - @test all(isapprox.(sinb[:, l], sin.(m .* θ .- n .* ν); atol)) + for (l, m) in enumerate(m_modes) + expected = exp.(-im .* (m .* θ .- n .* ν)) + @test all(isapprox.(basis[l, :], expected; atol)) end - # The FourierTransform constructor stores exactly the n=0, ν=0 basis. ft = FourierTransform(N, mpert, mlow) - @test ft.cslth == compute_fourier_coefficients(N, mpert, mlow, 1, 1, 1; n_2D=0, ν=zeros(N))[1] - @test ft.snlth == compute_fourier_coefficients(N, mpert, mlow, 1, 1, 1; n_2D=0, ν=zeros(N))[2] + @test ft.basis == compute_fourier_coefficients(N, m_modes, 0, zeros(N)) end @testset "3D basis shapes" begin mtheta, mpert, mlow = 8, 5, -2 nzeta, npert, nlow = 6, 3, -1 - cosb, sinb = compute_fourier_coefficients(mtheta, mpert, mlow, nzeta, npert, nlow) - @test size(cosb) == (mtheta * nzeta, mpert * npert) - @test size(sinb) == (mtheta * nzeta, mpert * npert) + basis = compute_fourier_coefficients(mtheta, mlow:(mlow+mpert-1), nzeta, nlow:(nlow+npert-1)) + @test size(basis) == (mpert * npert, mtheta * nzeta) end end diff --git a/test/runtests_utilities.jl b/test/runtests_utilities.jl index 43fa8d78..bd18bc40 100644 --- a/test/runtests_utilities.jl +++ b/test/runtests_utilities.jl @@ -29,50 +29,4 @@ GeneralizedPerturbedEquilibrium.Utilities.get_complex_coeffs!(out, fc, 10, 1) @test out[3] == c2 # Mode 2 is at index 3 (0-indexed mode) end - - @testset "FourierTransforms" begin - @info "Testing fourier_transform! and fourier_inverse_transform! from Utilities module" - using GeneralizedPerturbedEquilibrium.Utilities: fourier_transform!, fourier_inverse_transform! - - @testset "fourier_transform!" begin - mtheta, mpert = 4, 2 - gil = zeros(mtheta, mpert) - gij = [1.0 2.0 3.0 4.0; 5.0 6.0 7.0 8.0; 9.0 10.0 11.0 12.0; 13.0 14.0 15.0 16.0] - cs = zeros(mtheta, mpert) - cs[:, 1] .= 1.0 # First mode: all ones → output is column sum of gij - - fourier_transform!(gil, gij, cs; row_offset=0, col_offset=0) - - # gil[:, 1] = gij * cs[:, 1] = sum over rows of gij per column → columns of gij summed - # Column 1 of gij: [1,5,9,13] → sum 28; column 2: [2,6,10,14] → 32; etc. No, mul! does (gil block) = gij * cs. - # So gil[i, l] = sum_j gij[i,j] * cs[j,l]. For cs[:,1]=1: gil[:,1] = gij * [1,1,1,1]' = row sums of gij. - # Rows of gij: [1,2,3,4] sum 10, [5,6,7,8] sum 26, [9,10,11,12] sum 42, [13,14,15,16] sum 58. - @test gil[:, 1] == [10.0, 26.0, 42.0, 58.0] - @test gil[:, 2] == [0.0, 0.0, 0.0, 0.0] - - # Test with column offset: fill second block - gil2 = zeros(mtheta, 4) - cs2 = zeros(mtheta, 2) - cs2[:, 1] .= 1.0 - fourier_transform!(gil2, gij, cs2; col_offset=2) - @test gil2[:, 3] == [10.0, 26.0, 42.0, 58.0] - @test all(gil2[:, 1:2] .== 0) && all(gil2[:, 4] .== 0) - end - - @testset "fourier_inverse_transform!" begin - mtheta, mpert = 4, 2 - gil = zeros(mtheta, mpert) - gil[:, 1] = [1.0, 2.0, 3.0, 4.0] - cs = zeros(mtheta, mpert) - cs[:, 1] .= 1.0 - gll = zeros(mpert, mpert) - - fourier_inverse_transform!(gll, gil, cs) - - # gll = (1/mtheta) * cs' * gil → gll[1,1] = (1/4) * sum_i cs[i,1]*gil[i,1] = 10/4 - expected = 10.0 / mtheta - @test isapprox(gll[1, 1], expected) - @test isapprox(gll[2, 1], 0.0) - end - end end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 27c92a77..52a73494 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -11,22 +11,20 @@ @test vac.ν == Float64[] @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.m_modes == [1] + @test vac.n_modes == [1] @test vac.mtheta == 1 @test vac.nzeta == 1 @test vac.force_wv_symmetry == true end @testset "keyword constructor" begin - vac = VacuumInput(mtheta=32, mpert=3, mlow=1, nlow=2, npert=2, nzeta=1) + vac = VacuumInput(mtheta=32, m_modes=[1, 2, 3], n_modes=[2, 3], nzeta=1) @test vac.mtheta == 32 - @test vac.mpert == 3 - @test vac.mlow == 1 - @test vac.nlow == 2 - @test vac.npert == 2 + @test length(vac.m_modes) == 3 + @test vac.m_modes[1] == 1 + @test vac.n_modes[1] == 2 + @test length(vac.n_modes) == 2 @test vac.nzeta == 1 end end @@ -56,10 +54,8 @@ z=[0.0, 0.1, 0.0, -0.1, 0.0], ν=zeros(5), mtheta=5, - mpert=1, - mlow=1, - nlow=1, - npert=1, + m_modes=[1], + n_modes=[1], nzeta=1 ) surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(inputs) @@ -80,10 +76,8 @@ nzeta_in=1, ν=zeros(5), mtheta=8, - mpert=1, - mlow=0, - nlow=0, - npert=1, + m_modes=[0], + n_modes=[0], nzeta=1 ) surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(inputs) @@ -360,16 +354,14 @@ end @testset "compute_vacuum_response" begin - _make_inputs(; mtheta=128, mtheta_eq=17, mpert=2, nlow=1, npert=1) = VacuumInput( + _make_inputs(; mtheta=128, mtheta_eq=17, m_modes=1:2, n_modes=[1]) = VacuumInput( 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, - mpert=mpert, - nlow=nlow, - npert=npert, + m_modes=collect(Int, m_modes), + n_modes=collect(Int, n_modes), nzeta=1, mtheta=mtheta ) @@ -380,13 +372,13 @@ wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta - num_modes = inputs.mpert * inputs.npert + num_modes = length(inputs.m_modes) * length(inputs.n_modes) @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 size(grri) == (2 * numpoints, num_modes) + @test size(grre) == (2 * numpoints, num_modes) @test all(isfinite, grri) @test all(isfinite, grre) @test size(plasma_pts) == (numpoints, 3) @@ -403,9 +395,9 @@ wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta - num_modes = inputs.mpert * inputs.npert + num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test size(grri) == (2 * numpoints, 2 * num_modes) + @test size(grri) == (2 * numpoints, num_modes) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # plasma_pts layout: col1=R, col2=0, col3=Z @@ -415,12 +407,12 @@ end @testset "edge: single poloidal mode mpert=1" begin - inputs = _make_inputs(mpert=1, npert=1) + inputs = _make_inputs(m_modes=[1], n_modes=[1]) wall_settings = WallShapeSettings(shape="nowall") wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (1, 1) @test all(isfinite, wv) - @test size(grri, 2) == 2 # 2 * num_modes with num_modes=1 + @test size(grri, 2) == 1 end @testset "edge: small mtheta" begin @@ -429,7 +421,7 @@ wall_settings = WallShapeSettings(shape="nowall") wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (2, 2) - @test size(grri) == (32, 4) # 2*16, 2*2 + @test size(grri) == (32, 2) # 2*mtheta, num_modes=2 @test size(plasma_pts) == (16, 3) end @@ -441,8 +433,8 @@ wv, grri, grre, pp, wp = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta - num_modes = inputs.mpert * inputs.npert - vac = (wv=zeros(ComplexF64, num_modes, num_modes), grri=zeros(2 * numpoints, 2 * num_modes), grre=zeros(2 * numpoints, 2 * num_modes), + num_modes = length(inputs.m_modes) * length(inputs.n_modes) + vac = (wv=zeros(ComplexF64, num_modes, num_modes), grri=zeros(ComplexF64, 2 * numpoints, num_modes), grre=zeros(ComplexF64, 2 * numpoints, num_modes), plasma_pts=zeros(numpoints, 3), wall_pts=zeros(numpoints, 3)) compute_vacuum_response!(vac, inputs, wall_settings) @@ -476,18 +468,16 @@ end # 3D vacuum: nzeta > 1, full (m,n) coupling, PlasmaGeometry3D, WallGeometry3D - # Kernel requires mtheta, nzeta >= PATCH_DIM (23 for default KernelParams3D(11, 20, 5)) + # Kernel requires mtheta, nzeta >= PATCH_DIM (23 for default PATCH_RAD=11) @testset "Vacuum.jl (3D)" begin - _make_3d_inputs(; mtheta=32, mtheta_eq=17, mpert=2, nlow=0, npert=2, nzeta=32) = VacuumInput( + _make_3d_inputs(; mtheta=32, mtheta_eq=17, m_modes=1:2, n_modes=0:1, nzeta=32) = VacuumInput( 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, - mpert=mpert, - nlow=nlow, - npert=npert, + m_modes=collect(Int, m_modes), + n_modes=collect(Int, n_modes), nzeta=nzeta, mtheta=mtheta ) @@ -520,28 +510,19 @@ z=vec(Z), mtheta_in=mtheta_in, nzeta_in=nzeta_in, - mlow=1, - mpert=mpert, - nlow=nlow, - npert=npert, + m_modes=collect(1:mpert), + n_modes=collect(nlow:(nlow+npert-1)), mtheta=mtheta, nzeta=nzeta ) end @testset "VacuumInput nzeta > 1" begin - vac = VacuumInput(mtheta=32, nzeta=24, mpert=2, npert=2) + vac = VacuumInput(mtheta=32, nzeta=24, m_modes=[1, 2], n_modes=[1, 2]) @test vac.nzeta == 24 @test vac.mtheta == 32 end - @testset "KernelParams3D" begin - params = GeneralizedPerturbedEquilibrium.Vacuum.KernelParams3D(11, 20, 5) - @test params.PATCH_RAD == 11 - @test params.RAD_DIM == 20 - @test params.INTERP_ORDER == 5 - end - @testset "PlasmaGeometry3D" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs) @@ -604,12 +585,12 @@ wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta - num_modes = inputs.mpert * inputs.npert + num_modes = length(inputs.m_modes) * length(inputs.n_modes) @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 size(grri) == (2 * numpoints, num_modes) + @test size(grre) == (2 * numpoints, num_modes) @test all(isfinite, grri) @test all(isfinite, grre) @test size(plasma_pts) == (numpoints, 3) @@ -627,13 +608,13 @@ wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta - num_modes = inputs.mpert * inputs.npert + num_modes = length(inputs.m_modes) * length(inputs.n_modes) @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 size(grri) == (2 * numpoints, num_modes) + @test size(grre) == (2 * numpoints, num_modes) @test all(isfinite, grri) @test all(isfinite, grre) @test size(plasma_pts) == (numpoints, 3) @@ -648,9 +629,9 @@ wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta - num_modes = inputs.mpert * inputs.npert + num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test size(grri) == (2 * numpoints, 2 * num_modes) + @test size(grri) == (2 * numpoints, num_modes) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # Wall and plasma should differ (conformal wall offset from plasma) @@ -658,6 +639,85 @@ @test isapprox(wv, wv', rtol=1e-12) end + # Field-periodic (layer-2) reduction: an nfp-periodic boundary makes the boundary-integral + # operators block-circulant, so the reduced per-residue-class solve must reproduce the full + # torus result (the n_stride=1 bridge case spans several residue classes mod nfp). + @testset "compute_vacuum_response 3D field-periodic reduction" begin + # Build one field period of a 3D boundary on the full-torus ζ spacing; expand_field_periods + # tiles it by rigid rotation, so reduced and full paths share an identical full-torus surface. + _make_periodic_period(; mtheta, nzeta_p, nfp, R0=1.7, a=0.3, b=0.3, ε=0.04) = begin + θ = range(; start=0, length=mtheta, step=2π/mtheta) + X = zeros(mtheta, nzeta_p) + Y = zeros(mtheta, nzeta_p) + Z = zeros(mtheta, nzeta_p) + for j in 1:nzeta_p + ζ = (j - 1) * 2π / (nzeta_p * nfp) # period 0 of the full-torus grid + for (i, θi) in enumerate(θ) + R = R0 + a * cos(θi) + ε * cos(θi + nfp * ζ) + X[i, j] = R * cos(ζ) + Y[i, j] = R * sin(ζ) + Z[i, j] = b * sin(θi) + ε * sin(θi + nfp * ζ) + end + end + return vec(X), vec(Y), vec(Z) + end + + mtheta, nzeta_p, nfp = 24, 8, 3 # full torus nzeta = 24 ≥ PATCH_DIM (23) + m_modes = collect(-1:1) + n_modes = [1, 2, 3, 4] # residues {1, 2, 0, 1} mod nfp exercise grouping + Xp, Yp, Zp = _make_periodic_period(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=nfp) + + inputs_red = VacuumInput( + x=Xp, y=Yp, z=Zp, + mtheta_in=mtheta, nzeta_in=nzeta_p, + m_modes=m_modes, n_modes=n_modes, + mtheta=mtheta, nzeta=nzeta_p, + force_wv_symmetry=false, nfp=nfp + ) + wall_settings = WallShapeSettings(shape="nowall") + + # Reduced (block-circulant) path + wv_red, _, _, plasma_pts_red, _ = compute_vacuum_response(inputs_red, wall_settings) + + # Full-torus reference: pre-expand so nfp=1 forces the dense path + inputs_full = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(inputs_red) + @test inputs_full.nfp == 1 + @test inputs_full.nzeta == nzeta_p * nfp + wv_full, _, _, _, _ = compute_vacuum_response(inputs_full, wall_settings) + + num_modes = length(m_modes) * length(n_modes) + @test size(wv_red) == (num_modes, num_modes) + @test eltype(wv_red) == ComplexF64 + @test all(isfinite, wv_red) + # plasma points are the full torus + @test size(plasma_pts_red, 1) == mtheta * nzeta_p * nfp + + # Reduced == dense full-torus result (block-circulant solve is an exact reorganization) + @test isapprox(wv_red, wv_full; rtol=1e-6, atol=1e-7) + + # Modes in different residue classes mod nfp do not couple in the reduced result + classes = mod.(n_modes, nfp) + mpert = length(m_modes) + for in1 in eachindex(n_modes), in2 in eachindex(n_modes) + if classes[in1] != classes[in2] + rows = ((in1-1)*mpert+1):(in1*mpert) + cols = ((in2-1)*mpert+1):(in2*mpert) + @test all(iszero, wv_red[rows, cols]) + end + end + + # Hermitian symmetrization still works through the reduced path + inputs_sym = VacuumInput( + x=Xp, y=Yp, z=Zp, + mtheta_in=mtheta, nzeta_in=nzeta_p, + m_modes=m_modes, n_modes=n_modes, + mtheta=mtheta, nzeta=nzeta_p, + force_wv_symmetry=true, nfp=nfp + ) + wv_sym, _, _, _, _ = compute_vacuum_response(inputs_sym, wall_settings) + @test isapprox(wv_sym, wv_sym', rtol=1e-12) + end + @testset "Kernel3D laplace_kernel" begin G, K = GeneralizedPerturbedEquilibrium.Vacuum.laplace_kernel(1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 1.0, 0.0, 0.0) # Kernel returns 1/|r_obs - r_src| (4π factor applied elsewhere in BIE)