Utilities Module

Utilities

Shared mathematical and computational utilities for GPEC modules.

This module provides common functionality used across multiple GPEC modules, including efficient Fourier transforms, numerical integration, and other mathematical utilities.

Submodules

• FourierTransforms: Efficient Fourier transforms with pre-computed basis functions

FourierCoefficients

Lightweight container for Fourier coefficients without spline interpolation. Use this when you only need to access coefficients at original grid points.

Fields

• xs::Vector{Float64}: Radial coordinates
• mband::Int: Number of Fourier modes (0 to mband inclusive)
• nqty::Int: Number of quantities
• cos_coeffs::Array{Float64,3}: Cosine coefficients (npsi × nmodes × nqty)
• sin_coeffs::Array{Float64,3}: Sine coefficients (npsi × nmodes × nqty)

FourierCoefficients(xs, ys, fs, mband)

Compute Fourier coefficients via FFT without creating splines.

Arguments

• xs::Vector{Float64}: Radial coordinates
• ys::Vector{Float64}: Poloidal coordinates (periodic domain)
• fs::Array{Float64,3}: Function values (npsi × ntheta × nqty)
• mband::Int: Number of Fourier modes to retain

empty_FourierCoefficients()

Create an empty FourierCoefficients for initialization purposes.

get_complex_coeff(fc::FourierCoefficients, ipsi, mode, qty) -> ComplexF64

Get normalized complex FFT coefficient at grid point.

The FourierCoefficients internally stores:

• cos_coeffs = 2 * real(FFT) / ntheta (for m > 0)
• sin_coeffs = -2 * imag(FFT) / ntheta (for m > 0)

This function returns the normalized FFT coefficient (FFT/ntheta):

• c[0] = cos_coeffs[0] (DC component)
• c[m] = cos_coeffs[m]/2 - i*sin_coeffs[m]/2 (for m > 0)

get_complex_coeffs!(out, fc::FourierCoefficients, ipsi, qty)

Fill vector with normalized complex FFT coefficients for modes 0:mband.

FourierTransforms

Utility module for efficient Fourier transforms using pre-computed basis functions.

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θ + nqaδ), sin(mθ + nqaδ)) 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

Example

using GeneralizedPerturbedEquilibrium.Utilities.FourierTransforms

# For PerturbedEquilibrium (no phase shift)
ft = FourierTransform(mtheta, mpert, mlow)

# For Vacuum (with n*qa*delta phase)
ft_vac = FourierTransform(mthvac, mpert, mlow; n=1, qa=2.5, delta=delta_array)

# Forward transform: theta-space → Fourier modes
theta_data = randn(mtheta, 10)  # 10 different theta-space functions
modes = ft(theta_data)          # Complex modes [mpert, 10]

# Inverse transform: Fourier modes → theta-space
reconstructed = inverse(ft, modes)

FourierTransform

Callable struct for efficient Fourier transforms with pre-computed basis functions.

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:

# 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, inverse

(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:

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.

(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

Returns

• Complex Fourier coefficients
  • If input is Vector: returns Vector{ComplexF64} of length mpert
  • If input is Matrix: returns Matrix{ComplexF64} of size (mpert, n)

Formula

For each mode l (corresponding to mode number m = mlow + l - 1):

mode[l] = Σᵢ data[i] * (cos(m*θᵢ + phase) + im*sin(m*θᵢ + phase))

This is equivalent to computing:

real_part[l] = Σᵢ data[i] * cos(m*θᵢ + phase)
imag_part[l] = Σᵢ data[i] * sin(m*θᵢ + phase)
mode[l] = complex(real_part[l], imag_part[l])

Examples

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)

FourierTransform(mtheta, mpert, mlow; n=0, qa=0.0, delta=zeros(mtheta))

Construct a FourierTransform object with pre-computed basis functions.

Arguments

• mtheta::Int: Number of poloidal grid points
• mpert::Int: Number of Fourier modes
• mlow::Int: Lowest mode number

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

Returns

FourierTransform object ready for forward and inverse transforms.

Examples

# 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)

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

• 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

Keyword 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)

Returns

• 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]

Notes

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*θᵢ)

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 4π^2 / size(cs, 1) is 2π / mtheta * 2π / nzeta (nzeta = 1 for 2D)
• 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

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)

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

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)

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)

Formula

For each theta point i:

data[i] = (2π/mtheta) * Σₗ modes[l] * (cos(m*θᵢ + phase) + im*sin(m*θᵢ + phase))

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)]

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

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))

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.

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

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

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:

Re{mode[l]} = Σᵢ (data_real[i]*cos - data_imag[i]*sin)
Im{mode[l]} = Σᵢ (data_real[i]*sin + data_imag[i]*cos)

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

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