diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index a86b0059..f4090503 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -470,60 +470,3 @@ end # Initialize function for OdeState with relevant parameters for array initialization OdeState(numpert_total::Int, numsteps_init::Int, numunorms_init::Int, msing::Int) = OdeState(; numpert_total, numsteps_init, numunorms_init, msing) - -# Below here are debug output structs used for benchmarking and unit testing - - - - -""" - VacuumBenchmarkInputs - -A struct to hold all inputs required for vacuum benchmarking between Fortran and Julia implementations. - -## Fields - - - `wv_block::Matrix{ComplexF64}` - Vacuum response matrix block - - `mpert::Int` - Number of poloidal modes - - `mtheta_eq::Int` - Number of poloidal grid points in input equilibrium (corresponds to `mtheta_eq` in VacuumInput) - - `mthvac::Int` - Number of poloidal grid points in vacuum calculations (corresponds to `mtheta` in VacuumInput) - - `complex_flag::Bool` - Flag indicating if complex arithmetic is used - - `kernelsign::Float64` - Sign of the kernel for vacuum calculation - - `wall_flag::Bool` - Flag indicating presence of wall - - `farwall_flag::Bool` - Flag indicating presence of far wall - - `grri::Matrix{Float64}` - Green's function response matrix - - `xzpts::Matrix{Float64}` - Coordinate points on plasma boundary [R, Z] - - `ahg_file::String` - Filename for AHG data - - `dir_path::String` - Directory path for input/output files - - `vac_inputs::Vacuum.VacuumInput` - VacuumInput struct for Julia vacuum code - - `wall_settings::Vacuum.WallShapeSettings` - Wall shape settings - - `n::Int` - Toroidal mode number - - `ipert_n::Int` - Index of perturbed toroidal mode - - `psifac::Float64` - Normalized flux coordinate -""" -@kwdef struct VacuumBenchmarkInputs - # Vacuum computation parameters - wv_block::Matrix{ComplexF64} - mpert::Int - mtheta_eq::Int - mthvac::Int - complex_flag::Bool - kernelsign::Float64 - wall_flag::Bool - farwall_flag::Bool - grri::Matrix{Float64} - xzpts::Matrix{Float64} - ahg_file::String - dir_path::String - - # VacuumInput struct for Julia code - vac_inputs::Vacuum.VacuumInput - - # Wall settings - wall_settings::Vacuum.WallShapeSettings - - # Additional context - n::Int - ipert_n::Int - psifac::Float64 -end diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 9dd14a14..e0cfa5e1 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -28,7 +28,7 @@ function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibr for ipert_n in 1:intr.npert # Set VACUUM run parameters and boundary shape n = ipert_n - 1 + intr.nlow - vac_inputs = compute_vacuum_inputs(intr.psilim, n, equil, intr, ctrl) + vac_inputs = Vacuum.VacuumInput(equil, intr.psilim, ctrl.mthvac, intr.mpert, intr.mlow, n; force_wv_symmetry=ctrl.force_wv_symmetry) fill!(vac.grri, 0.0) fill!(vac.grre, 0.0) fill!(vac.xzpts, 0.0) @@ -36,20 +36,6 @@ function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibr # Compute vacuum energy matrix and both Green's functions wv_block, vac.grri, vac.grre, vac.xzpts = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) - # Output data for unit testing and benchmarking - if intr.debug_settings.output_benchmark_data - @info "Outputting top level vacuum debug data for n = $n" - farwall_flag = intr.wall_settings.shape == "nowall" ? true : false - benchmark_inputs = VacuumBenchmarkInputs( - wv_block, intr.mpert, equil.control.mtheta, ctrl.mthvac, - true, vac_inputs.kernelsign, false, - farwall_flag, vac.grri, vac.xzpts, "ahg2msc_dcon.out", intr.dir_path, - vac_inputs, intr.wall_settings, - n, ipert_n, intr.psilim - ) - @save "vacuum_response_inputs.jld2" benchmark_inputs - end - # Equation 126 in Chance 1997 - scale by (m - n*q)(m' - n*q) singfac = collect(intr.mlow:intr.mhigh) .- (n * intr.qlim) @inbounds for ipert in 1:intr.mpert @@ -129,66 +115,6 @@ function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibr return vac end -""" - compute_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl) - -Compute the necessary input paramters to the VACUUM code for computing the vacuum response matrix. -Performs the same function as `free_write_msc` in the Fortran code, except we create a data struct -to pass in memory to the Julia VACUUM code instead of writing to a file. Computed quantities include -the r, z, and ν values at the plasma boundary, as well as mode numbers and number of poloidal points. - -### Arguments - - - `psifac`: Flux surface value at the plasma boundary (Float64) - - `n`: Toroidal mode number (Int) - - `equil`: Plasma equilibrium data (Equilibrium.PlasmaEquilibrium) - - `intr`: Internal ForceFreeStates parameters - - `ctrl`: ForceFreeStates control parameters -""" -function compute_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl) - - # Allocations - theta_norm = equil.rzphi_ys - mtheta = equil.config.mtheta + 1 - angle = zeros(Float64, mtheta) - r = zeros(Float64, mtheta) - z = zeros(Float64, mtheta) - ν = zeros(Float64, mtheta) - rfac = zeros(Float64, mtheta) - - # Compute r, z, and ν at the plasma boundary - qa = equil.profiles.q_spline(psifac) - for itheta in 1:mtheta - # Evaluate geometric quantities at (ψ, θ) - r2 = equil.rzphi_rsquared((psifac, theta_norm[itheta])) - offset = equil.rzphi_offset((psifac, theta_norm[itheta])) - nu_val = equil.rzphi_nu((psifac, theta_norm[itheta])) - - rfac[itheta] = sqrt(r2) - angle[itheta] = 2π * (theta_norm[itheta] + offset) - ν[itheta] = nu_val - end - r .= equil.ro .+ rfac .* cos.(angle) - z .= equil.zo .+ rfac .* sin.(angle) - - # Invert values for n < 0 - if n < 0 - ν .= -ν - n = -n - end - - return Vacuum.VacuumInput(; - r=reverse(r), - z=reverse(z), - ν=reverse(ν), - mlow=intr.mlow, - mpert=intr.mpert, - n=n, - mtheta=ctrl.mthvac, - force_wv_symmetry=ctrl.force_wv_symmetry - ) -end - """ free_compute_wv_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) @@ -231,7 +157,7 @@ function free_compute_wv_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium for ipert_n in 1:intr.npert # Compute vacuum matrix n = ipert_n - 1 + intr.nlow - vac_inputs = compute_vacuum_inputs(intr.psilim, n, ctrl, equil, intr) + vac_inputs = Vacuum.VacuumInput(equil, intr.psilim, intr.mtheta, intr.mpert, intr.mlow, n; force_wv_symmetry=ctrl.force_wv_symmetry) wv_block, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) # Apply singular factor scaling diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 99d25c70..924b8ba2 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -28,28 +28,32 @@ Calculates the coupling between forcing modes and resonant surfaces, which deter the effectiveness of external perturbations at rational surfaces where q = m/n. Implements GPEC algorithm from gpout.f: -1. For each forcing mode, apply unit field and compute plasma response -2. For each singular surface, evaluate field jump and compute coupling quantities -3. Calculate diagnostic island widths and overlap parameters + + 1. For each forcing mode, apply unit field and compute plasma response + 2. For each singular surface, evaluate field jump and compute coupling quantities + 3. Calculate diagnostic island widths and overlap parameters ## Arguments -- `state`: Output state structure to store results -- `equil`: Equilibrium solution with q-profile and flux surfaces -- `ForceFreeStates_results`: ForceFreeStates stability calculation results -- `vac_data`: Vacuum response data including Green's functions -- `ffs_intr`: ForceFreeStates internal state with singular surface data -- `intr`: PerturbedEquilibrium internal state with permeability matrix -- `ctrl`: Control parameters + + - `state`: Output state structure to store results + - `equil`: Equilibrium solution with q-profile and flux surfaces + - `ForceFreeStates_results`: ForceFreeStates stability calculation results + - `vac_data`: Vacuum response data including Green's functions + - `ffs_intr`: ForceFreeStates internal state with singular surface data + - `intr`: PerturbedEquilibrium internal state with permeability matrix + - `ctrl`: Control parameters ## Modifies + Populates the following fields in `state`: -- `resonant_flux[msing, numpert_total]`: Normalized resonant flux Φ_r/A -- `resonant_current[msing, numpert_total]`: Resonant current density -- `island_width_sq[msing, numpert_total]`: Square of island half-width (w/2)² -- `penetrated_field[msing, numpert_total]`: Normal field at resonant surface -- `delta_prime[msing, numpert_total]`: Tearing stability parameter Δ' -- `island_half_width[msing]`: Dimensional island half-width -- `chirikov_parameter[msing]`: Island overlap metric + + - `resonant_flux[msing, numpert_total]`: Normalized resonant flux Φ_r/A + - `resonant_current[msing, numpert_total]`: Resonant current density + - `island_width_sq[msing, numpert_total]`: Square of island half-width (w/2)² + - `penetrated_field[msing, numpert_total]`: Normal field at resonant surface + - `delta_prime[msing, numpert_total]`: Tearing stability parameter Δ' + - `island_half_width[msing]`: Dimensional island half-width + - `chirikov_parameter[msing]`: Island overlap metric Note: numpert_total = mpert × npert handles all (m,n) mode combinations """ @@ -113,14 +117,13 @@ function compute_singular_coupling_metrics!( end # Get vacuum calculation parameters - mtheta_eq = length(equil.rzphi_ys) # Equilibrium poloidal grid size mtheta = vac_data.mthvac # Vacuum poloidal grid size mlow = ffs_intr.mlow nlow = ffs_intr.nlow nhigh = ffs_intr.nhigh # Perturbed equilibrium calculations always use nowall - wall_settings = Vacuum.WallShapeSettings(shape="nowall") + wall_settings = Vacuum.WallShapeSettings(; shape="nowall") if ctrl.verbose println(" Processing toroidal modes n = $nlow:$nhigh") @@ -156,16 +159,8 @@ function compute_singular_coupling_metrics!( end # Compute Green's functions at this surface for this n - vac_input = Vacuum.create_vacuum_input_at_psi( - equil, - sing_surf.psifac, - mtheta_eq, - mtheta, - mpert, - mlow, - nn - ) - grri, grre = Vacuum.compute_greens_functions_only(vac_input, wall_settings) + vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, mpert, mlow, nn) + _, grri, grre, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; green_only=true) # Store in singular surface struct (overwrites for each n) ffs_intr.sing[s].grri = grri @@ -311,7 +306,7 @@ Interpolate magnetic field at arbitrary flux surface. Interpolates displacement from ForceFreeStates solution and converts to field using ideal MHD relation from flux coordinates: - b^ψ = i * χ₁ * (m - n*q) * ξ_ψ +b^ψ = i * χ₁ * (m - n*q) * ξ_ψ This matches the field reconstruction in FieldReconstruction.jl. """ @@ -367,10 +362,10 @@ end Compute effective current density coefficient at given flux surface. Implements GPEC's j_c calculation (gpout.f line 560-578): - j_c = χ₁² * q / (μ₀ * integral) +j_c = χ₁² * q / (μ₀ * integral) where the integral is computed via flux surface integration: - integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ +integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ ## GPEC Formula @@ -394,9 +389,10 @@ j_c(ising)=1.0/j_c(ising)*chi1**2*sq%f(4)/mu0 Uses trapezoidal rule integration around the flux surface with metric quantities from the equilibrium bicubic spline (rzphi). The integrand includes: -- jac: Jacobian of flux coordinates -- |∇ψ|: Flux gradient magnitude (delpsi) -- sqreqb: Magnetic field quantity (F² + χ₁²|∇ψ|²)/(2πR)² + + - jac: Jacobian of flux coordinates + - |∇ψ|: Flux gradient magnitude (delpsi) + - sqreqb: Magnetic field quantity (F² + χ₁²|∇ψ|²)/(2πR)² """ function compute_current_density( equil::Equilibrium.PlasmaEquilibrium, @@ -436,8 +432,8 @@ function compute_current_density( r2 = equil.rzphi_rsquared((psi, theta)) # rfac² deta = equil.rzphi_offset((psi, theta)) # angle offset jac = equil.rzphi_jac((psi, theta)) # Jacobian - r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0,1)) # ∂(rfac²)/∂theta - deta_y = equil.rzphi_offset((psi, theta); deriv=(0,1)) # ∂(deta)/∂theta + r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1)) # ∂(rfac²)/∂theta + deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1)) # ∂(deta)/∂theta rfac = sqrt(abs(r2)) fy_rfac2 = r2_y @@ -493,11 +489,13 @@ Apply Green's function to Fourier mode coefficients. Simplified version that extracts only plasma surface rows (first mtheta rows). ## Arguments -- `green`: Green's function matrix [2*mtheta, 2*mpert] -- `mode_coeffs`: Complex Fourier coefficients [mpert] + + - `green`: Green's function matrix [2*mtheta, 2*mpert] + - `mode_coeffs`: Complex Fourier coefficients [mpert] ## Returns -- Potential at plasma surface theta points [mtheta] + + - Potential at plasma surface theta points [mtheta] """ function apply_green_function_simple( green::Matrix{Float64}, @@ -509,7 +507,7 @@ function apply_green_function_simple( # Pack complex to real/imag format packed = zeros(Float64, 2 * mpert) for i in 1:mpert - packed[2*i - 1] = real(mode_coeffs[i]) + packed[2*i-1] = real(mode_coeffs[i]) packed[2*i] = imag(mode_coeffs[i]) end @@ -530,21 +528,25 @@ end Compute surface inductance matrix from Green's functions at flux surface. This implements the full GPEC gpvacuum_flxsurf algorithm (gpvacuum.f line 299-331): -1. For each mode, apply Green's functions to unit flux -2. Compute surface current from potential jump: kax = (chi - che) / μ₀ -3. Fourier transform to mode space -4. Return inductance: L_surf = flux * inv(current) + + 1. For each mode, apply Green's functions to unit flux + 2. Compute surface current from potential jump: kax = (chi - che) / μ₀ + 3. Fourier transform to mode space + 4. Return inductance: L_surf = flux * inv(current) ## Arguments -- `grri`: Interior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) -- `grre`: Exterior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) -- `ffs_intr`: ForceFreeStates internal state -- `intr`: PerturbedEquilibrium internal state with mode arrays + + - `grri`: Interior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) + - `grre`: Exterior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) + - `ffs_intr`: ForceFreeStates internal state + - `intr`: PerturbedEquilibrium internal state with mode arrays ## Returns + Surface inductance matrix [numpert_total, numpert_total] ## GPEC Reference + ```fortran DO i=1,mpert vbwp_mn=0; vbwp_mn(i)=1.0 @@ -650,7 +652,7 @@ end Compute singular flux from resonant current using surface inductance. Uses surface inductance matrix at the singular surface: - Φ = L_surf * I +Φ = L_surf * I where L_surf is computed from Green's functions stored in sing_surf. @@ -718,7 +720,7 @@ end Compute flux surface area at given ψ. Implements GPEC's area calculation (gpout.f line 568): - area = ∫ jac * |∇ψ| dθ +area = ∫ jac * |∇ψ| dθ where the integral is computed around the flux surface. @@ -740,8 +742,9 @@ area(ising)=area(ising)-jac*delpsi(mthsurf)/mthsurf ! trapezoidal rule ## Implementation Uses trapezoidal rule integration around the flux surface with: -- jac: Jacobian of flux coordinates from rzphi -- |∇ψ|: Flux gradient magnitude (delpsi) from metric tensor + + - jac: Jacobian of flux coordinates from rzphi + - |∇ψ|: Flux gradient magnitude (delpsi) from metric tensor """ function compute_surface_area( equil::Equilibrium.PlasmaEquilibrium, @@ -772,8 +775,8 @@ function compute_surface_area( r2 = equil.rzphi_rsquared((psi, theta)) jac = equil.rzphi_jac((psi, theta)) deta = equil.rzphi_offset((psi, theta)) - r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0,1)) - deta_y = equil.rzphi_offset((psi, theta); deriv=(0,1)) + r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1)) + deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1)) # Compute rfac rfac = sqrt(abs(r2)) diff --git a/src/Utilities/FourierTransforms.jl b/src/Utilities/FourierTransforms.jl index 920261bd..4d3db17d 100644 --- a/src/Utilities/FourierTransforms.jl +++ b/src/Utilities/FourierTransforms.jl @@ -43,7 +43,7 @@ export transform!, inverse_transform! export fourier_transform!, fourier_inverse_transform! """ - compute_fourier_coefficients(mtheta, mpert, mlow; n=0, qa=0.0, delta=zeros(mtheta)) + compute_fourier_coefficients(mtheta, mpert, mlow; n=0, ν=zeros(mtheta)) Compute Fourier basis function coefficients for transforms between theta-space and mode-space. @@ -56,68 +56,46 @@ Compute Fourier basis function coefficients for transforms between theta-space a # Keyword Arguments - `n::Int=0`: Toroidal mode number (default: 0, no toroidal coupling) -- `qa::Float64=0.0`: Safety factor at boundary (default: 0.0) -- `delta::Vector{Float64}=zeros(mtheta)`: Toroidal phase shift array (default: no shift) +- `ν::Vector{Float64}=zeros(mtheta)`: Toroidal angle offset array (default: no offset, n*ν = 0) # Returns -- `cslth::Matrix{Float64}`: Cosine coefficients [mtheta, mpert] where `cslth[i,l] = cos(m*θᵢ + n*qa*δᵢ)` -- `snlth::Matrix{Float64}`: Sine coefficients [mtheta, mpert] where `snlth[i,l] = sin(m*θᵢ + n*qa*δᵢ)` +- `cos_mn_basis::Matrix{Float64}`: Cosine coefficients [mtheta, mpert] where `cos_mn_basis[i,l] = cos(m*θᵢ - n*νᵢ)` +- `sin_mn_basis::Matrix{Float64}`: Sine coefficients [mtheta, mpert] where `sin_mn_basis[i,l] = sin(m*θᵢ - n*νᵢ)` # Notes -The theta grid is uniform: `θᵢ = 2π*(i-1)/mtheta` for `i = 1:mtheta` +The theta grid is uniform: `θᵢ = 2π*i/mtheta` for `i = 0:mtheta-1` Mode numbers are: `m = mlow, mlow+1, ..., mlow+mpert-1` -When `n=0, qa=0, delta=0` (default), this reduces to simple harmonic basis: -- `cslth[i,l] = cos(m*θᵢ)` -- `snlth[i,l] = sin(m*θᵢ)` - -For vacuum calculations with toroidal coupling, the phase includes `n*qa*delta`: -- `cslth[i,l] = cos(m*θᵢ + n*qa*δᵢ)` -- `snlth[i,l] = sin(m*θᵢ + n*qa*δᵢ)` +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*θᵢ)` """ function compute_fourier_coefficients( mtheta::Int, mpert::Int, mlow::Int; n::Int=0, - qa::Float64=0.0, - delta::Vector{Float64}=zeros(Float64, mtheta) + ν::Vector{Float64}=zeros(Float64, mtheta) ) # Validate inputs @assert mtheta > 0 "mtheta must be positive" @assert mpert >= 0 "mpert must be non-negative" - @assert length(delta) == mtheta "delta must have length mtheta" + @assert length(ν) == mtheta "ν must have length mtheta" # Handle edge case: mpert = 0 returns empty arrays - if mpert == 0 - return zeros(Float64, mtheta, 0), zeros(Float64, mtheta, 0) - end + mpert == 0 && return zeros(Float64, mtheta, 0), zeros(Float64, mtheta, 0) # Uniform theta grid: [0, 2π) - theta_grid = range(0, 2π, length=mtheta+1)[1:end-1] - - # Compute toroidal phase shift (zero if n=0 or qa=0) - nqdelta = n * qa .* delta - - # Allocate coefficient matrices - cslth = zeros(Float64, mtheta, mpert) - snlth = zeros(Float64, mtheta, mpert) - - # Compute basis functions for each mode and theta point - for l in 1:mpert - m = mlow + l - 1 # Mode number - for i in 1:mtheta - # Total argument: m*θ + n*qa*δ - arg = theta_grid[i] * m + nqdelta[i] - cslth[i, l] = cos(arg) - snlth[i, l] = sin(arg) - end - end + θ_grid = range(; start=0, length=mtheta, step=2π/mtheta) + + # Compute sin(mθ - nν) and cos(mθ - nν) for all modes and theta points + sin_mn_basis = sin.((mlow .+ (0:(mpert-1))') .* θ_grid .- n .* ν) + cos_mn_basis = cos.((mlow .+ (0:(mpert-1))') .* θ_grid .- n .* ν) - return cslth, snlth + return cos_mn_basis, sin_mn_basis end """ @@ -186,7 +164,7 @@ Construct a FourierTransform object with pre-computed basis functions. ft = FourierTransform(480, 40, -20) # With toroidal phase (for Vacuum calculations) -ft_vac = FourierTransform(480, 40, -20; n=1, qa=2.5, delta=delta_array) +ft_vac = FourierTransform(480, 40, -20; n=1, ν=ν_array) ``` """ function FourierTransform( @@ -194,11 +172,10 @@ function FourierTransform( mpert::Int, mlow::Int; n::Int=0, - qa::Float64=0.0, - delta::Vector{Float64}=zeros(Float64, mtheta) + ν::Vector{Float64}=zeros(Float64, mtheta) ) - cslth, snlth = compute_fourier_coefficients(mtheta, mpert, mlow; n, qa, delta) - return FourierTransform(mtheta, mpert, mlow, cslth, snlth) + cos_mn_basis, sin_mn_basis = compute_fourier_coefficients(mtheta, mpert, mlow; n, ν) + return FourierTransform(mtheta, mpert, mlow, cos_mn_basis, sin_mn_basis) end """ diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index cfcde7bf..bc6a1c39 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -9,17 +9,15 @@ Struct holding plasma boundary and mode data as provided from ForceFreeStates na # Fields -- `r::Vector{Float64}`: Plasma boundary R-coordinate as a function of poloidal angle -- `z::Vector{Float64}`: Plasma boundary Z-coordinate as a function of poloidal angle -- `ν::Vector{Float64}`: Free parameter in specifying toroidal angle, ϕ = 2πζ + ν(ψ, θ), on theta grid -- `mlow::Int`: Lower poloidal mode number for spectral representation -- `mhigh::Int`: Upper poloidal mode number (mhigh = mlow + mpert - 1) -- `mpert::Int`: Number of poloidal modes (mhigh - mlow + 1) -- `n::Int`: Toroidal mode number -- `qa::Float64`: Safety factor at plasma boundary -- `mtheta_eq::Int`: Number of equilibrium poloidal grid points (input grid resolution) -- `mtheta::Int`: Number of vacuum calculation poloidal grid points -- `force_wv_symmetry::Bool`: Boolean flag to enforce symmetry in the vacuum response matrix + - `r::Vector{Float64}`: Plasma boundary R-coordinate as a function of poloidal angle + - `z::Vector{Float64}`: Plasma boundary Z-coordinate as a function of poloidal angle + - `ν::Vector{Float64}`: Free parameter in specifying toroidal angle, ϕ = 2πζ + ν(ψ, θ), on theta grid + - `mlow::Int`: Lower poloidal mode number for spectral representation + - `mhigh::Int`: Upper poloidal mode number (mhigh = mlow + mpert - 1) + - `mpert::Int`: Number of poloidal modes (mhigh - mlow + 1) + - `n::Int`: Toroidal mode number + - `mtheta::Int`: Number of vacuum calculation poloidal grid points + - `force_wv_symmetry::Bool`: Boolean flag to enforce symmetry in the vacuum response matrix """ @kwdef struct VacuumInput r::Vector{Float64} = Float64[] @@ -29,11 +27,76 @@ Struct holding plasma boundary and mode data as provided from ForceFreeStates na mhigh::Int = 0 mpert::Int = 0 n::Int = 0 - qa::Float64 = 1.0 - mtheta_eq::Int = 1 mtheta::Int = 1 force_wv_symmetry::Bool = true - # NOTE: kernelsign parameter deprecated - compute_vacuum_response now computes both grri and grre +end + +""" + VacuumInput( + equil::Equilibrium.PlasmaEquilibrium, + ψ::Float64, + mtheta::Int, + mpert::Int, + mlow::Int, + n::Int, + force_wv_symmetry::Bool = true + ) -> VacuumInput + +Constructor to create a VacuumInput struct for computing Green's functions at arbitrary flux surface. +Extracts plasma geometry from equilibrium at the given flux surface and packages it into VacuumInput format. + +## Arguments + + - `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) + +## Returns + +VacuumInput structure ready for compute_vacuum_response() + +## Usage + +This is used to compute Green's functions at singular surfaces: + +```julia +vac_input = VacuumInput(equil, sing_surf.psifac, mtheta, mpert, mlow, n; force_wv_symmetry=true) +``` +""" +function VacuumInput( + equil::Equilibrium.PlasmaEquilibrium, + ψ::Float64, + mtheta::Int, + mpert::Int, + mlow::Int, + n::Int; + force_wv_symmetry::Bool=true +) + # Extract plasma surface geometry at this psi + r, z, ν = extract_plasma_surface_at_psi(equil, ψ) + + # TODO: this was in Free.jl - is this general for GPEC too? + # Invert values for n < 0 + if n < 0 + ν .= -ν + n = -n + end + + return VacuumInput(; + r=reverse(r), + z=reverse(z), + ν=reverse(ν), + mlow=mlow, + mhigh=mlow + mpert - 1, + mpert=mpert, + n=n, + mtheta=mtheta, + force_wv_symmetry=force_wv_symmetry + ) end """ @@ -46,16 +109,12 @@ of length `mtheta`, where `mtheta` is the number of poloidal grid points and θ - `x::Vector{Float64}`: Plasma surface R-coordinate on VACUUM theta grid - `z::Vector{Float64}`: Plasma surface Z-coordinate on VACUUM theta grid - - `delta::Vector{Float64}`: Toroidal angle offset δ = -ν/(n*qa) for vacuum phase factor - - `dx_dtheta::Vector{Float64}`: Derivative dR/dθ at plasma surface - - `dz_dtheta::Vector{Float64}`: Derivative dZ/dθ at plasma surface + - `ν::Vector{Float64}`: Magnetic toroidal angle offset from geometric toroidal angle """ struct PlasmaGeometry x::Vector{Float64} z::Vector{Float64} - delta::Vector{Float64} - dx_dtheta::Vector{Float64} - dz_dtheta::Vector{Float64} + ν::Vector{Float64} end """ @@ -66,14 +125,14 @@ Struct holding wall geometry data for vacuum calculations. Arrays are of length # Fields - - `nowall::Bool`: Boolean flag indicating if there is no wall - - `x::Vector{Float64}`: Wall R-coordinates - - `z::Vector{Float64}`: Wall Z-coordinates + - `nowall::Bool`: Boolean flag indicating if there is no wall + - `x::Vector{Float64}`: Wall R-coordinates + - `z::Vector{Float64}`: Wall Z-coordinates """ -@kwdef struct WallGeometry - nowall::Bool = true - x::Vector{Float64} = Float64[] - z::Vector{Float64} = Float64[] +struct WallGeometry + nowall::Bool + x::Vector{Float64} + z::Vector{Float64} end """ @@ -83,29 +142,29 @@ Struct containing input settings for vacuum wall geometry. # Fields -- `shape::String`: String selecting wall shape. Options are: - - + `"nowall"`: No wall - + `"conformal"`: Wall conformal to plasma surface at distance `a` - + `"elliptical"`: Elliptical wall - + `"dee"`: Dee-shaped wall - + `"mod_dee"`: Modified Dee-shaped wall - + `"filepath"`: Custom wall shape from the file you specify - -- `a::Float64`: Distance of wall from plasma in units of major radius (conformal), or shape parameter (others) -- `aw::Float64`: Half-thickness parameter for Dee-shaped walls -- `bw::Float64`: Elongation parameter for wall shapes -- `cw::Float64`: Offset of the center of the wall from the major radius -- `dw::Float64`: Triangularity parameter for wall shapes -- `tw::Float64`: Sharpness of the corners of the wall (try 0.05 as initial value) -- `equal_arc_wall::Bool`: Flag to enforce equal arc length distribution of nodes on the wall - (recommended unless wall is very close to plasma) + - `shape::String`: String selecting wall shape. Options are: + + + `"nowall"`: No wall + + `"conformal"`: Wall conformal to plasma surface at distance `a` + + `"elliptical"`: Elliptical wall + + `"dee"`: Dee-shaped wall + + `"mod_dee"`: Modified Dee-shaped wall + + `"filepath"`: Custom wall shape from the file you specify + + - `a::Float64`: Distance of wall from plasma in units of major radius (conformal), or shape parameter (others) + - `aw::Float64`: Half-thickness parameter for Dee-shaped walls + - `bw::Float64`: Elongation parameter for wall shapes + - `cw::Float64`: Offset of the center of the wall from the major radius + - `dw::Float64`: Triangularity parameter for wall shapes + - `tw::Float64`: Sharpness of the corners of the wall (try 0.05 as initial value) + - `equal_arc_wall::Bool`: Flag to enforce equal arc length distribution of nodes on the wall + (recommended unless wall is very close to plasma) """ -@kwdef struct WallShapeSettings +@kwdef struct WallShapeSettings # Core shape selection shape::String = "nowall" - + # Standard geometric parameters for Dee/Mod-Dee a::Float64 = 0.3 aw::Float64 = 0.05 @@ -113,7 +172,7 @@ Struct containing input settings for vacuum wall geometry. cw::Float64 = 0.0 dw::Float64 = 0.5 tw::Float64 = 0.05 - + # Algorithmic options equal_arc_wall::Bool = true end @@ -121,93 +180,76 @@ end """ initialize_plasma_surface(inputs::VacuumInput) -> PlasmaGeometry -Initialize the plasma surface geometry based on the provided vacuum inputs. - -This function performs functionality from `readahg`, `arrays`, and `funint` in the -original Fortran VACUUM code. It returns a `PlasmaGeometry` struct containing -the necessary plasma surface data for vacuum calculations. - -# Process - -1. Interpolate the input plasma boundary arrays onto the mtheta grid -2. Compute derivatives of the plasma boundary with respect to poloidal angle θ - using periodic cubic spline differentiation -3. Compute trigonometric basis functions needed for Fourier calculations +Initialize the plasma surface geometry based on the provided vacuum inputs. +We interpolate the input plasma boundary arrays from the inputs struct onto the mtheta grid. # Arguments -- `inputs::VacuumInput`: Struct containing plasma boundary data and calculation parameters + - `inputs::VacuumInput`: Struct containing plasma boundary data # Returns -- `PlasmaGeometry`: Struct containing plasma surface coordinates, derivatives, and basis functions + - `PlasmaGeometry`: Struct containing plasma surface coordinates, derivatives, and basis functions """ -function initialize_plasma_surface(inputs::VacuumInput) +function PlasmaGeometry(inputs::VacuumInput) # Interpolate arrays from input onto mtheta grid - mtheta = inputs.mtheta - x_plasma = interp_to_new_grid(inputs.r, mtheta) - z_plasma = interp_to_new_grid(inputs.z, mtheta) - ν = interp_to_new_grid(inputs.ν, mtheta) - - # Compute delta from ν for vacuum phase factor - # delta = -ν/qa for use in phase: cos(m*θ + n*qa*δ) = cos(m*θ - n*ν) - delta = -ν ./ inputs.qa - - # Plasma boundary theta derivative: length mth with θ = [0, 1) - θ_grid = range(; start=0, length=mtheta, step=2π/mtheta) - dx_dtheta = periodic_cubic_deriv(θ_grid, x_plasma) - dz_dtheta = periodic_cubic_deriv(θ_grid, z_plasma) - - return PlasmaGeometry( - x_plasma, - z_plasma, - delta, - dx_dtheta, - dz_dtheta - ) + θ_in = range(0.0, 2π; length=length(inputs.r)) # VacuumInput uses [0, 2π] grid + θ_out = range(; start=0, length=inputs.mtheta, step=2π/inputs.mtheta) # VACUUM uses [0, 2π) grid + x = cubic_interp(θ_in, inputs.r; bc=PeriodicBC()).(θ_out) # no endpoint handling needed! + z = cubic_interp(θ_in, inputs.z; bc=PeriodicBC()).(θ_out) + ν = cubic_interp(θ_in, inputs.ν; bc=PeriodicBC()).(θ_out) + + return PlasmaGeometry(x, z, ν) end """ - initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_settings::WallShapeSettings) -> WallGeometry + WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_settings::WallShapeSettings) -> WallGeometry -Initialize the wall geometry based on the provided vacuum inputs and wall shape settings. +Constructor to initialize the wall geometry based on the provided vacuum inputs and wall shape settings. -This performs functionality similar to portions of the `arrays` function in the original -Fortran VACUUM code. It returns a `WallGeometry` struct containing the necessary wall +This performs functionality similar to portions of the `arrays` function in the original +Fortran VACUUM code. It returns a `WallGeometry` struct containing the necessary wall surface data for vacuum calculations. # Arguments -- `inputs::VacuumInput`: Struct containing vacuum calculation parameters -- `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry (used for reference) -- `wall_settings::WallShapeSettings`: Struct specifying wall shape and parameters + - `inputs::VacuumInput`: Struct containing vacuum calculation parameters + - `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry (used for reference) + - `wall_settings::WallShapeSettings`: Struct specifying wall shape and parameters # Returns -- `WallGeometry`: Struct containing wall surface coordinates and derivatives + - `WallGeometry`: Struct containing wall surface coordinates and derivatives # Notes -- Supports multiple wall shapes: nowall, conformal, elliptical, dee, mod_dee, from_file -- Optionally redistributes wall points to equal arc length spacing if `equal_arc_wall=true` + - Supports multiple wall shapes: nowall, conformal, elliptical, dee, mod_dee, from_file + - Optionally redistributes wall points to equal arc length spacing if `equal_arc_wall=true` """ -function initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_settings::WallShapeSettings) - - # Basic wall flags - nowall = wall_settings.shape == "nowall" +function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_settings::WallShapeSettings) - # All of these arrays are of length mtheta with θ = [0, 1) + # Output wall coordinate arrays mtheta = inputs.mtheta - - # Get wall shape from form_wall - # Plasma surface coordinates + x_wall = zeros(mtheta) + z_wall = zeros(mtheta) + + if wall_settings.shape == "nowall" + @info "Using no wall" + return WallGeometry( + true, + x_wall, + z_wall + ) + end + + # Compute plasma surface quantities x_plasma = plasma_surf.x z_plasma = plasma_surf.z # Output wall coordinate arrays x_wall = zeros(Float64, mtheta) - z_wall = zeros(Float64, mtheta) + z_wall = zeros(Float64, mtheta) # Common geometric parameters xmin = minimum(x_plasma) @@ -217,16 +259,14 @@ function initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_ r_minor = 0.5 * (xmax - xmin) r_major = 0.5 * (xmax + xmin) - + # Destructuring settings for readability (; aw, bw, cw, dw, tw, a) = wall_settings wcentr = 0.0 # Initialize - if wall_settings.shape == "nowall" - @info "Using no wall" - elseif wall_settings.shape == "conformal" + if wall_settings.shape == "conformal" dx = a * r_minor - @info "Calculating conformal wall shape $((@sprintf "%.2e" dx)) m from plasma surface." + @info "Calculating conformal wall shape $((@sprintf "%.2e" dx)) m from plasma surface." wcentr = r_major centerstack_min = min(0.1, 0.1 * minimum(x_plasma)) # Avoid wall crossing R=0 axis for i in 1:mtheta @@ -248,9 +288,9 @@ function initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_ zrad = 0.5 * (zmax - zmin) zh = sqrt(abs(zrad^2 - r_minor^2)) - zmuw = log((a/zh) + sqrt((a/zh)^2 + 1)) - bw_eff = (zh * cosh(zmuw)) / a - + zmuw = log((a/zh) + sqrt((a/zh)^2 + 1)) + bw_eff = (zh * cosh(zmuw)) / a + for i in 1:mtheta the = (i - 1) * (2π / mtheta) x_wall[i] = r_major + a * cos(the) @@ -277,8 +317,13 @@ function initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_ else filepath = wall_settings.shape - !isfile(filepath) && @error "ERROR: Wall geometry file $filepath does not exist. - Please set the wall shape parameter to a valid file path or a built-in shape (nowall, conformal, elliptical, dee, mod_dee)." + if !isfile(filepath) + error( + "Wall geometry file $filepath does not exist. " * + "Please set the wall shape parameter to a valid file path or a built-in shape " * + "(nowall, conformal, elliptical, dee, mod_dee)." + ) + end wcentr = 0.0 open(wall_settings.shape, "r") do io @@ -286,7 +331,7 @@ function initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_ wcentr = parse(Float64, readline(io)) readline(io) # Skip header/comment line - (npots0 < mtheta) && @error "ERROR: $filename contains fewer points ($npots0) than mtheta ($mtheta)." + npots0 < mtheta && error("Wall geometry file $filepath contains fewer points ($npots0) than mtheta ($mtheta).") for i in 1:mtheta line = split(readline(io)) @@ -297,20 +342,14 @@ function initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_ end end - # Optional: Re-parameterization - if wall_settings.equal_arc_wall && (wall_settings.shape != "nowall") + # Optional: Re-parameterization for equal arc length spacing of wall points + if wall_settings.equal_arc_wall @info "Re-distributing wall points to equal arc length spacing (assumes closed, toroidal wall)." - x_wall, z_wall, _, _, _ = distribute_to_equal_arc_grid(x_wall, z_wall, mtheta) + x_wall, z_wall = distribute_to_equal_arc_grid(x_wall, z_wall) end - if any(x_wall .<= 0.0) && !nowall - # to add support for x<0 walls, be sure to carefully replicate Chance's fortran code x<0 handling in the kernel function to account for the additional singularities associated with this - error("Wall R-coordinates contain non-physical values (R <= 0). Check wall geometry.") - end + # To add support for x<0 walls, be sure to carefully replicate Chance's fortran code x<0 handling in the kernel function to account for the additional singularities associated with this + any(x_wall .<= 0.0) && error("Wall R-coordinates contain non-physical values (R <= 0). Check wall geometry.") - return WallGeometry( - nowall=nowall, - x=x_wall, - z=z_wall - ) + return WallGeometry(false, x_wall, z_wall) end diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index bce1b7a7..1da6bc1d 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/src/Vacuum/Kernel2D.jl @@ -45,29 +45,25 @@ const GAUSSIANPOINTS32 = [ ] """ - kernel!(grad_greenfunction_mat, greenfunction_mat, x_obspoints, z_obspoints, x_sourcepoints, z_sourcepoints, j1, j2, isgn, iops, inputs; xwall=nothing, zwall=nothing) + kernel!(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. The residue calculation needs to be updated for open walls.** # Arguments - - `grad_greenfunction_mat`: Gradient Green's function matrix (output) - - `greenfunction_mat`: Green's function matrix (output) - - `x_obspoints`: Observer x coordinates (R coordinates) - - `z_obspoints`: Observer z coordinates (Z coordinates) - - `x_sourcepoints`: Source x coordinates (R coordinates) - - `z_sourcepoints`: Source z coordinates (Z coordinates) - - `j1/j2`: Block index for observer/source (1=plasma, 2=wall) + - `grad_greenfunction`: Gradient Green's function matrix (output) + - `greenfunction`: Green's function matrix (output) + - `observer`: Observer geometry struct (PlasmaGeometry or WallGeometry) + - `source`: Source geometry struct (PlasmaGeometry or WallGeometry) - `n`: Toroidal mode number # Returns -Modifies `grad_greenfunction_mat` and `greenfunction_mat` in place. -Note that greenfunction_mat is zeroed each time this function is called, -but grad_greenfunction_mat is not since it fills a different block of the +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 (2 * mtheta, 2 * mtheta) depending on the source/observer. # Notes @@ -77,87 +73,72 @@ but grad_greenfunction_mat is not since it fills a different block of the - Implements analytical singularity removal following Chance 1997 """ function kernel!( - grad_greenfunction_mat::Matrix{Float64}, - greenfunction_mat::Matrix{Float64}, - x_obspoints::Vector{Float64}, - z_obspoints::Vector{Float64}, - x_sourcepoints::Vector{Float64}, - z_sourcepoints::Vector{Float64}, - j1::Int, - j2::Int, + grad_greenfunction::Matrix{Float64}, + greenfunction::Matrix{Float64}, + observer::Union{PlasmaGeometry,WallGeometry}, + source::Union{PlasmaGeometry,WallGeometry}, n::Int ) - # These used to be function arguments, but can just set inside here based on j1/j2 - plasma_plasma_block = j1 == 1 && j2 == 1 # previously iops - plasma_is_source = j2 == 1 # previously iopw - isgn = plasma_is_source ? -1 : 1 - - mtheta = length(x_obspoints) + mtheta = length(observer.x) dtheta = 2π / mtheta theta_grid = range(; start=0, length=mtheta, step=dtheta) - # Zero out greenfunction_mat at start of each kernel call (matches Fortran behavior) - fill!(greenfunction_mat, 0.0) + # Take a view of the corresponding block of the grad_greenfunction + col_index = (source isa PlasmaGeometry ? 1 : 2) + row_index = (observer isa PlasmaGeometry ? 1 : 2) + grad_greenfunction_block = view( + grad_greenfunction, + ((row_index-1)*mtheta+1):(row_index*mtheta), + ((col_index-1)*mtheta+1):(col_index*mtheta) + ) - if mtheta != length(z_obspoints) || mtheta != length(x_sourcepoints) || mtheta != length(z_sourcepoints) - error("Length of input arrays (xobs, zobs, xsource, zsce) are different. All length should be the same") - end + # 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 # S₁ᵢ in Chance 1997, eq.(78) log_correction_0=16.0*dtheta*(log(2*dtheta)-68.0/15.0)/15.0 log_correction_1=128.0*dtheta*(log(2*dtheta)-8.0/15.0)/45.0 log_correction_2=4.0*dtheta*(7.0*log(2*dtheta)-11.0/15.0)/45.0 - # Used for Z'_θ and X'_θ in eq.(51) - # Close the loop for periodic BC by appending first point at the end - theta_closed = vcat(collect(theta_grid), theta_grid[end] + dtheta) - x_closed = vcat(x_sourcepoints, x_sourcepoints[1]) - z_closed = vcat(z_sourcepoints, z_sourcepoints[1]) - spline_x = cubic_interp(theta_closed, x_closed; bc=PeriodicBC()) - spline_z = cubic_interp(theta_closed, z_closed; bc=PeriodicBC()) - # Create derivative views once, reuse for all evaluations (avoids allocation per call) + # Precompute composite Simpson's 1/3 rule weights, excluding singular points + # Note we set to 4 for even/2 for odd since we index from 1 while the formula assumes indexing from 0 + nsrc = mtheta - 3 + simpson_weights = dtheta / 3 .* [(k == 1 || k == nsrc) ? 1 : (iseven(k) ? 4 : 2) for k in 1:nsrc] + + # Set up periodic splines used for off-grid Gaussian quadrature points + spline_x = cubic_interp(theta_grid, source.x; bc=PeriodicBC(; endpoint=:exclusive)) + spline_z = cubic_interp(theta_grid, source.z; bc=PeriodicBC(; endpoint=:exclusive)) d1_spline_x = deriv1(spline_x) d1_spline_z = deriv1(spline_z) - theta_vec = collect(theta_grid) - dx_dtheta = d1_spline_x.(theta_vec) - dz_dtheta = d1_spline_z.(theta_vec) # Loop through observer points for j in 1:mtheta - # Initialize variables - x_obs=x_obspoints[j] - z_obs=z_obspoints[j] - theta_obs=theta_grid[j] - grad_green_0 = 0.0 # simpson integral for coupling_0 (𝒥 ∇'𝒢⁰∇'ℒ) - # Workspace = view of appropriate row of grad_greenfunction_mat for this observer point - grad_green_work = @view(grad_greenfunction_mat[(j1-1)*mtheta+j, (j2-1)*mtheta .+ (1:mtheta)]) + # Get observer coordinates + x_obs, z_obs, theta_obs = observer.x[j], observer.z[j], theta_grid[j] # Obtain nonsingular region (endpoints at j+2 and j-2, so exclude j-1, j, and j+1) - nonsing_src_indices = mod1.((j+2):(j+mtheta-2), mtheta) # mod1 ensures isrc is in [1, mtheta] - - # Compute composite Simpson's 1/3 rule weights (https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_1/3_rule) - # Note we set to 4 for even/2 for odd since we index from 1 while the formula assumes indexing from 0 - nsrc = length(nonsing_src_indices) - simpson_weights = dtheta / 3 .* [(k == 1 || k == nsrc) ? 1 : (iseven(k) ? 4 : 2) for k in 1:nsrc] + nonsing_idx = mod1.((j+2):(j+mtheta-2), mtheta) # mod1 ensures isrc is in [1, mtheta] # Perform Simpson integration for nonsingular source points - for (isrc, wsimpson) in zip(nonsing_src_indices, simpson_weights) - x_source=x_sourcepoints[isrc] - z_source=z_sourcepoints[isrc] - - # G_n is 2pi𝒢ⁿ; coupling_n is 𝒥 ∇'𝒢ⁿ∇'ℒ; coupling_0 is 𝒥 ∇'𝒢ⁿ∇'ℒ for n=0 - G_n, coupling_n, coupling_0 = green(x_obs, z_obs, x_source, z_source, dx_dtheta[isrc], dz_dtheta[isrc], n) + for (isrc, wsimpson) in zip(nonsing_idx, simpson_weights) + dx_dtheta, dz_dtheta = d1_spline_x(theta_grid[isrc]), d1_spline_z(theta_grid[isrc]) + G_n, gradG_n, gradG_0 = green(x_obs, z_obs, source.x[isrc], source.z[isrc], dx_dtheta, dz_dtheta, n) # Sum contributions to Green's function matrices using Simpson weight - grad_green_work[isrc] += isgn * coupling_n * wsimpson - greenfunction_mat[j, isrc] += G_n * wsimpson - grad_green_0 += coupling_0 * wsimpson + if populate_greenfunction + greenfunction[j, isrc] += G_n * wsimpson + end + grad_greenfunction_block[j, isrc] += gradG_n * wsimpson + # Subtract regular integral component of δⱼᵢK⁰ in eq. 83 + grad_greenfunction_block[j, j] -= gradG_0 * wsimpson end # Perform Gaussian quadrature for singular points (source = obs point) # Get indices of the singularity region, [j-2, j-1, j, j+1, j+2] - js = mod.(j .+ ((mtheta-3):(mtheta+1)), mtheta) .+ 1 + sing_idx = mod1.(j .+ ((mtheta-2):(mtheta+2)), mtheta) # Integrate region of length 2 * dtheta on left/right of singularity for region in ["left", "right"] gauss_xleft = theta_obs - (region == "left" ? 2 * dtheta : 0) @@ -171,10 +152,7 @@ function kernel!( dx_dtheta_gauss = d1_spline_x(theta_gauss0) z_gauss = spline_z(theta_gauss0) dz_dtheta_gauss = d1_spline_z(theta_gauss0) - G_n, coupling_n, coupling_0 = green(x_obs, z_obs, x_gauss, z_gauss, dx_dtheta_gauss, dz_dtheta_gauss, n) - - # Add logarithm to G_n to analytically isolate the singularity (first type), Chance eq.(75) - G_n_nonsingular = plasma_plasma_block ? G_n + log((theta_obs-theta_gauss[ig])^2)/x_obs : G_n + G_n, gradG_n, gradG_0 = green(x_obs, z_obs, x_gauss, z_gauss, dx_dtheta_gauss, dz_dtheta_gauss, n) # Redefine hardcoded Gaussian weights on the interval [-1, 1] to physical interval with length 2 * dtheta wgauss = GAUSSIANWEIGHTS[ig] * dtheta @@ -188,49 +166,61 @@ function kernel!( A2_minus = (pgauss^2-1)*pgauss*(pgauss-2)/24.0 * wgauss # First type of singularity: 𝒢ⁿ, occurs plasma as source only (see RHS of Chance eqs. 26/27) - if plasma_is_source - greenfunction_mat[j, js[1]] += G_n_nonsingular * A2_minus - greenfunction_mat[j, js[2]] += G_n_nonsingular * A1_minus - greenfunction_mat[j, js[3]] += G_n_nonsingular * A0 - greenfunction_mat[j, js[4]] += G_n_nonsingular * A1_plus - greenfunction_mat[j, js[5]] += G_n_nonsingular * A2_plus + if populate_greenfunction + if observer isa PlasmaGeometry + # Remove singular behavior by adding on leading-order term, Chance eq.(75) + G_n += log((theta_obs - theta_gauss[ig])^2) / x_obs + end + greenfunction[j, sing_idx[1]] += G_n * A2_minus + greenfunction[j, sing_idx[2]] += G_n * A1_minus + greenfunction[j, sing_idx[3]] += G_n * A0 + greenfunction[j, sing_idx[4]] += G_n * A1_plus + greenfunction[j, sing_idx[5]] += G_n * A2_plus end - # Second type of singularity: 𝒦ⁿ - # Eq. 86: 𝒦ⁿαᵢ - δⱼᵢK⁰ - grad_green_work[js[1]] += isgn * coupling_n * A2_minus - grad_green_work[js[2]] += isgn * coupling_n * A1_minus - grad_green_work[js[3]] += isgn * coupling_n * A0 - grad_green_work[js[4]] += isgn * coupling_n * A1_plus - grad_green_work[js[5]] += isgn * coupling_n * A2_plus + # Second type of singularity: 𝒦ⁿ (Eq. 86: 𝒦ⁿαᵢ - δⱼᵢK⁰) + grad_greenfunction_block[j, sing_idx[1]] += gradG_n * A2_minus + grad_greenfunction_block[j, sing_idx[2]] += gradG_n * A1_minus + grad_greenfunction_block[j, sing_idx[3]] += gradG_n * A0 + grad_greenfunction_block[j, sing_idx[4]] += gradG_n * A1_plus + grad_greenfunction_block[j, sing_idx[5]] += gradG_n * A2_plus # Subtract off the diverging singular n=0 component - grad_green_work[j] -= isgn * coupling_0 * wgauss + grad_greenfunction_block[j, j] -= gradG_0 * wgauss end end - # Set residue based on logic similar to Table I of Chance 1997 + existing δⱼᵢ in eq. 69 - # Would need to pass in wall geometry to generalize this to open walls - is_closed_toroidal = true - if is_closed_toroidal - residue = (j1 == 2.0) ? 0.0 : (j2 == 1 ? 2.0 : -2.0) # Chance eq. 89 - else - # TODO: this line can be gotten rid of if we are never doing open walls - residue = (j1 == j2) ? 2.0 : 0.0 # Chance eq. 90 - end - # Subtract regular integral component of δⱼᵢK⁰ in eq. 83 and add residue value in eq. 89/90 - grad_green_work[j] = grad_green_work[j] - isgn * grad_green_0 + residue - # Subtract off analytic singular integral from Chance eq.(75) if plasma-plasma block - if plasma_plasma_block - greenfunction_mat[j, js[1]] -= log_correction_2 / x_obs - greenfunction_mat[j, js[2]] -= log_correction_1 / x_obs - greenfunction_mat[j, js[3]] -= log_correction_0 / x_obs - greenfunction_mat[j, js[4]] -= log_correction_1 / x_obs - greenfunction_mat[j, js[5]] -= log_correction_2 / x_obs + if populate_greenfunction && observer isa PlasmaGeometry + greenfunction[j, sing_idx[1]] -= log_correction_2 / x_obs + greenfunction[j, sing_idx[2]] -= log_correction_1 / x_obs + greenfunction[j, sing_idx[3]] -= log_correction_0 / x_obs + greenfunction[j, sing_idx[4]] -= log_correction_1 / x_obs + greenfunction[j, sing_idx[5]] -= log_correction_2 / x_obs end end + + # Normals need to point outward from vacuum region. In VACUUM clockwise θ convention, normal points + # out of vacuum for wall but inward for plasma, so we multiply by -1 for wall sources + if source isa PlasmaGeometry + grad_greenfunction_block .*= -1 + end + + # Add analytic singular integral (second type) from Table I of Chance 1997 + existing δⱼᵢ in eq. 69 + # Would need to pass in wall geometry to generalize this to open walls + is_closed_toroidal = true + if is_closed_toroidal # Chance eq. 89 + residue = (observer isa WallGeometry) ? 0.0 : (source isa PlasmaGeometry ? 2.0 : -2.0) + else # Chance eq. 90 + # TODO: this line can be gotten rid of if we are never doing open walls + residue = (typeof(observer) == typeof(source)) ? 2.0 : 0.0 + end + # Add residue value from eq. 89/90 to block diagonal + @inbounds for i in 1:mtheta + grad_greenfunction_block[i, i] += residue + end + # Since we computed 2π𝒢, divide by 2π to get 𝒢 - greenfunction_mat ./= 2π + greenfunction ./= 2π end ############################################################# diff --git a/src/Vacuum/MathUtils.jl b/src/Vacuum/MathUtils.jl index 0df35dae..40fc0689 100644 --- a/src/Vacuum/MathUtils.jl +++ b/src/Vacuum/MathUtils.jl @@ -1,242 +1,49 @@ -# NOTE: fourier_transform! and fourier_inverse_transform! have been moved to -# src/Utilities/FourierTransforms.jl and are imported at the module level in Vacuum.jl - -# Returns the array of derivatives at all x points, I think this acts like difspl -# in the Fortran but need to check/consolidate spline routines later -function periodic_cubic_deriv(theta, vals) - # Close the loop for periodic BC by appending first point at the end - theta_closed = vcat(collect(theta), theta[end] + (theta[2] - theta[1])) - vals_closed = vcat(vals, vals[1]) - spline = cubic_interp(theta_closed, vals_closed; bc=PeriodicBC()) - d1 = deriv1(spline) - return d1.(collect(theta)) -end - """ - lagrange1d(ax, af, m, nl, x, iop) + distribute_to_equal_arc_grid(xin, zin) -Perform Lagrange interpolation and optionally compute its derivative. -Replaces `lagp`, `lagpe4`, and `lag` from Fortran code. +Given a set of points (xin, zin) that define a closed curve in 2D, redistribute these points to be equally spaced +in terms of arc length along the curve. This is useful for ensuring that points on a wall or plasma surface are +uniformly distributed in space rather than in the parameter θ. This performs the same function as eqarcw in the +Fortran code. We now use FastInterpolations instead of a manual Lagrange interpolation. # Arguments - - `ax::Vector{Float64}`: Array of x-coordinates for the interpolation points - - `af::Vector{Float64}`: Array of y-coordinates (function values) for the interpolation points - - `m::Int`: Number of interpolation points - - `nl::Int`: Number of points to use for the local interpolation (polynomial degree + 1) - - `x::Float64`: The x-value at which to evaluate the interpolated function and/or its derivative - - `iop::Int`: Flag controlling output (0 = value only, 1 = value and derivative) + - `xin::Vector{Float64}`: x-coordinates of the original points defining the curve (endpoint not included) + - `zin::Vector{Float64}`: z-coordinates of the original points defining the curve (endpoint not included) # Returns - - `f::Float64`: The interpolated function value at `x` - - `df::Float64`: The interpolated function derivative at `x` (0.0 if iop=0) - -# Notes - - - Uses local Lagrange interpolation with `nl` points centered around `x` - - Automatically adjusts interpolation window to stay within array bounds -""" -function lagrange1d(ax::Vector{Float64}, af::Vector{Float64}, m::Int, nl::Int, x::Float64, iop::Int) - - # --- Error fix: Initialize f and df internally to 0.0 --- - f::Float64 = 0.0 - df::Float64 = 0.0 - - jn = findfirst(i -> ax[i] >= x, 1:m) - jn = (jn === nothing) ? m : jn - jn = max(jn - 1, 1) - if jn < m && abs(ax[jn+1] - x) < abs(x - ax[jn]) - jn += 1 - end - - # Determine the range of indices for interpolation - jnmm = floor(Int, (nl - 0.1) / 2) - jnpp = floor(Int, (nl + 0.1) / 2) - - nll = jn - jnmm - nlr = jn + jnpp - - # Adjust for even nl when ax[jn] > x (Window shift) - if (nl % 2 == 0) && (ax[jn] > x) - nll -= 1 - nlr -= 1 # <--- Shifting the window (was nlr += 1) - end - - # Clamp indices to valid array bounds - if nlr > m - nlr = m - nll = nlr - nl + 1 - elseif nll < 1 - nll = 1 - nlr = nl - end - - # Compute function value f - for i in nll:nlr - alag = 1.0 - for j in nll:nlr - (i == j) && continue - alag *= (x - ax[j]) / (ax[i] - ax[j]) - end - f += alag * af[i] - end - - # --- Error fix: Use the 'iop' argument --- - (iop == 0) && return f, df # df is returned as 0.0 - - # Compute derivative df - for i in nll:nlr - slag = 0.0 - for id in nll:nlr - (id == i) && continue - alag = 1.0 - for j in nll:nlr - (j == i) && continue - alag *= (j != id) ? ((x - ax[j]) / (ax[i] - ax[j])) : (1.0 / (ax[i] - ax[id])) - end - slag += alag - end - df += slag * af[i] - end - - return f, df # Return value and derivative -end - + - `xout::Vector{Float64}`: x-coordinates of the redistributed points, equally spaced in arc length + - `zout::Vector{Float64}`: z-coordinates of the redistributed points, equally spaced in arc length """ - interp_to_new_grid(vecin, mtheta; dx0=0.0, dx1=0.0) +function distribute_to_equal_arc_grid(xin::Vector{Float64}, zin::Vector{Float64}) -Resample the input array `vecin` using a periodic cubic spline to an output array of length `mtheta`. + mtheta = length(xin) + dθ = 2π / mtheta + θ_grid = range(; start=0, length=mtheta, step=dθ) -This function unifies the Fortran functions `trans`, `transdx`, and `transdxx` into a single -function with optional offset parameters. - -# Arguments - - - `vecin::Vector{Float64}`: Input array to be resampled - - `mtheta::Int`: Desired length of the output array - - `dx0::Float64`: Global offset added to all x-coordinates (default 0, applied as `x += dx0 / mtheta_in`) - - `dx1::Float64`: Fine offset added to each index (default 0, applied as `ai = (i-1) + dx1`) - -# Returns - - - `vecout::Vector{Float64}`: The resampled output array (length `mtheta`) - -# Notes - - - If `mtheta == length(vecin)`, returns the input vector unchanged - - Uses periodic cubic spline interpolation for resampling - - Input grid is normalized to [0, 1] for interpolation -""" -function interp_to_new_grid(vecin::Vector{Float64}, mtheta::Int; dx0=0.0, dx1=0.0) - - # Initialize - mtheta_in = length(vecin) - - # If mtheta == mtheta_in, just return the input vector - if mtheta == mtheta_in - return vecin - end - - # Enforce periodicity: FastInterpolations requires y[1] ≈ y[end] for PeriodicBC - # For truly periodic data (poloidal angle), force endpoints to match exactly - vecin_periodic = copy(vecin) - vecin_periodic[end] = vecin_periodic[1] - - # Input grids are from [0, 1] inclusive, while no interpolants will fall outside of this, the bc still impacts the end behavior - θin = collect(range(0.0, 1.0; length=mtheta_in)) - spline = cubic_interp(θin, vecin_periodic; bc=PeriodicBC()) - - # Interpolate to new grid with optional offsets - vecout = zeros(mtheta) - for i in 1:mtheta - x = (i - 1 + dx1) / mtheta + dx0 / mtheta_in - x = x % 1.0 # This is for periodicity in the case of dx1/dx0 ≠ 0 - vecout[i] = spline(x) - end - return vecout -end - -""" - distribute_to_equal_arc_grid(xin, zin, mw1) - -Perform arc length re-parameterization of a 2D curve. - -Takes an input curve defined by `(xin, zin)` coordinates and re-samples it such that -the new points `(xout, zout)` are equally spaced in arc length along the curve. - -# Arguments - - - `xin::Vector{Float64}`: Array of x-coordinates of the input curve - - `zin::Vector{Float64}`: Array of z-coordinates of the input curve - - `mw1::Int`: Number of points in the input and output curves - -# Returns - - - `xout::Vector{Float64}`: Array of x-coordinates of the arc-length re-parameterized curve - - `zout::Vector{Float64}`: Array of z-coordinates of the arc-length re-parameterized curve - - `ell::Vector{Float64}`: Array of cumulative arc lengths for the input curve - - `thgr::Vector{Float64}`: Array of re-parameterized 'theta' values corresponding to equal arc lengths - - `thlag::Vector{Float64}`: Array of normalized 'theta' values for the input curve (0 to 1) - -# Notes - - - Uses Lagrange interpolation for calculating arc length and resampling - - Ensures uniform spacing in arc length for improved numerical stability -""" -function distribute_to_equal_arc_grid(xin::Vector{Float64}, zin::Vector{Float64}, mtheta::Int) - - # Temporary arrays for interpolation and arc-length calculation - theta_in = zeros(Float64, mtheta) # Normalized input parameter [0, 1) - theta_out = zeros(Float64, mtheta) # New parameter distribution for equal spacing - xout = zeros(Float64, mtheta) # Uniformly spaced R-coordinates - zout = zeros(Float64, mtheta) # Uniformly spaced Z-coordinates - - # Define initial normalized parameter theta_in - dt = 1.0 / mtheta - theta_in .= range(; start=0, length=mtheta, step=dt) # θ ∈ [0, 1) - # we need a closed loop for arc length calculation - mtheta1 = mtheta + 1 - xin1 = vcat(xin, xin[1]) - zin1 = vcat(zin, zin[1]) - theta_in1 = vcat(theta_in, [1.0]) - ell = zeros(Float64, mtheta1) # Cumulative arc length of closed loop + # Build periodic splines for derivatives on the closed loop + spline_x = cubic_interp(θ_grid, xin; bc=PeriodicBC(; endpoint=:exclusive)) + spline_z = cubic_interp(θ_grid, zin; bc=PeriodicBC(; endpoint=:exclusive)) # Calculate cumulative arc length using numerical integration - # We use a mid-point derivative approximation to find the path length - for iw in 2:mtheta1 - # Evaluate derivative at the midpoint of the interval - theta = (theta_in1[iw] + theta_in1[iw-1]) / 2.0 - - # Calculate dx/dt and dz/dt using Lagrange interpolation (order 3) - _, d_xin = lagrange1d(theta_in1, xin1, mtheta1, 3, theta, 1) - _, d_zin = lagrange1d(theta_in1, zin1, mtheta1, 3, theta, 1) - - # Instantaneous speed (ds/dt) - ds_dt = sqrt(d_xin^2 + d_zin^2) - - # Accumulate length: ds = (ds/dt) * dt - ell[iw] = ell[iw-1] + ds_dt * dt - end - - # Re-parameterize based on equal arc-length segments - ell_targets = collect(range(0; step=ell[end]/mtheta, length=mtheta)) # [0, Length) for open loop result - for i in 2:mtheta - # Find the value of 'theta_in' that corresponds to the target arc length 's' - f_th, _ = lagrange1d(ell, theta_in1, mtheta1, 3, ell_targets[i], 0) - theta_out[i] = f_th - end - - # Interpolate the original (x,z) data at the new parameter points to get (xout, zout) - # Chance interpolates in theta_out to get xin, zin but this introduces small errors in arc lengths + arc_length = zeros(mtheta + 1) for i in 1:mtheta - f_x, _ = lagrange1d(ell, xin1, mtheta1, 3, ell_targets[i], 0) - f_z, _ = lagrange1d(ell, zin1, mtheta1, 3, ell_targets[i], 0) - xout[i] = f_x - zout[i] = f_z + # Use a mid-point derivative approximation + theta_mid = θ_grid[i] + dθ / 2.0 + d_xin = spline_x(theta_mid; deriv=1) + d_zin = spline_z(theta_mid; deriv=1) + # Accumulate length: ds = (ds/dt) * dt + ds_dθ = sqrt(d_xin^2 + d_zin^2) + arc_length[i+1] = arc_length[i] + ds_dθ * dθ end - return xout, zout, ell, theta_out, theta_in + # Re-parameterize based on equal arc-length segments by interpolate the original (x,z) data at the equal arc length points + arc_length_targets = range(; start=0, length=mtheta, step=arc_length[end]/mtheta) + # arc_length is an uneven grid so we have to specify the period explicitly + xout = cubic_interp(arc_length[1:(end-1)], xin; bc=PeriodicBC(; endpoint=:exclusive, period=arc_length[end])).(arc_length_targets) + zout = cubic_interp(arc_length[1:(end-1)], zin; bc=PeriodicBC(; endpoint=:exclusive, period=arc_length[end])).(arc_length_targets) + return xout, zout end # Helper functions for compute_vacuum_field diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index a2696d88..1b848b53 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -22,7 +22,7 @@ export mscvac, set_surface_params, VacuumInput, compute_vacuum_response export compute_vacuum_field export kernel! export WallShapeSettings -export extract_plasma_surface_at_psi, create_vacuum_input_at_psi, compute_greens_functions_only +export extract_plasma_surface_at_psi, create_vacuum_input_at_psi # ====================================================================== # Legacy fortran vacuum module interface @@ -257,61 +257,61 @@ It computes both interior (grri) and exterior (grre) Green's functions for GPEC - The vacuum response includes plasma-plasma and plasma-wall coupling effects - For n=0 modes with closed walls, a regularization factor is added to prevent singularities """ -function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; green_only=false) # Initialization and allocations - (; mtheta, mpert, mlow, n, qa, force_wv_symmetry) = inputs - plasma_surf = initialize_plasma_surface(inputs) - wall = initialize_wall(inputs, plasma_surf, wall_settings) + (; mtheta, mpert, mlow, n, force_wv_symmetry) = inputs + plasma_surf = PlasmaGeometry(inputs) + wall = WallGeometry(inputs, plasma_surf, wall_settings) # Compute Fourier basis coefficients using FourierTransforms utility # We only need the coefficient arrays for the existing fourier_transform! functions - cos_ln_basis, sin_ln_basis = compute_fourier_coefficients(mtheta, mpert, mlow; n=n, qa=qa, delta=plasma_surf.delta) + cos_mn_basis, sin_mn_basis = compute_fourier_coefficients(mtheta, mpert, mlow; n=n, ν=plasma_surf.ν) # Allocate arrays for both Green's functions grri = zeros(2 * mtheta, 2 * mpert) # Interior (kernelsign=-1) grre = zeros(2 * mtheta, 2 * mpert) # Exterior (kernelsign=+1) - grad_greenfunction_mat = zeros(2 * mtheta, 2 * mtheta) - greenfunction_temp = zeros(mtheta, mtheta) + grad_green = zeros(2 * mtheta, 2 * mtheta) + green_temp = zeros(mtheta, mtheta) + + # Fourier transforms offsets into grri/grre: first mtheta rows are plasma as observer, second are wall + # First mpert columns are real (cosine), second mpert are imaginary (sine) + PLASMA_ROW_OFFSET = 0 + WALL_ROW_OFFSET = mtheta + COS_COL_OFFSET = 0 + SIN_COL_OFFSET = mpert # Plasma–Plasma block - j1, j2 = 1, 1 - kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, plasma_surf.x, plasma_surf.z, j1, j2, n) + kernel!(grad_green, green_temp, plasma_surf, plasma_surf, n) - # Fourier transform plasma-plasma block - # Populate both grri and grre with the same right-hand side - fourier_transform!(grri, greenfunction_temp, cos_ln_basis, 0, 0) - fourier_transform!(grri, greenfunction_temp, sin_ln_basis, 0, mpert) - fourier_transform!(grre, greenfunction_temp, cos_ln_basis, 0, 0) - fourier_transform!(grre, greenfunction_temp, sin_ln_basis, 0, mpert) + # Fourier transform obs=plasma, src=plasma block + fourier_transform!(grri, green_temp, cos_mn_basis, PLASMA_ROW_OFFSET, COS_COL_OFFSET) + fourier_transform!(grri, green_temp, sin_mn_basis, PLASMA_ROW_OFFSET, SIN_COL_OFFSET) + fourier_transform!(grre, green_temp, cos_mn_basis, PLASMA_ROW_OFFSET, COS_COL_OFFSET) + fourier_transform!(grre, green_temp, sin_mn_basis, PLASMA_ROW_OFFSET, SIN_COL_OFFSET) - !wall.nowall && begin + if !wall.nowall # Plasma–Wall block - j1, j2 = 1, 2 - kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, wall.x, wall.z, j1, j2, n) - + kernel!(grad_green, green_temp, plasma_surf, wall, n) # Wall–Wall block - j1, j2 = 2, 2 - kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, wall.x, wall.z, j1, j2, n) - + kernel!(grad_green, green_temp, wall, wall, n) # Wall–Plasma block - j1, j2 = 2, 1 - kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, plasma_surf.x, plasma_surf.z, j1, j2, n) - - # Fourier transform wall blocks into both grri and grre - fourier_transform!(grri, greenfunction_temp, cos_ln_basis, mtheta, 0) - fourier_transform!(grri, greenfunction_temp, sin_ln_basis, mtheta, mpert) - fourier_transform!(grre, greenfunction_temp, cos_ln_basis, mtheta, 0) - fourier_transform!(grre, greenfunction_temp, sin_ln_basis, mtheta, mpert) + kernel!(grad_green, green_temp, wall, plasma_surf, n) + + # Fourier transform obs=wall, src=plasma block + fourier_transform!(grri, green_temp, cos_mn_basis, WALL_ROW_OFFSET, COS_COL_OFFSET) + fourier_transform!(grri, green_temp, sin_mn_basis, WALL_ROW_OFFSET, SIN_COL_OFFSET) + fourier_transform!(grre, green_temp, cos_mn_basis, WALL_ROW_OFFSET, COS_COL_OFFSET) + fourier_transform!(grre, green_temp, sin_mn_basis, WALL_ROW_OFFSET, SIN_COL_OFFSET) end # Add cn0 to make grdgre nonsingular for n=0 modes cn0 = 1.0 # expose to user if anyone ever actually tries to use this - (abs(n) <= 1e-5 && !wall.nowall && wall.is_closed_toroidal) && begin + (n == 0 && !wall.nowall && wall.is_closed_toroidal) && begin @warn "Adding $cn0 to diagonal of grdgre to regularize n=0 mode; this may affect accuracy of results." mth12 = wall.nowall ? mtheta : 2 * mtheta for i in 1:mth12, j in 1:mth12 - grad_greenfunction_mat[i, j] += cn0 + grad_green[i, j] += cn0 end end @@ -320,47 +320,47 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe # grre: exterior potential (kernelsign=+1) # Make copies for each kernelsign - grad_greenfunction_mat_interior = copy(grad_greenfunction_mat) - grad_greenfunction_mat_exterior = copy(grad_greenfunction_mat) + grad_green_interior = copy(grad_green) + grad_green_exterior = copy(grad_green) # Apply kernelsign transformations - apply_kernelsign!(grad_greenfunction_mat_interior, -1.0, mtheta) # Interior - apply_kernelsign!(grad_greenfunction_mat_exterior, +1.0, mtheta) # Exterior (no-op) + apply_kernelsign!(grad_green_interior, -1.0, mtheta) # Interior + apply_kernelsign!(grad_green_exterior, +1.0, mtheta) # Exterior (no-op) # Invert the vacuum response system for both cases # If plasma only, lower blocks will be empty if wall.nowall - @views grri[1:mtheta, :] .= grad_greenfunction_mat_interior[1:mtheta, 1:mtheta] \ grri[1:mtheta, :] - @views grre[1:mtheta, :] .= grad_greenfunction_mat_exterior[1:mtheta, 1:mtheta] \ grre[1:mtheta, :] + @views grri[1:mtheta, :] .= grad_green_interior[1:mtheta, 1:mtheta] \ grri[1:mtheta, :] + @views grre[1:mtheta, :] .= grad_green_exterior[1:mtheta, 1:mtheta] \ grre[1:mtheta, :] else - grri .= grad_greenfunction_mat_interior \ grri - grre .= grad_greenfunction_mat_exterior \ grre + grri .= grad_green_interior \ grri + grre .= grad_green_exterior \ grre end # There's some logic that computes xpass/zpass and chiwc/chiws here, might eventually be needed? - # Perform inverse Fourier transforms to get response matrix components (eq. 115-118 of Chance 2007) - arr = zeros(mpert, mpert) - aii = zeros(mpert, mpert) - ari = zeros(mpert, mpert) - air = zeros(mpert, mpert) - fourier_inverse_transform!(arr, grre, cos_ln_basis, 0, 0) - fourier_inverse_transform!(aii, grre, sin_ln_basis, 0, mpert) - fourier_inverse_transform!(ari, grre, sin_ln_basis, 0, 0) - fourier_inverse_transform!(air, grre, cos_ln_basis, 0, mpert) - - # Final form of vacuum response matrix (eq. 114 of Chance 2007) - wv = complex.(arr .+ aii, air .- ari) - - # Force symmetry of response matrix if desired - force_wv_symmetry && hermitianpart!(wv) - - # Create xzpts array + wv = zeros(mpert, mpert) xzpts = zeros(Float64, inputs.mtheta, 4) - @views xzpts[:, 1] .= plasma_surf.x - @views xzpts[:, 2] .= plasma_surf.z - @views xzpts[:, 3] .= wall.x - @views xzpts[:, 4] .= wall.z + if !green_only + # Perform inverse Fourier transforms to get response matrix components (eq. 115-118 of Chance 2007) + arr, aii, ari, air = ntuple(_ -> zeros(mpert, mpert), 4) + fourier_inverse_transform!(arr, grre, cos_mn_basis, PLASMA_ROW_OFFSET, COS_COL_OFFSET) + fourier_inverse_transform!(aii, grre, sin_mn_basis, PLASMA_ROW_OFFSET, SIN_COL_OFFSET) + fourier_inverse_transform!(ari, grre, sin_mn_basis, PLASMA_ROW_OFFSET, COS_COL_OFFSET) + fourier_inverse_transform!(air, grre, cos_mn_basis, PLASMA_ROW_OFFSET, SIN_COL_OFFSET) + + # Final form of vacuum response matrix (eq. 114 of Chance 2007) + wv = complex.(arr .+ aii, air .- ari) + + # Force symmetry of response matrix if desired + force_wv_symmetry && hermitianpart!(wv) + + # Create xzpts array + @views xzpts[:, 1] .= plasma_surf.x + @views xzpts[:, 2] .= plasma_surf.z + @views xzpts[:, 3] .= wall.x + @views xzpts[:, 4] .= wall.z + end return wv, grri, grre, xzpts end diff --git a/src/Vacuum/VacuumFromEquilibrium.jl b/src/Vacuum/VacuumFromEquilibrium.jl index 296621a6..0b0481dd 100644 --- a/src/Vacuum/VacuumFromEquilibrium.jl +++ b/src/Vacuum/VacuumFromEquilibrium.jl @@ -1,272 +1,62 @@ """ -Functions to extract vacuum calculation inputs from equilibrium at arbitrary flux surfaces. - -This allows computing Green's functions at interior flux surfaces (e.g., singular surfaces) -by extracting the plasma geometry from the equilibrium bicubic spline. -""" - -# Note: Equilibrium module is imported in parent scope via ..Equilibrium -# Import FourierTransforms for coefficient calculation -using ..Utilities.FourierTransforms: compute_fourier_coefficients, fourier_transform! - -""" - extract_plasma_surface_at_psi( - equil::Equilibrium.PlasmaEquilibrium, - psi::Float64, - mtheta_eq::Int - ) -> (r, z, delta, qa) + extract_plasma_surface_at_psi(equil::Equilibrium.PlasmaEquilibrium, ψ::Float64) -> (r, z, ν) Extract plasma surface geometry from equilibrium at specified flux coordinate. Evaluates equilibrium bicubic spline around the flux surface to get R, Z coordinates -and computes the toroidal angle offset delta for vacuum calculations. +and computes the toroidal angle offset ν for vacuum calculations. ## Arguments -- `equil`: Equilibrium solution with rzphi bicubic spline -- `psi`: Normalized flux coordinate (0 at axis, 1 at edge) -- `mtheta_eq`: Number of poloidal grid points for surface extraction + + - `equil`: Equilibrium solution with rzphi bicubic spline + - `ψ`: Normalized flux coordinate (0 at axis, 1 at edge) ## Returns + Tuple of: -- `r::Vector{Float64}`: R-coordinates around flux surface [mtheta_eq] -- `z::Vector{Float64}`: Z-coordinates around flux surface [mtheta_eq] -- `delta::Vector{Float64}`: Toroidal angle offset δ = -ν/qa [mtheta_eq] -- `qa::Float64`: Safety factor at this surface + + - `r::Vector{Float64}`: R-coordinates around flux surface [mtheta] + - `z::Vector{Float64}`: Z-coordinates around flux surface [mtheta] + - `ν::Vector{Float64}`: Toroidal angle offset ν [mtheta] ## Implementation The equilibrium bicubic spline `rzphi` stores: -- rzphi.f[1] = r² (or rfac²) -- rzphi.f[2] = deta (angle offset) -- rzphi.f[3] = dphi (toroidal angle) -- rzphi.f[4] = jac (Jacobian) + + - rzphi.f[1] = r² (or rfac²) + - rzphi.f[2] = offset (straight-fieldline poloidal angle offset from geometric poloidal angle) + - rzphi.f[3] = ν (toroidal angle offset from geometric toroidal angle) + - rzphi.f[4] = jac (Jacobian) From these we compute: -- rfac = √(rzphi.f[1]) -- eta = 2π*(θ + rzphi.f[2]) -- R = R₀ + rfac*cos(eta) -- Z = Z₀ + rfac*sin(eta) -- delta = rzphi.f[3] / qa -## Reference + - r_minor = √(rzphi.f[1]) + - θ_cyl = 2π*(θ_sfl + rzphi.f[2]) + - R = R₀ + r_minor * cos(θ_cyl) + - Z = Z₀ + r_minor * sin(θ_cyl) + +## Reference # Allocate output arrays + Matches GPEC's ahg_write and gpvacuum_flxsurf approach (gpvacuum.f line 291-296) """ -function extract_plasma_surface_at_psi( - equil::Equilibrium.PlasmaEquilibrium, - psi::Float64, - mtheta_eq::Int -) - # Get magnetic axis location - ro = equil.ro - zo = equil.zo - - # Get safety factor at this surface - qa = equil.profiles.q_spline(psi) +function extract_plasma_surface_at_psi(equil::Equilibrium.PlasmaEquilibrium, ψ::Float64) # Allocate output arrays - r = zeros(Float64, mtheta_eq) - z = zeros(Float64, mtheta_eq) - delta = zeros(Float64, mtheta_eq) + mtheta = length(equil.rzphi_ys) + r_minor = zeros(mtheta) + θ_cyl = zeros(mtheta) + ν = zeros(mtheta) # Evaluate equilibrium around the flux surface - twopi = 2π - for itheta in 0:mtheta_eq-1 - # Theta coordinate normalized to [0, 1) - theta = itheta / mtheta_eq - - # Evaluate bicubic splines at (psi, theta) - # New API uses separate interpolants for each component - r2 = equil.rzphi_rsquared((psi, theta)) # r² or rfac² - deta = equil.rzphi_offset((psi, theta)) # angle offset - dphi = equil.rzphi_nu((psi, theta)) # toroidal angle offset (nu) - - rfac = sqrt(abs(r2)) - - # Compute R, Z coordinates - eta = twopi * (theta + deta) - r[itheta+1] = ro + rfac * cos(eta) - z[itheta+1] = zo + rfac * sin(eta) - - # Toroidal angle offset: delta = -ν/qa where ϕ = 2πζ + ν - # In GPEC: delta = dphi (already normalized) - delta[itheta+1] = dphi - end - - return r, z, delta, qa -end - -""" - create_vacuum_input_at_psi( - equil::Equilibrium.PlasmaEquilibrium, - psi::Float64, - mtheta_eq::Int, - mtheta::Int, - mpert::Int, - mlow::Int, - n::Int - ) -> VacuumInput - -Create VacuumInput structure for computing Green's functions at arbitrary flux surface. - -Extracts plasma geometry from equilibrium and packages it into VacuumInput format. - -## Arguments -- `equil`: Equilibrium solution -- `psi`: Normalized flux coordinate -- `mtheta_eq`: Number of equilibrium poloidal points -- `mtheta`: Number of vacuum calculation poloidal points -- `mpert`: Number of perturbing poloidal modes -- `mlow`: Lowest poloidal mode number -- `n`: Toroidal mode number - -## Returns -VacuumInput structure ready for compute_vacuum_response() - -## Usage -This is used to compute Green's functions at singular surfaces: -```julia -vac_input = create_vacuum_input_at_psi(equil, sing_surf.psifac, ...) -grri, grre = compute_greens_functions_only(vac_input, wall_settings) -``` -""" -function create_vacuum_input_at_psi( - equil::Equilibrium.PlasmaEquilibrium, - psi::Float64, - mtheta_eq::Int, - mtheta::Int, - mpert::Int, - mlow::Int, - n::Int -) - # Extract plasma surface geometry at this psi - r, z, delta, qa = extract_plasma_surface_at_psi(equil, psi, mtheta_eq) - - # Convert delta to ν for VacuumInput - # Relationship: delta = -ν/qa, so ν = -delta*qa - ν = -delta .* qa - - # Create VacuumInput structure - mhigh = mlow + mpert - 1 - - return VacuumInput( - r = r, - z = z, - ν = ν, - mlow = mlow, - mhigh = mhigh, - mpert = mpert, - n = n, - qa = qa, - mtheta_eq = mtheta_eq, - mtheta = mtheta, - force_wv_symmetry = true - ) -end - -""" - compute_greens_functions_only( - inputs::VacuumInput, - wall_settings::WallShapeSettings - ) -> (grri, grre) - -Compute only Green's functions without full vacuum response matrix. - -This is a lightweight version of compute_vacuum_response() that skips -the wv matrix calculation and only returns the Green's functions needed -for surface inductance calculations. - -## Arguments -- `inputs`: Vacuum input with plasma geometry -- `wall_settings`: Wall shape settings - -## Returns -Tuple of: -- `grri::Matrix{Float64}`: Interior Green's function [2*mtheta, 2*mpert] -- `grre::Matrix{Float64}`: Exterior Green's function [2*mtheta, 2*mpert] - -## Performance -Much faster than full compute_vacuum_response() since it skips: -- wv matrix construction -- Eigenvalue calculations -- Mode coupling matrices - -## Reference -Mimics GPEC's gpvacuum_flxsurf which only computes Green's functions -for surface inductance (gpvacuum.f line 279-283) -""" -function compute_greens_functions_only( - inputs::VacuumInput, - wall_settings::WallShapeSettings -) - # Reuse initialization from compute_vacuum_response - (; mtheta, mpert, mlow, n, qa) = inputs - plasma_surf = initialize_plasma_surface(inputs) - wall = initialize_wall(inputs, plasma_surf, wall_settings) - - # Compute Fourier basis coefficients - cslth, snlth = compute_fourier_coefficients(mtheta, mpert, mlow; n=n, qa=qa, delta=plasma_surf.delta) - - # Allocate arrays for both Green's functions - grri = zeros(2 * mtheta, 2 * mpert) - grre = zeros(2 * mtheta, 2 * mpert) - grad_greenfunction_mat = zeros(2 * mtheta, 2 * mtheta) - greenfunction_temp = zeros(mtheta, mtheta) - - # Plasma–Plasma block - j1, j2 = 1, 1 - kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, plasma_surf.x, plasma_surf.z, j1, j2, n) - - # Fourier transform plasma-plasma block - fourier_transform!(grri, greenfunction_temp, cslth, 0, 0) - fourier_transform!(grri, greenfunction_temp, snlth, 0, mpert) - fourier_transform!(grre, greenfunction_temp, cslth, 0, 0) - fourier_transform!(grre, greenfunction_temp, snlth, 0, mpert) - - !wall.nowall && begin - # Plasma–Wall block - j1, j2 = 1, 2 - kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, wall.x, wall.z, j1, j2, n) - - # Wall–Wall block - j1, j2 = 2, 2 - kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, wall.x, wall.z, j1, j2, n) - - # Wall–Plasma block - j1, j2 = 2, 1 - kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, plasma_surf.x, plasma_surf.z, j1, j2, n) - - # Fourier transform wall blocks - fourier_transform!(grri, greenfunction_temp, cslth, mtheta, 0) - fourier_transform!(grri, greenfunction_temp, snlth, mtheta, mpert) - fourier_transform!(grre, greenfunction_temp, cslth, mtheta, 0) - fourier_transform!(grre, greenfunction_temp, snlth, mtheta, mpert) - end - - # Add cn0 for n=0 modes if needed - cn0 = 1.0 - (abs(n) <= 1e-5 && !wall.nowall && wall.is_closed_toroidal) && begin - mth12 = wall.nowall ? mtheta : 2 * mtheta - for i in 1:mth12, j in 1:mth12 - grad_greenfunction_mat[i, j] += cn0 - end - end - - # Apply kernelsign transformations and invert to get Green's functions - # grri: interior (kernelsign=-1), grre: exterior (kernelsign=+1) - grad_greenfunction_mat_interior = copy(grad_greenfunction_mat) - grad_greenfunction_mat_exterior = copy(grad_greenfunction_mat) - - apply_kernelsign!(grad_greenfunction_mat_interior, -1.0, mtheta) # Interior - apply_kernelsign!(grad_greenfunction_mat_exterior, +1.0, mtheta) # Exterior - - # Invert the system for both cases - if wall.nowall - @views grri[1:mtheta, :] .= grad_greenfunction_mat_interior[1:mtheta, 1:mtheta] \ grri[1:mtheta, :] - @views grre[1:mtheta, :] .= grad_greenfunction_mat_exterior[1:mtheta, 1:mtheta] \ grre[1:mtheta, :] - else - grri .= grad_greenfunction_mat_interior \ grri - grre .= grad_greenfunction_mat_exterior \ grre + for (i, θ_sfl) in enumerate(equil.rzphi_ys) + # Get minor radius, geometric poloidal angle, and toroidal angle offset + r_minor[i] = sqrt(equil.rzphi_rsquared((ψ, θ_sfl))) + θ_cyl[i] = 2π * (θ_sfl + equil.rzphi_offset((ψ, θ_sfl))) + ν[i] = equil.rzphi_nu((ψ, θ_sfl)) end - return grri, grre + # Compute R and Z using the computed cylindrical coordinates + r = equil.ro .+ r_minor .* cos.(θ_cyl) + z = equil.zo .+ r_minor .* sin.(θ_cyl) + return r, z, ν end diff --git a/test/runtests_vacuum_julia.jl b/test/runtests_vacuum_julia.jl index b62abd01..427ab5e4 100644 --- a/test/runtests_vacuum_julia.jl +++ b/test/runtests_vacuum_julia.jl @@ -24,27 +24,23 @@ using LinearAlgebra @testset "initialize_plasma_surface" begin inputs = VacuumInput( - r=[1.0, 1.1, 1.2, 1.1], - z=[0.0, 0.1, 0.0, -0.1], - ν=zeros(4), - mtheta=4, + r=[1.0, 1.1, 1.2, 1.1, 1.0], + z=[0.0, 0.1, 0.0, -0.1, 0.0], + ν=zeros(5), + mtheta=5, mpert=1, mlow=1, n=1 ) - plasma_surf = JPEC.Vacuum.initialize_plasma_surface(inputs) - @test length(plasma_surf.x) == 4 - @test length(plasma_surf.z) == 4 - @test length(plasma_surf.delta) == 4 - @test length(plasma_surf.dx_dtheta) == 4 - @test length(plasma_surf.dz_dtheta) == 4 - @test !any(isnan, plasma_surf.dx_dtheta) - @test !any(isnan, plasma_surf.dz_dtheta) + plasma_surf = JPEC.Vacuum.PlasmaGeometry(inputs) + @test length(plasma_surf.x) == 5 + @test length(plasma_surf.z) == 5 + @test length(plasma_surf.ν) == 5 end @testset "initialize_wall" begin inputs = VacuumInput(mtheta=16) - plasma_surf = JPEC.Vacuum.initialize_plasma_surface( + plasma_surf = JPEC.Vacuum.PlasmaGeometry( VacuumInput( r=1.7 .+ 0.3 .* cos.(range(0, 2pi, length=16)), z=0.3 .* sin.(range(0, 2pi, length=16)), @@ -55,19 +51,19 @@ using LinearAlgebra # Test "nowall" wall_settings = WallShapeSettings(shape="nowall") - wall_geo = JPEC.Vacuum.initialize_wall(inputs, plasma_surf, wall_settings) + wall_geo = JPEC.Vacuum.WallGeometry(inputs, plasma_surf, wall_settings) @test wall_geo.nowall == true # Test "conformal" wall wall_settings = WallShapeSettings(shape="conformal", a=0.2) - wall_geo = JPEC.Vacuum.initialize_wall(inputs, plasma_surf, wall_settings) + wall_geo = JPEC.Vacuum.WallGeometry(inputs, plasma_surf, wall_settings) @test wall_geo.nowall == false @test length(wall_geo.x) == 16 @test !any(isnan, wall_geo.x) # Test that wall with R <= 0 throws an error # Create a plasma surface very close to R=0 - plasma_surf_near_zero = JPEC.Vacuum.initialize_plasma_surface( + plasma_surf_near_zero = JPEC.Vacuum.PlasmaGeometry( VacuumInput( r=0.05 .+ 0.03 .* cos.(range(0, 2pi, length=16)), z=0.03 .* sin.(range(0, 2pi, length=16)), @@ -78,13 +74,13 @@ using LinearAlgebra # Use a "dee" wall shape with parameters that will produce R < 0 # Setting cw (offset) to a large negative value will shift the wall left past R=0 wall_settings_negative = WallShapeSettings(shape="dee", cw=-1.5, a=0.1) - @test_throws ErrorException JPEC.Vacuum.initialize_wall(inputs, plasma_surf_near_zero, wall_settings_negative) + @test_throws ErrorException JPEC.Vacuum.WallGeometry(inputs, plasma_surf_near_zero, wall_settings_negative) # Test that conformal wall R-coordinates are clamped by centerstack_min # With a very large 'a' parameter, a conformal wall would naturally go to R < 0, # but it should be clamped to centerstack_min = min(0.1, 0.1 * minimum(x_plasma)) wall_settings_large_a = WallShapeSettings(shape="conformal", a=10.0, equal_arc_wall=false) - wall_geo_clamped = JPEC.Vacuum.initialize_wall(inputs, plasma_surf_near_zero, wall_settings_large_a) + wall_geo_clamped = JPEC.Vacuum.WallGeometry(inputs, plasma_surf_near_zero, wall_settings_large_a) expected_min = min(0.1, 0.1 * minimum(plasma_surf_near_zero.x)) @test all(wall_geo_clamped.x .>= expected_min) @test any(wall_geo_clamped.x .<= expected_min + 1e-10) # At least one point should be at the minimum @@ -95,7 +91,7 @@ using LinearAlgebra theta = range(0, step=2pi/10, length=10) xin_circ = cos.(theta) zin_circ = sin.(theta) - xout_circ, zout_circ, _, _, _ = JPEC.Vacuum.distribute_to_equal_arc_grid(xin_circ, zin_circ, 10) + xout_circ, zout_circ = JPEC.Vacuum.distribute_to_equal_arc_grid(xin_circ, zin_circ) # The points should still be on the unit circle @test all(r -> isapprox(r, 1.0, atol=1e-9), sqrt.(xout_circ .^ 2 + zout_circ .^ 2)) end @@ -160,47 +156,6 @@ using LinearAlgebra # ---------------------------------------------------------------------- @testset "VacuumInternals.jl - Interpolation/Grid" begin - @testset "lagrange1d" begin - ax = [1.0, 2.0, 3.0, 4.0] - af = [2.0, 4.0, 6.0, 8.0] # f(x) = 2x - m = 4 - nl = 3 - - # Test value - f, df = JPEC.Vacuum.lagrange1d(ax, af, m, nl, 2.5, 0) - @test isapprox(f, 5.0, atol=1e-12) - - # Test value and derivative - f, df = JPEC.Vacuum.lagrange1d(ax, af, m, nl, 2.5, 1) - @test isapprox(f, 5.0, atol=1e-12) - @test isapprox(df, 2.0, atol=1e-12) - end - - @testset "interp_to_new_grid" begin - # Test upsampling a sine wave - mtheta_in = 17 - theta_in = collect(range(0, 1, length=mtheta_in)) - vecin = sin.(2π .* theta_in) - - mtheta_out = 33 - vecout = JPEC.Vacuum.interp_to_new_grid(vecin, mtheta_out) - - theta_out = (0:(mtheta_out-1)) ./ mtheta_out - expected_out = sin.(2π .* theta_out) - - @test isapprox(vecout, expected_out, atol=1e-2) - - # Test no-op - vecout_noop = JPEC.Vacuum.interp_to_new_grid(vecin, mtheta_in) - @test vecout_noop == vecin - end - - @testset "periodic_cubic_deriv" begin - theta = range(0, 2pi, length=101)[1:(end-1)] - vals = sin.(theta) - derivs = JPEC.Vacuum.periodic_cubic_deriv(theta, vals) - @test isapprox(derivs, cos.(theta), atol=1e-3) - end @testset "_create_pickup_grid" begin R_grid = [1, 2]