diff --git a/docs/development/architecture.md b/docs/development/architecture.md index b9ca1a43f..c2db18e62 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -46,13 +46,15 @@ GPEC consists of **seven main modules** organized in `src/`: 4. **Vacuum** (`src/Vacuum/`) - Vacuum field calculations and Green's functions - Computes vacuum response matrices for ideal MHD analysis - - Calculates both **interior** (grri) and **exterior** (grre) Green's functions + - Solves the exterior boundary-integral system for the vacuum energy matrix `wv`, and optionally + (`compute_Iv=true`) the interior system as well to build the surface-current matrix `I_v` + (Park 2007 eq. 21b). The interior/exterior Green's functions themselves are internal scratch. - Main functions: - - `compute_vacuum_response()` - Pure Julia implementation + - `compute_vacuum_response()` / `compute_vacuum_response!()` - allocating and in-place entry points - Key files: - - `VacuumStructs.jl` - Data structures - - `VacuumInternals.jl` - Core algorithms - - `VacuumFromEquilibrium.jl` - Integration with equilibrium data + - `DataTypes.jl` - Data structures (`VacuumInput`, `PlasmaGeometry`, `WallGeometry`) + - `Kernel2D.jl` / `Kernel3D.jl` - Single-/double-layer kernel assembly + - `Field.jl` - Vacuum field and potential evaluation off the surface - Status: **Pure Julia implementation complete and available** 5. **ForceFreeStates** (`src/ForceFreeStates/`) - Ideal MHD stability analysis (DCON-style) @@ -130,8 +132,7 @@ The complete GPEC analysis pipeline: 2. **Vacuum Response**: - Initialize plasma and wall surfaces from equilibrium - - Compute vacuum response matrices (wv, grri, grre) - - Calculate both interior and exterior Green's functions + - Compute the vacuum energy matrix `wv` (and `I_v` when `compute_Iv=true`) - Pure Julia implementation 3. **Stability Analysis** (ForceFreeStates): @@ -169,11 +170,9 @@ The complete GPEC analysis pipeline: ### Stability - `SingType` - Singular surface data including: - - Rational surface location (ψ, q = m/n) - - Δ' (tearing stability parameter) - - Eigenmode structure at singular surface - - Green's functions (grri, grre) at interior singular surfaces - - Surface inductance + - Rational surface location (ψ, ρ, q = m/n, dq/dψ) + - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix` + - Asymptotic solution bases at the inner-layer boundaries ### Perturbed Equilibrium - `PerturbedEquilibriumControl` - User-facing TOML configuration parameters diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md index 46ca5f44d..633df587e 100644 --- a/docs/src/vacuum.md +++ b/docs/src/vacuum.md @@ -74,7 +74,7 @@ wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings( ) # Compute vacuum response matrix -wv, grri, xzpts = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings) +wv, _, _, _ = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings) ``` ### Vacuum Field Calculation at Observation Points diff --git a/docs/src/workflow.md b/docs/src/workflow.md index d872dcea6..0aa3c8433 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -59,8 +59,7 @@ The single `gpec.toml` file supplies user-selected options to every module. The **Outputs**: - `wv` — Vacuum response matrix (scaled by the singular factor (m - nq)(m' - nq), see Chance 1997) -- `grri` — Interior Green's function matrix (plasma boundary → plasma boundary) -- `grre` — Exterior Green's function matrix (plasma boundary → wall) +- `I_v` — Vacuum surface-current matrix Iᵛ when `compute_Iv=true` (otherwise zeros); PerturbedEquilibrium inverts this to surface inductance `L` **Key references**: [Chance et al. (1997)](citations.md#Vacuum-Module), [Chance et al. (2007)](citations.md#Vacuum-Module) @@ -182,5 +181,5 @@ All results are written to a single HDF5 file (default: `gpec.h5`). The file is | `locstab/` | Local stability: Mercier criterion, shear | | `integration/` | ODE integration results: energy matrices, eigenvalues | | `singular/` | Per-surface data: ψ_s, m/n, Δ', small solution coefficients | -| `vacuum/` | Vacuum response matrices: wv, grri, grre | +| `vacuum/` | Vacuum response matrices: wv (and I_v when computed for PE) | | `perturbed/` | Perturbed equilibrium: ξ, b in mode space, island diagnostics | diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index d7569a437..0d86696a5 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -11,8 +11,6 @@ A mutable struct holding data related to the singular surfaces in the equilibriu - `n::Vector{Int}` - Toroidal mode number(s) - `q::Float64` - Safety factor (= m/n) - `q1::Float64` - Derivative of safety factor with respect to ψ - - `grri::Array{ComplexF64,2}` - Interior Green's function at this surface [mthvac, mpert] - - `grre::Array{ComplexF64,2}` - Exterior Green's function at this surface [mthvac, mpert] - `delta_prime::Vector{ComplexF64}` - **STUB (not physically valid)**. Per-surface ca-based Δ' estimate retained for future work / debugging only. The physically valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`, computed via the STRIDE global BVP (Glasser 2018 PoP 25, 032501). Do not use this field for tearing-stability analysis; do not expect agreement with `delta_prime_matrix`. - `delta_prime_col::Matrix{ComplexF64}` - **STUB (not physically valid)**. Per-surface ca-based Δ' column retained for future work / debugging only. Shape (numpert_total × n_res_modes); `delta_prime_col[j, i] = (ca_r[j,ipert_res_i,2] - ca_l[j,ipert_res_i,2]) / (4π²·psio)`. The diagonal element matches the (also stubbed) `delta_prime[i]`. Only populated for the Riccati/parallel FM paths. The physically valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`; this field exists for future development on intra-surface coupling diagnostics, not for production use. """ @@ -23,8 +21,6 @@ A mutable struct holding data related to the singular surfaces in the equilibriu n::Vector{Int} = Int[] q::Float64 = 0.0 q1::Float64 = 0.0 - grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 0, 0) - grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 0, 0) delta_prime::Vector{ComplexF64} = ComplexF64[] delta_prime_col::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) ua_left::Array{ComplexF64,3} = Array{ComplexF64}(undef, 0, 0, 0) # asymptotic basis at left inner-layer boundary @@ -422,8 +418,6 @@ Populated in `Free.jl`. - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N): power-normalized and invariant to the working (Jacobian) coordinate; et = ep + ev per mode - `n_tor_idx::Vector{Int}` - 0-based toroidal mode number index of each sorted eigenvalue (numpert_total). Needed in `write_imas` - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero - - `grri::Array{ComplexF64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × numpert_total) - - `grre::Array{ComplexF64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × numpert_total) - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z) - `wall_pts::Array{Float64, 3}` - Cartesian coordinates of wall points, shape (mthvac * nzvac) × 3 for (x, y, z) """ @@ -441,8 +435,6 @@ Populated in `Free.jl`. et::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total) n_tor_idx::Vector{Int} = zeros(Int, numpert_total) vacuum_eigenvalue::Float64 = NaN - grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) - grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) plasma_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) wall_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) end diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 52acb9ba2..0ce2af361 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -175,7 +175,7 @@ q-window minimum. # Compute raw vacuum matrix at the actual scan psi (singfac NOT applied; free_compute_total applies it analytically) vac_inputs = Vacuum.VacuumInput(equil, psi_array[i], ctrl.mthvac, ctrl.nzvac, intr.mlow:intr.mhigh, intr.nlow:intr.nhigh) - wv, _, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) + wv, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) @views wv_array[i, :, :] .= wv end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index f6b1c7b41..dd69d6273 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -461,12 +461,10 @@ function main_from_inputs( @info "Computing free boundary energies ($wall_desc)" end vac_data = free_run!(odet, ctrl, equil, ffit, intr) - if real(vac_data.et[1]) < 0 - if ctrl.verbose + if ctrl.verbose + if real(vac_data.et[1]) < 0 @warn "Free-boundary mode unstable for n = $nstring" - end - else - if ctrl.verbose + else @info "All free-boundary modes stable for n = $nstring" end end @@ -554,6 +552,30 @@ function main_from_inputs( slayer=slayer_result) end + # No perturbed equilibrium calculations if vacuum data is not available + if !ctrl.vac_flag + @warn "Vacuum data not available. Skipping perturbed equilibrium calculations. Set vac_flag=true in [ForceFreeStates] section." + return + end + + # No perturbed equilibrium calculations if free-boundary mode is unstable + if real(vac_data.et[1]) < 0 + @warn "Since a free-boundary mode is unstable, perturbed equilibrium calculations will not run." + return + end + + # No perturbed equilibrium calculations if wall shape is not nowall + if intr.wall_settings.shape !== "nowall" + @warn "Perturbed equilibrium calculations are only supported for nowall cases - skipping perturbed equilibrium calculations." + return + end + + # Perturbed equilibrium calculations still need multi-n support added + if intr.npert > 1 + @warn "Perturbed equilibrium calculations are only supported for single-n cases - skipping perturbed equilibrium calculations." + return + end + # ---------------------------------------------------------------- # Perturbed Equilibrium # ---------------------------------------------------------------- @@ -672,17 +694,11 @@ function main_from_inputs( end slayer_result = _run_slayer_stage(pe_file) - # ---------------------------------------------------------------- - # Done - # ---------------------------------------------------------------- @info "\n$_BANNER\n GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER" - # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found - return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, vac_data=ctrl.vac_flag ? vac_data : nothing, slayer=slayer_result) - end """ @@ -693,10 +709,6 @@ This combines the functionality of several pieces of the Fortran code in `ode_ou primarily `ode_output_open` and the various `bin_euler` writes that occur throughout the integration. Some parameters are only dumped in their respective flags are true, e.g. vacuum data if `vac_flag` is true. - -### TODOs - -Combine spline unpacking if possible, too many extra lines """ function write_outputs_to_HDF5( ctrl::ForceFreeStatesControl, diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 5f1fe83b5..6fe7465a9 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -1,5 +1,7 @@ module PerturbedEquilibrium +const μ0 = 4π * 1e-7 + # Imports using HDF5 using Printf @@ -78,7 +80,7 @@ function compute_perturbed_equilibrium( state = PerturbedEquilibriumState() - # Step 0: Initialize mode arrays for convenient indexing + # Initialize mode arrays for convenient indexing initialize_mode_arrays!(intr, ffs_intr) # Load forcing data. On the gpec.h5 replay path the caller preloads @@ -113,27 +115,16 @@ function compute_perturbed_equilibrium( end end - # Step 2: Compute plasma response + # Compute plasma response if ctrl.compute_response - if vac_data === nothing - @warn "Vacuum data not available. Skipping plasma response calculation. Set vac_flag=true in [ForceFreeStates] section." - else - compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit) - end + compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit) end - # Step 3: Compute singular coupling metrics + # Compute singular coupling metrics if ctrl.compute_singular_coupling - if vac_data === nothing - @warn "Vacuum data not available. Skipping singular coupling calculation. Set vac_flag=true in [ForceFreeStates] section." - else - compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl) - end + compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl) end - # Step 4: Output eigenmode fields (integrated into HDF5 output) - # This is handled by write_outputs_to_HDF5 in main() - return state end diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index ba2827246..1a0954daa 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -25,40 +25,31 @@ function compute_plasma_response!( metric::MetricData, ffit::FourFitVars ) - if ctrl.verbose - @info "Computing plasma response (wt0-based inductance)" - end + ctrl.verbose && @info "Computing plasma response" # Build flux matrix from ForceFreeStates eigenmodes [mode × eigenmode] - flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, vac_data, ffs_intr) + flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, ffs_intr) - # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default) - plasma_inductance = calc_plasma_inductance(vac_data, ffs_intr, equil.psio) + # Compute plasma inductance + plasma_inductance = calc_plasma_inductance(ffs_intr, vac_data.wt0, equil.psio) - # Surface inductance L from Green's functions at psilim. - # Requires a 2D (nzvac=1) vacuum response so rows are theta points only, + # Surface inductance L from vacuum surface-current matrix at psilim nn = ffs_intr.nlow vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, vac_data.mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn]) - wall_nowall = Vacuum.WallShapeSettings(; shape="nowall") - _, grri_2d_raw, grre_2d_raw, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall) - grri_2d = Matrix{ComplexF64}(grri_2d_raw) - grre_2d = Matrix{ComplexF64}(grre_2d_raw) - ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν - surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac) - permeability = calc_permeability(plasma_inductance, surface_inductance) - - # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats). - # Λ (plasma inductance) is not Hermitian — its anti-Hermitian part is the dissipative/torque response — so the adjoint matters. - L_inv = inv(surface_inductance) - reluctance = L_inv * (plasma_inductance' - surface_inductance) * L_inv + _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall; compute_Iv=true) + surface_inductance = calc_surface_inductance(I_v) - # Store permeability in internal state for singular coupling / field reconstruction. - # These consumers operate on the physical control-surface flux Φ_x, so the internal - # copy stays in flux space; only the stored/output quantities are conformed to fields below. + # Compute permeability P = Λ·L⁻¹ and store in internal state for singular coupling / field reconstruction. + permeability = plasma_inductance / surface_inductance intr.plasma_response = permeability + # Compute reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ + # Λ is not Hermitian — its anti-Hermitian part is the dissipative/torque response — so the adjoint matters. + L_inv = inv(surface_inductance) + reluctance = L_inv * (plasma_inductance' - surface_inductance) * L_inv + # Conform the control-surface matrices to the coordinate-invariant root-area-weighted - # field (b̃) space for output (issue #233 / Pharr 2026). Store the b̃→b̄ operator S = Σ/√A + # field (b̃) space for output (Pharr 2026). Store the b̃→b̄ operator S = Σ/√A # and the scalar surface area A so users can recover the area-weighted field (b̄ = S·b̃) or # flux (Φ = A·b̄) — see Utils.jl output docs. rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs_intr) @@ -70,11 +61,9 @@ function compute_plasma_response!( state.rootarea_to_area_weight = rootarea_to_area_weight state.surface_area = surface_area - # Forcing and response on the control surface. Flux Φ appears only as a brief internal bridge: - # forcing arrives as Φ_x, the field reconstruction below consumes Φ_tot, and the b̃ spectra are - # formed via the conform operator R = S·A (Φ = R·b̃). + # Compute actual flux from external flux and permiability, Φ = P Φ^x forcing_flux = map_forcing_to_eigenmodes(intr.forcing_modes, ffs_intr) - response_flux = compute_plasma_response_vector(permeability, forcing_flux) + response_flux = permeability * forcing_flux # Output forcing/response in the three Pharr field representations (all tesla): # b̃ (root-area-weighted) = R⁻¹·Φ, b (bare) = Σ⁻¹·b̃, b̄ (area-weighted) = S·b̃. @@ -119,7 +108,5 @@ function compute_plasma_response!( state.b_n_modes = b_n_modes state.xi_n_modes = xi_n_modes - if ctrl.verbose - @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))" - end + ctrl.verbose && @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))" end diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index a4b798a6d..9faf6a483 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -11,283 +11,85 @@ Reference: [Park Phys. Plasmas 2009 056115] using ..Utilities.FourierTransforms """ - extract_boundary_displacements( - equil::Equilibrium.PlasmaEquilibrium, - ForceFreeStates_results::OdeState, - intr::ForceFreeStatesInternal - )::NamedTuple - -Extract eigenmode displacements and equilibrium quantities at the plasma boundary. - -This function extracts the data needed to compute the normal magnetic field at the -plasma surface from ForceFreeStates eigenmode solutions. - -## What's extracted: + build_flux_matrix(equil, ffs_results, intr) -> Matrix{ComplexF64} - 1. **Boundary displacement**: ξ_ψ from `u_store[:, :, 1, end]` +Calculate the vacuum poloidal flux matrix at the plasma boundary: +bwp_mn[i,j] = 1im * (dΨ/dρ) * (m - n * q) * ξ_ψ(i,j) - + This is the radial (normal) component of the eigenmode displacement - + At the last radial integration point (plasma edge) - + Dimensions: [numpert_total, numpert_total] - - 2. **Flux surface spacing**: dΨ/dρ at boundary - - + From equilibrium bicubic spline evaluation - + Needed to convert displacement to magnetic field - - 3. **Safety factor**: q at boundary - - + Used to compute singular factors (m - n*q) - + Identifies resonant surfaces - -## Arguments +Arguments: - - `equil`: Equilibrium solution containing flux surfaces and q-profile - - `ForceFreeStates_results`: ODE integration results containing u_store with eigenmodes - - `intr`: ForceFreeStates internal state with boundary location (psilim) - -## Returns + - `equil`: Plasma equilibrium + - `ffs_results`: ForceFreeStates eigenmode ODE results + - `intr`: ForceFreeStates internal state with mode information -Named tuple with: +Returns: - - `ξ_psi_boundary`: Boundary displacement [numpert_total, numpert_total] - - `dPsi_drho`: Flux surface spacing at boundary (scalar) - - `q_boundary`: Safety factor at boundary (scalar) - - `psi_boundary`: Normalized flux at boundary (scalar) + - `bwp_mn`: Complex flux matrix [mode_i, eigenmode_j] """ -function extract_boundary_displacements( - equil::Equilibrium.PlasmaEquilibrium, - ForceFreeStates_results::OdeState, - intr::ForceFreeStatesInternal -) - # Extract boundary displacement (normal component) - # u_store dimensions: [numpert_total, numpert_total, 2, numsteps] - # Index 1 in 3rd dimension is ξ_ψ (radial displacement) - # Last index in 4th dimension is the boundary - ξ_psi_boundary = ForceFreeStates_results.u_store[:, :, 1, ForceFreeStates_results.step] - - # Get boundary location in normalized flux coordinates - psi_boundary = ForceFreeStates_results.psi_store[ForceFreeStates_results.step] - - # Evaluate equilibrium quantities at boundary - # Safety factor at boundary - q_boundary = ForceFreeStates_results.q_store[ForceFreeStates_results.step] - - # FFS ODE integrates in ψ (normalized flux), so bwp_mn = chi1·singfac·2πi·ξ_ψ - # where chi1 = 2π·psio (Fortran idcon.f: chi1 = twopi*psio) - # Combined flux factor = chi1·2π = (2π)²·psio (gpeq.f: bwp_mn = chi1·singfac·twopi·ifac·xsp) - dPsi_drho = (2π)^2 * equil.psio +function build_flux_matrix(equil::Equilibrium.PlasmaEquilibrium, ffs_results::OdeState, intr::ForceFreeStatesInternal)::Matrix{ComplexF64} - return ( - ξ_psi_boundary=ξ_psi_boundary, - dPsi_drho=dPsi_drho, - q_boundary=q_boundary, - psi_boundary=psi_boundary - ) -end - -""" - compute_normal_magnetic_field( - boundary_data::NamedTuple, - intr::ForceFreeStatesInternal - )::Matrix{ComplexF64} - -Compute normal magnetic field at plasma boundary from eigenmode displacements. - -This is the key step that converts eigenmode displacements to magnetic field perturbations -at the plasma surface. From the ideal MHD constraint [Park Phys. Plasmas 2009 056115 eq. 4]: - - B_n = i * (dΨ/dρ) * (m - n*q) * ξ_ψ - -where ξ_ψ is the radial displacement eigenfunction. - -## Physical Interpretation [Park Phys. Plasmas 2007 052110 Section II]: - - - ξ_ψ[i,j]: Displacement of mode i due to eigenmode j - - singfac[i] = m[i] - n*q: Singular factor measuring distance from rational surface - - dΨ/dρ: Converts displacement to flux perturbation (poloidal flux gradient) - - Factor of i: Phase relationship for oscillating fields in complex representation - -## Arguments - - - `boundary_data`: Output from extract_boundary_displacements() - + ξ_psi_boundary: Boundary displacement [numpert_total, numpert_total] - + dPsi_drho: Flux surface spacing at boundary (scalar) - + q_boundary: Safety factor at boundary (scalar) - + psi_boundary: Normalized flux at boundary (scalar) - - `intr`: ForceFreeStates internal state with mode arrays (mlow, mhigh, nlow, etc.) - -## Returns + # Extract boundary displacements and equilibrium quantities + ξ_psi_boundary = ffs_results.u_store[:, :, 1, ffs_results.step] + q_boundary = ffs_results.q_store[ffs_results.step] + dPsi_drho = (2π)^2 * equil.psio - - `bwp_mn[numpert_total, numpert_total]`: Normal magnetic field matrix where bwp_mn[i,j] - is the normal field of Fourier mode i in response to eigenmode j -""" -function compute_normal_magnetic_field( - boundary_data::NamedTuple, - intr::ForceFreeStatesInternal -)::Matrix{ComplexF64} - - numpert_total = intr.numpert_total - bwp_mn = zeros(ComplexF64, numpert_total, numpert_total) - - # Extract boundary data - ξ_psi = boundary_data.ξ_psi_boundary - dPsi_drho = boundary_data.dPsi_drho - q_boundary = boundary_data.q_boundary - - # Compute singular factor for each Fourier mode [Park Phys. Plasmas 2009 056115 eq. 4] - # singfac[i] = m[i] - n*q_boundary measures distance from rational surface - # Mode indexing: modes are ordered as (m, n) pairs - # Linear index i corresponds to: m = (i-1) % mpert + mlow, n = (i-1) ÷ mpert + nlow - singfac = zeros(Float64, numpert_total) - for i in 1:numpert_total - m_mode = (i - 1) % intr.mpert + intr.mlow - n_mode = (i - 1) ÷ intr.mpert + intr.nlow - singfac[i] = m_mode - n_mode * q_boundary - end + # Compute singular factor for each Fourier mode + singfac = vec((intr.mlow:intr.mhigh) .- q_boundary .* (intr.nlow:intr.nhigh)') - # Compute normal magnetic field [Park Phys. Plasmas 2009 056115 eq. 4] - # bwp_mn[i,j] = i * (dΨ/dρ) * singfac[i] * ξ_ψ[i,j] - for i in 1:numpert_total - for j in 1:numpert_total - bwp_mn[i, j] = 1im * dPsi_drho * singfac[i] * ξ_psi[i, j] - end - end + # Compute normal magnetic field at plasma boundary [Park Phys. Plasmas 2009 056115 eq. 4] + bwp_mn = 1im * dPsi_drho .* singfac .* ξ_psi_boundary return bwp_mn end """ - build_flux_matrix( - equil::Equilibrium.PlasmaEquilibrium, - ForceFreeStates_results::OdeState, - vac_data::VacuumData, - intr::ForceFreeStatesInternal - )::Matrix{ComplexF64} - -Build vacuum poloidal flux matrix from ForceFreeStates eigenmode solutions. - -This extracts the vacuum flux response for each eigenmode at the plasma boundary. -In GPEC, this comes from `bwp_mn` (boundary normal field) computed from eigenmode -displacements. + calc_plasma_inductance(ffs_intr, wt0, psio) -> Matrix{ComplexF64} -The flux matrix relates eigenmode displacements to vacuum poloidal flux: +Compute the plasma inductance Λ from the displacement-space energy matrix ``W_0 = W_p + W_v`` +(Fortran `gpresp.f`, `resp_induct_flag=TRUE`): - 1. Extract eigenmode displacement at plasma boundary from u_store - 2. Compute normal magnetic field: B_ψ = i×(dΨ/dρ)×(m - n×q)×ξ_ψ - 3. Result is flux[mode_i, eigenmode_j] = bwp_mn[i,j] + 1. Normalize `wt0` to SI units (`ψ₀²/(μ₀·2)`, as in `idcon.f`). + 2. Convert ξ → Φ via ``Λ⁻¹ = 2·T₁·W₀·T₂`` with ``T_{ii} ∝ 1/(χ₁ s_i 2π)`` and + ``s_i = m_i - n·q_lim``; the ``s_i s_j`` factor undoes the `wv` scaling in `free_run!`. + 3. Invert to obtain Λ ## Arguments - - `equil`: Equilibrium solution containing flux surfaces and q-profile - - `ForceFreeStates_results`: ForceFreeStates ODE integration results containing eigenmodes - - `vac_data`: Vacuum response data from free boundary calculation - - `intr`: ForceFreeStates internal state with mode information + - `ffs_intr`: ForceFreeStates internal state (`mlow`, `mhigh`, `nlow`, `nhigh`, `qlim`) + - `wt0`: Total energy matrix `wp + wv` [numpert_total × numpert_total] + - `psio`: Toroidal flux [Wb/rad] ## Returns - - `flxmats[numpert_total, numpert_total]`: Complex flux matrix where flxmats[i,j] is the - vacuum flux of mode i in response to eigenmode j -""" -function build_flux_matrix( - equil::Equilibrium.PlasmaEquilibrium, - ForceFreeStates_results::OdeState, - vac_data::VacuumData, - intr::ForceFreeStatesInternal -)::Matrix{ComplexF64} - - # Step 1: Extract boundary displacements and equilibrium quantities - boundary_data = extract_boundary_displacements(equil, ForceFreeStates_results, intr) - - # Step 2: Compute normal magnetic field at plasma boundary - # This is the actual implementation of GPEC's bwp_mn calculation - # bwp_mn[i,j] = i * (dΨ/dρ) * (m[i] - n*q_boundary) * ξ_ψ[i,j] - flxmats = compute_normal_magnetic_field(boundary_data, intr) - - return flxmats -end - + - Plasma inductance Λ matrix [numpert_total × numpert_total] """ - calc_plasma_inductance( - vac_data::VacuumData, - ffs_intr::ForceFreeStatesInternal, - psio::Float64 - )::Matrix{ComplexF64} - -Calculate plasma inductance matrix Λ using the wt0-based energy formula -(matches Fortran `gpresp_induct` with `resp_induct_flag=TRUE`). - - Λ = inv(2·t₁·wt0·t₂) - -where t₁ = im/(χ₁·s_i·2π), t₂ = -im/(χ₁·s_j·2π), s_i = m_i - n·q_lim. +function calc_plasma_inductance(ffs_intr::ForceFreeStatesInternal, wt0::Matrix{ComplexF64}, psio::Float64)::Matrix{ComplexF64} -Note: `vac_data.wt0` already contains singfac² factors (s_i·s_j) baked into the -vacuum term via the scaling in `free_run!`. The t₁/t₂ factors divide by s_i·s_j, -correctly recovering the properly-normalized inductance. - -## Arguments - - - `vac_data`: Vacuum data containing wt0 (total energy matrix before eigenvector sorting) - - `ffs_intr`: ForceFreeStates internal state with mode info (mlow, mpert, nlow, qlim) - - `psio`: Total toroidal flux [Wb/rad] from equilibrium (equil.psio) - -## Returns - - - Plasma inductance matrix Lambda [numpert_total × numpert_total] -""" -function calc_plasma_inductance( - vac_data::VacuumData, - ffs_intr::ForceFreeStatesInternal, - psio::Float64 -)::Matrix{ComplexF64} + # Singular factors s_i = m_i - n*qlim (same as Fortran: mfac(i) - nn*qlim) + singfac = vec((ffs_intr.mlow:ffs_intr.mhigh) .- ffs_intr.qlim .* (ffs_intr.nlow:ffs_intr.nhigh)') - mpert = ffs_intr.numpert_total - chi1 = 2π * psio # = Fortran's chi1 = twopi*psio - n = ffs_intr.nlow - qlim = ffs_intr.qlim # q at psilim + # Convert to metric units + wt0_norm = wt0 .* (psio^2 / (μ0 * 2)) - # Singular factors s_i = m_i - n*qlim (same as Fortran: mfac(i) - nn*qlim) - s = [((i-1) % ffs_intr.mpert + ffs_intr.mlow) - n * qlim for i in 1:mpert] - - # Fortran idcon_norm: wt0 = wt0/(mu0*2)*psio^2 - # Julia's vac_data.wt0 is raw wp+wv; Fortran additionally scales by psio^2/(mu0*2) - mu0 = 4π * 1e-7 - wt0_norm = vac_data.wt0 .* (psio^2 / (mu0 * 2)) - - # Build temp2[i,j] = 2·t1_i·wt0[i,j]·t2_j (matches Fortran gpresp_induct) - temp2 = Matrix{ComplexF64}(undef, mpert, mpert) - for i in 1:mpert, j in 1:mpert - t1 = im / (chi1 * s[i] * 2π) - t2 = -im / (chi1 * s[j] * 2π) - temp2[i, j] = 2 * t1 * wt0_norm[i, j] * t2 - end + # Convert from displacement to flux space using 1 / (singfac * chi1 * 2π) factor + # The factor of 2 comes from the 1/2 + inverse in eq. 31 of Park Phys. Plasmas 2007 + chi1 = 2π * psio + wt0_norm .*= 2 ./ (2π * chi1)^2 ./ (singfac' .* singfac) - return inv(temp2) + return inv(wt0_norm) end """ - calc_permeability( - plasma_inductance::Matrix{ComplexF64}, - surface_inductance::Matrix{ComplexF64} - )::Matrix{ComplexF64} - -Calculate permeability matrix P = Λ·L⁻¹ (matches Fortran `gpresp_permeab`). + calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} -## Arguments - - - `plasma_inductance`: Plasma inductance matrix Lambda - - `surface_inductance`: Surface inductance matrix L - -## Returns - - - Permeability matrix P = Lambda * L^{-1} [mpert, mpert] +Surface inductance from the vacuum surface-current matrix, Park 2007 eq. 7: `Φ_x = L·I_v`. +Columns of `I_v` are driven by unit flux harmonics which are consistent since the eq. 3 +weight `1/(J|∇ψ|)` cancels the Jacobian in the vacuum solver's source density — so `Φ_x = 𝕀` +and `L = I_v⁻¹`, carrying the `μ₀(2π)²` normalization. """ -function calc_permeability( - plasma_inductance::Matrix{ComplexF64}, - surface_inductance::Matrix{ComplexF64} -)::Matrix{ComplexF64} - # P = Lambda * L^{-1} (right-division solves for P s.t. P*L = Lambda) - return plasma_inductance / surface_inductance +function calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} + μ₀ = 4π * 1e-7 + return inv(I_v) .* (μ₀ * (2π)^2) end """ @@ -404,50 +206,13 @@ function map_forcing_to_eigenmodes( numpert_total = intr.mpert * intr.npert forcing_vector = zeros(ComplexF64, numpert_total) - # Create mode index map: (m,n) -> linear index + # Combine all forcing modes into a single vector for forcing_mode in forcing_modes - # Find matching mode in eigenmode basis - for i in 1:numpert_total - # Calculate m and n for this index - # Using 0-based indexing converted to 1-based: - # m = (i-1) % mpert + mlow - # n = (i-1) ÷ mpert + nlow - m_mode = (i - 1) % intr.mpert + intr.mlow - n_mode = (i - 1) ÷ intr.mpert + intr.nlow - - if m_mode == forcing_mode.m && n_mode == forcing_mode.n - forcing_vector[i] = forcing_mode.amplitude - break - end + i = (forcing_mode.m - intr.mlow) + (forcing_mode.n - intr.nlow) * intr.mpert + 1 + if 1 <= i <= numpert_total + forcing_vector[i] = forcing_mode.amplitude end end return forcing_vector end - -""" - compute_plasma_response_vector( - permeability::Matrix{ComplexF64}, - forcing_vector::Vector{ComplexF64} - )::Vector{ComplexF64} - -Compute plasma response to external forcing. - -Response = Permeability * Forcing - -## Arguments - - - `permeability`: Permeability matrix - - `forcing_vector`: External forcing in eigenmode basis - -## Returns - - - Plasma response vector in eigenmode basis -""" -function compute_plasma_response_vector( - permeability::Matrix{ComplexF64}, - forcing_vector::Vector{ComplexF64} -)::Vector{ComplexF64} - - return permeability * forcing_vector -end diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index d9cd1ebc9..e2f64b3b1 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -40,9 +40,6 @@ function _hermite_cubic_val(u_a, u_b, du_a, du_b, psi_a, psi_b, psi) return @. h00 * u_a + h * h10 * du_a + h01 * u_b + h * h11 * du_b end -# Reflect a periodic theta-space vector θ → -θ (the theta reversal in gpvacuum_flxsurf). -_reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) - """ compute_singular_coupling_metrics!( state::PerturbedEquilibriumState, @@ -99,7 +96,6 @@ function compute_singular_coupling_metrics!( end chi1 = 2π * equil.psio - twopi = 2π mtheta = vac_data.mthvac wall_settings = Vacuum.WallShapeSettings(; shape="nowall") @@ -173,26 +169,18 @@ function compute_singular_coupling_metrics!( continue end - # Compute Green's functions at this surface for this n (once per pair) + # Surface-current matrix at this surface for this n (once per pair) vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn]) - _, grri_raw, grre_raw, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings) - grri = Matrix{ComplexF64}(grri_raw) - grre = Matrix{ComplexF64}(grre_raw) - - # Get ν on the vacuum theta grid (same ν used in the vacuum Fourier basis computation) - ν_vac = Vacuum.PlasmaGeometry(vac_input).ν + _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_Iv=true) + L_surf = calc_surface_inductance(I_v) - # Precompute L_surf; only the (m_res, m_res) diagonal element is needed for singflx - L_surf = compute_surface_inductance_from_greens(grri, grre, ffs_intr, nn, ν_vac) + # Only the (m_res, m_res) diagonal element is needed for singflx m_idx = m_res - mlow + 1 L_mm = L_surf[m_idx, m_idx] j_c = compute_current_density(equil, sing_surf.psifac) area = compute_surface_area(equil, sing_surf.psifac) - # Matches Fortran gpout_resp: shear = m*dq/dψ / q² = n*dq/dψ / q (since m=n*q). - # Uses abs(nn) because island_half_width = sqrt(abs(island_width_sq)), so the sign - # of shear only affects the sign of C_island_width_sq, not the physical island width. - shear = abs(nn) * sing_surf.q1 / sing_surf.q + shear = nn * sing_surf.q1 / sing_surf.q # m*dq/dψ / q² = n*dq/dψ / q (since m=n*q) # Evaluate bwp1_mn = ∂b^ψ/∂ψ at lpsi and rpsi using permeability-weighted eigenstates. # Matches Fortran gpout_resp: evaluate bwp1_mn at lpsi/rpsi via gpeq_sol @@ -265,14 +253,14 @@ function compute_singular_coupling_metrics!( # - resonant (shielding) current: j_c already integrates jac·|∇ψ| over the surface, so the # Jacobian weighting is carried inside j_c — no separate area factor needed. # - resonant flux → field: Φ^r/A^r [T], invariant [Park 2008; Pharr 2026]. - state.C_delta_prime[row, :] = jump_vec ./ (twopi * chi1) - state.C_resonant_current[row, :] = jump_vec .* (-j_c / (twopi * m_res)) + state.C_delta_prime[row, :] = jump_vec ./ (2π * chi1) + state.C_resonant_current[row, :] = jump_vec .* (-j_c / (2π * m_res)) # Matches Fortran gpout_resp: singflx = L·fkaxmn, resonant area-weighted field = singflx/area, # islandhwids = 4·singflx/(2π·shear·q·chi1) - singflx_pre = (L_mm / (twopi * nn)) .* state.C_resonant_current[row, :] + singflx_pre = (L_mm / (2π * nn)) .* state.C_resonant_current[row, :] state.C_resonant_area_weighted_field[row, :] = singflx_pre ./ area if abs(shear) > 1e-10 - state.C_island_width_sq[row, :] = (4.0 / (twopi * shear * sing_surf.q * chi1)) .* singflx_pre + state.C_island_width_sq[row, :] = abs.(4.0 / (2π * shear * sing_surf.q * chi1)) .* singflx_pre end state.rational_psi[row] = sing_surf.psifac @@ -343,24 +331,6 @@ j_c = χ₁² * q / (μ₀ * integral) where the integral is computed via flux surface integration: integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ -## GPEC Formula - -```fortran -DO itheta=0,mthsurf - CALL bicube_eval(rzphi,respsi,theta(itheta),1) - rfac=SQRT(rzphi%f(1)) - jac=rzphi%f(4) - w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac - w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) - delpsi(itheta)=SQRT(w(1,1)**2+w(1,2)**2) - sqreqb(itheta)=(sq%f(1)**2+chi1**2*delpsi(itheta)**2)/(twopi*r(itheta))**2 - jcfun(itheta)=sqreqb(itheta)/(delpsi(itheta)**3) - j_c(ising)=j_c(ising)+jac*delpsi(itheta)*jcfun(itheta)/mthsurf -ENDDO -j_c(ising)=j_c(ising)-jac*delpsi(mthsurf)*jcfun(mthsurf)/mthsurf ! trapezoidal rule -j_c(ising)=1.0/j_c(ising)*chi1**2*sq%f(4)/mu0 -``` - ## Implementation Uses trapezoidal rule integration around the flux surface with metric quantities @@ -375,16 +345,12 @@ function compute_current_density( psi::Float64 )::Float64 # Physical constants - μ₀ = 4π * 1e-7 chi1 = 2π * equil.psio - twopi = 2π # Get equilibrium quantities at this surface F_tor = equil.profiles.F_spline(psi) # Toroidal field function (2π·R·B_tor in GPEC convention) q = equil.profiles.q_spline(psi) # Safety factor - ro = equil.ro - # Number of theta points for integration # Match GPEC's mthsurf (typically 101 points from theta=0 to theta=1) mthsurf = length(equil.rzphi_ys) - 1 @@ -407,7 +373,7 @@ function compute_current_density( delpsi = m.delpsi # flux gradient magnitude |∇ψ| # sqreqb = (F² + χ₁²|∇ψ|²) / (2πR)² where F = R·B_tor (Fortran sq%f(1)) - sqreqb = (F_tor^2 + chi1^2 * delpsi^2) / (twopi * m.r)^2 + sqreqb = (F_tor^2 + chi1^2 * delpsi^2) / (2π * m.r)^2 # Integrand function jcfun = sqreqb / (delpsi^3) @@ -426,105 +392,12 @@ function compute_current_density( # Trapezoidal rule end correction (subtract half of last point contribution) integral -= last_jac * last_delpsi * last_jcfun / mthsurf - # Final normalization: j_c = (1/integral) * χ₁² * q / μ₀ - j_c = (1.0 / integral) * chi1^2 * q / μ₀ + # Final normalization: j_c = (1/integral) * χ₁² * q / μ0 + j_c = (1.0 / integral) * chi1^2 * q / μ0 return j_c end -""" - compute_surface_inductance_from_greens( - grri::Matrix{ComplexF64}, - grre::Matrix{ComplexF64}, - ffs_intr::ForceFreeStatesInternal, - nn::Int, - ν::Vector{Float64} - )::Matrix{ComplexF64} - -Compute surface inductance matrix from Green's functions at flux surface. - -Implements the GPEC `gpvacuum_flxsurf` algorithm. - -The Julia vacuum code uses SFL Fourier basis `cos(m*θ - n*ν)` in the column transform, -so the row DFT must apply the matching toroidal phase correction `exp(-i*n*ν)` before -the DFT (matching Fortran `gpvacuum_flxsurf`'s `EXP(-ifac*nn*dphi)` phase correction). - -## Arguments - - - `grri`: Interior Green's function [mtheta, mpert] - - `grre`: Exterior Green's function [mtheta, mpert] - - `ffs_intr`: ForceFreeStates internal state - - `nn`: Toroidal mode number - - `ν`: Toroidal angle offset on the vacuum theta grid [mtheta] - -## Returns - -Surface inductance matrix [mpert × mpert] -""" -@with_pool pool function compute_surface_inductance_from_greens( - grri::Matrix{ComplexF64}, - grre::Matrix{ComplexF64}, - ffs_intr::ForceFreeStatesInternal, - nn::Int, - ν::Vector{Float64} -)::Matrix{ComplexF64} - mpert = ffs_intr.mpert - mtheta = length(ν) - μ₀ = 4π * 1e-7 - - ft = FourierTransforms.FourierTransform(mtheta, mpert, ffs_intr.mlow) - - flux_matrix = zeros!(pool, ComplexF64, mpert, mpert) - current_matrix = zeros!(pool, ComplexF64, mpert, mpert) - - kax = zeros!(pool, ComplexF64, mtheta) - grri_surf = @view grri[1:mtheta, :] - grre_surf = @view grre[1:mtheta, :] - - # Toroidal phase correction: exp(-i*n*ν) - phase = cis.(-nn .* ν) - - for i in 1:mpert - flux_matrix[i, i] = 1.0 - - # Complex grri/e stores exp(i(mθ-nν)) projection, need conjugate for exp(-i(mθ-nν)) - kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ₀ * (2π)^2) - - # Port of Fortran gpvacuum_flxsurf: apply toroidal phase, reverse theta, forward-DFT. - g_phased = kax .* phase - current_matrix[:, i] = ft(_reverse_theta(g_phased)) - end - - # Compute surface inductance: L_surf = flux * inv(current) = inv(current) - L_surf = zeros(ComplexF64, mpert, mpert) - - current_mag = maximum(abs.(current_matrix)) - - if current_mag < 1e-15 - @warn "Current matrix is all zeros! Cannot compute surface inductance." maxlog=1 - for i in 1:mpert - L_surf[i, i] = μ₀ * 1e-6 - end - else - try - regularization = 1e-12 * current_mag - current_reg = current_matrix + regularization * I - - L_surf = flux_matrix * inv(current_reg) - - # Hermitianize (matches Fortran: temp1 = 0.5*(temp1 + CONJG(TRANSPOSE(temp1)))) - L_surf = 0.5 * (L_surf + L_surf') - catch e - @warn "Surface inductance inversion failed: $e" maxlog=1 - for i in 1:mpert - L_surf[i, i] = μ₀ * 1e-6 - end - end - end - - return L_surf -end - """ compute_surface_area( equil::Equilibrium.PlasmaEquilibrium, diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index 06ee4a6a0..fc8bd0ed0 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -73,7 +73,8 @@ function VacuumInput( r, z, ν = extract_plasma_surface_at_psi(equil, ψ) # Remove the last point to go from the [0, 2π] grid to VACUUM's [0, 2π) grid - # and reverse the arrays for VACUUM's CW θ direction + # and reverse the arrays for VACUUM's CW θ direction (θ_VAC = -θ_GPEC). This handedness is why + # operators returned to GPEC (e.g. the surface-inductance current matrix) are conjugated. return VacuumInput(; x=reverse(r)[1:(end-1)], z=reverse(z)[1:(end-1)], diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index b060f59fc..f0cbd3ad6 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -24,50 +24,45 @@ export extract_plasma_surface_at_psi export PlasmaGeometry # Relative anti-Hermitian residual above which we warn that the vacuum grid should be refined. -const _WV_HERMITICITY_WARN_TOL = 1e-4 +const _HERMITICITY_WARN_TOL = 1e-4 """ - _symmetrize_vacuum_energy!(wv) + _warn_and_symmetrize!(mat, name) -Enforce Hermiticity of the vacuum energy matrix Wᵛ in place. - -Wᵛ is the generator of the vacuum magnetic energy, δW_v = ξ† Wᵛ ξ. Because that energy is a real -number for every perturbation ξ, the exact operator is Hermitian. The finite-resolution -boundary-integral quadrature (finite `mtheta`/`nzeta`) breaks exact Hermiticity, leaving a small -anti-Hermitian residual that is a pure discretization artifact and vanishes as the grid is refined. -We replace Wᵛ by its Hermitian part to restore this physical property, warning when the residual -is large enough that the vacuum grid should be refined. +Replace `mat` by its Hermitian part in place, warning first if the anti-Hermitian residual exceeds +`_HERMITICITY_WARN_TOL`. """ -function _symmetrize_vacuum_energy!(wv::AbstractMatrix) - - herm_norm = norm(wv + wv') +function _warn_and_symmetrize!(mat::AbstractMatrix, name::String) + herm_norm = norm(mat + mat') if herm_norm > 0 - # Relative anti-Hermitian residual ‖½(W−W†)‖/‖½(W+W†)‖ - rel_residual = norm(wv - wv') / herm_norm - if rel_residual > _WV_HERMITICITY_WARN_TOL - @warn "Vacuum energy matrix Wᵛ is non-Hermitian above tolerance $(rel_residual) > $(_WV_HERMITICITY_WARN_TOL) before " * + # Relative anti-Hermitian residual ‖½(M−M†)‖/‖½(M+M†)‖ + rel_residual = norm(mat - mat') / herm_norm + if rel_residual > _HERMITICITY_WARN_TOL + @warn "$name is non-Hermitian above tolerance $(rel_residual) > $(_HERMITICITY_WARN_TOL) before " * "symmetrization. Increase vacuum grid resolution to reduce it." end end - hermitianpart!(wv) + hermitianpart!(mat) end """ - _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) 2D (axisymmetric) vacuum response calculation. Each toroidal mode `n` decouples in 2D geometry, so the routine loops over `inputs.n_modes`, -building the double-/single-layer operators, solving the exterior and interior systems, and -filling the corresponding diagonal block of the response matrix and the matching column block -of the Green's functions. +building the double-/single-layer operators, solving the exterior system for `wv`, and +optionally the interior system to build `I_v` when `compute_Iv=true`. +Green's functions are internal scratch only. """ -@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) mpert = length(inputs.m_modes) + mlow = inputs.m_modes[1] num_points_surf = inputs.mtheta vac_data.wv .= 0 + compute_Iv && (vac_data.I_v .= 0) # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) @@ -75,13 +70,11 @@ of the Green's functions. # Loop over all decoupled toroidal modes for (idx_n, n) in enumerate(inputs.n_modes) - ft = FourierTransform(inputs.mtheta, mpert, inputs.m_modes[1]; n=n, ν=plasma_surf.ν) + ft = FourierTransform(inputs.mtheta, mpert, mlow; n=n, ν=plasma_surf.ν) - # Diagonal block of wv and matching column block of the Green's functions + # Diagonal block of wv (and I_v when requested) block_idx = ((idx_n-1)*mpert+1):(idx_n*mpert) wv_block = @view vac_data.wv[block_idx, block_idx] - grri_block = @view vac_data.grri[:, block_idx] - grre_block = @view vac_data.grre[:, block_idx] # Active rows for computation (plasma only if no wall, plasma+wall if wall present) num_points_total = wall.nowall ? num_points_surf : 2 * num_points_surf @@ -89,10 +82,7 @@ of the Green's functions. # Local work matrices grad_green = zeros!(pool, num_points_total, num_points_total) green_temp = zeros!(pool, num_points_surf, num_points_surf) - - # Views into output Green's function matrices for the active rows/columns - grre = @view grre_block[1:num_points_total, :] - grri = @view grri_block[1:num_points_total, :] + grre = zeros!(pool, ComplexF64, num_points_total, mpert) # Plasma–Plasma block compute_2D_kernel_matrices!(grad_green, green_temp, plasma_surf, plasma_surf, n) @@ -111,29 +101,45 @@ of the Green's functions. mul!(view(grre, (num_points_surf+1):num_points_total, :), green_temp, ft.basis') end - # Compute both Green's functions: exterior (kernelsign=+1) then interior (kernelsign=-1) - grri .= grre # start from same as exterior - grad_green_interior = similar!(pool, grad_green) - grad_green_interior .= grad_green - - # Solve exterior first, overwriting grad_green to save memory since we already have the interior kernel - ldiv!(lu!(grad_green), grre) - - # Interior flips the sign of the normal, but not the diagonal terms, so we multiply by -1 and add 2I to the diagonal - grad_green_interior .*= -1 - for i in 1:num_points_total - grad_green_interior[i, i] += 2.0 + if compute_Iv + # Copy RHS before exterior solve overwrites grre; keep a kernel copy for interior + grri = similar!(pool, grre) + grri .= grre + grad_green_interior = similar!(pool, grad_green) + grad_green_interior .= grad_green + + # Exterior operator D_ext = 2I + 𝒦 (Chance 1997 eq. 89); solve for the physical + # vacuum-outside potential grre = χ^(vo). Overwrites grad_green to save memory. + ldiv!(lu!(grad_green), grre) + + # Interior operator D_int = D_ext - 2I: the double-layer jump between the two one-sided + # boundary limits is 2I here, giving the vacuum-inside potential grri = χ^(vi). + for i in 1:num_points_total + grad_green_interior[i, i] -= 2.0 + end + ldiv!(lu!(grad_green_interior), grri) + + # Surface-current matrix, Park 2007 eq. 21b: μ₀I^v = χ^(vi) - χ^(vo) = grri - grre + # They are flipped because VACUUM builds the operators in its CW-θ frame while GPEC + # uses CCW-θ, flipping the outward-normal sign. + I_v_block = @view vac_data.I_v[block_idx, block_idx] + @views g_sum = grre[1:num_points_surf, :] .- grri[1:num_points_surf, :] + mul!(I_v_block, ft.basis, g_sum) + conj!(I_v_block) # Flip θ_VAC → -θ_VAC to get I^v in GPEC's CCW-θ frame. + I_v_block ./= num_points_surf + else + # Only need exterior system for wv + ldiv!(lu!(grad_green), grre) end - ldiv!(lu!(grad_green_interior), grri) # Project exterior kernel onto observer basis exp(-i*(mθ - nν)) and scale to get the response matrix mul!(wv_block, ft.basis, @view(grre[1:num_points_surf, :])) wv_block .*= 4π^2 / num_points_surf end - - # δW_v = ξ† Wᵛ ξ is real, so Wᵛ must be Hermitian; remove any residual from discretization - _symmetrize_vacuum_energy!(vac_data.wv) + # Remove any non-Hermitian residual from Hermitian matrices due to discretization + _warn_and_symmetrize!(vac_data.wv, "Wᵛ") + compute_Iv && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays @views begin @@ -147,7 +153,7 @@ of the Green's functions. end """ - _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) 3D (`inputs.nzeta > 1`) vacuum response via block-circulant field-period reduction. For `nfp == 1` the block-circulant assembly and residue-class loop are skipped in favour of a more efficient @@ -165,14 +171,18 @@ class needs one solve `wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \\ Ŝₖ)|_p block-row of the operators is built, so the kernel cost drops by `nfp` and the dense `O(N³)` factorization is replaced by per-class `O(M³)` solves (`M = N/nfp`). -Only `wv` is produced currently; `grri`/`grre` are returned zeroed and are not yet supported in the 3D path. +Only `wv` is produced currently; `I_v` is left zeroed when `compute_Iv=true` +(surface-current / inductance not yet supported in 3D). Extension point: per residue class, apply the per-period basis to `D̂ₖ⁻¹Ŝₖ` for the exterior columns and the interior variant `-D + 2I` for the interior columns, then scatter back into the `[2N × 2·num_modes]` arrays. """ -@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs fill!(vac_data.wv, 0) + compute_Iv && fill!(vac_data.I_v, 0) + + compute_Iv && @warn "compute_Iv=true is not supported for 3D vacuum response; I_v left as zeros" maxlog=1 # Full-torus geometry for source surface; observers are restricted to one field period full = expand_field_periods(inputs) @@ -245,11 +255,9 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end end - _symmetrize_vacuum_energy!(vac_data.wv) - - # Zero out the Green's function matrices (not tested in 3D yet) - fill!(vac_data.grri, 0) - fill!(vac_data.grre, 0) + # Remove any non-Hermitian residual from Hermitian matrices due to discretization + _warn_and_symmetrize!(vac_data.wv, "Wᵛ") + compute_Iv && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays vac_data.plasma_pts .= plasma_surf.r @@ -257,10 +265,10 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) + compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) -Allocate and return the vacuum response matrix and Green's functions for the given vacuum -inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the +Allocate and return the vacuum response matrix and optional surface-current matrix for the given +vacuum inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the output arrays for the full torus and forwards to the same 2D/3D workers. For performance-critical paths that already own preallocated storage (e.g. `ForceFreeStates.VacuumData`), prefer the in-place method to avoid extra heap allocations. @@ -268,28 +276,27 @@ in-place method to avoid extra heap allocations. # Returns - `wv`: complex vacuum response matrix (`num_modes × num_modes`). - - `grri`, `grre`: interior/exterior Green's functions (zeroed on the 3D nowall path). + - `I_v`: Vacuum surface-current matrix (`num_modes × num_modes`); zeros unless `compute_Iv=true`. - `plasma_pts`, `wall_pts`: surface coordinate arrays. """ -function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) num_points = inputs.mtheta * inputs.nzeta * inputs.nfp # mtheta for 2D num_modes = length(inputs.m_modes) * length(inputs.n_modes) vac = ( wv=zeros(ComplexF64, num_modes, num_modes), - grri=zeros(ComplexF64, 2 * num_points, num_modes), - grre=zeros(ComplexF64, 2 * num_points, num_modes), + I_v=zeros(ComplexF64, num_modes, num_modes), plasma_pts=zeros(num_points, 3), wall_pts=zeros(num_points, 3) ) - compute_vacuum_response!(vac, inputs, wall_settings) + compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv) - return vac.wv, vac.grri, vac.grre, vac.plasma_pts, vac.wall_pts + return vac.wv, vac.I_v, vac.plasma_pts, vac.wall_pts end """ - compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) In-place variant that computes the vacuum response and directly populates the arrays stored in `vac_data`. Dispatches on dimensionality only: 2D (`inputs.nzeta == 1`) routes to @@ -299,19 +306,18 @@ The `vac_data` argument is expected to provide the following writable fields wit sizes: - `wv::AbstractMatrix{ComplexF64}` – vacuum response matrix - - `grri::AbstractMatrix{ComplexF64}` – interior Green's functions - - `grre::AbstractMatrix{ComplexF64}` – exterior Green's functions - `plasma_pts::AbstractMatrix{Float64}` – plasma surface coordinates - `wall_pts::AbstractMatrix{Float64}` – wall surface coordinates + - `I_v::AbstractMatrix{ComplexF64}` – required when `compute_Iv=true` This is designed to work with `ForceFreeStates.VacuumData` but does not depend on its concrete type (duck-typed on field names only). """ -function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) if inputs.nzeta == 1 - _compute_vacuum_response_2d!(vac_data, inputs, wall_settings) + _compute_vacuum_response_2d!(vac_data, inputs, wall_settings; compute_Iv) else - _compute_vacuum_response_3d!(vac_data, inputs, wall_settings) + _compute_vacuum_response_3d!(vac_data, inputs, wall_settings; compute_Iv) end end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 7e8b51b46..92cc8afdb 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -368,7 +368,7 @@ @testset "nowall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -376,10 +376,8 @@ @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(grri) == (2 * numpoints, num_modes) - @test size(grre) == (2 * numpoints, num_modes) - @test all(isfinite, grri) - @test all(isfinite, grre) + @test size(I_v) == (num_modes, num_modes) + @test all(iszero, I_v) # compute_Iv=false by default @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -388,15 +386,29 @@ @test isapprox(wv, wv', rtol=1e-12) end + @testset "nowall compute_Iv=true" begin + inputs = _make_inputs() + wall_settings = WallShapeSettings(shape="nowall") + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) + + num_modes = length(inputs.m_modes) * length(inputs.n_modes) + @test size(wv) == (num_modes, num_modes) + @test all(isfinite, wv) + @test size(I_v) == (num_modes, num_modes) + @test all(isfinite, I_v) + @test !all(iszero, I_v) + @test isapprox(wv, wv', rtol=1e-12) + end + @testset "conformal wall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="conformal", a=0.5) - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test size(grri) == (2 * numpoints, num_modes) + @test all(iszero, I_v) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # plasma_pts layout: col1=R, col2=0, col3=Z @@ -408,38 +420,40 @@ @testset "edge: single poloidal mode mpert=1" begin inputs = _make_inputs(m_modes=[1], n_modes=[1]) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (1, 1) @test all(isfinite, wv) - @test size(grri, 2) == 1 + @test size(I_v) == (1, 1) end @testset "edge: small mtheta" begin # Keep mtheta_eq=17 so boundary has enough points for periodic spline inputs = _make_inputs(mtheta=16, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (2, 2) - @test size(grri) == (32, 2) # 2*mtheta, num_modes=2 + @test size(I_v) == (2, 2) @test size(plasma_pts) == (16, 3) end @testset "in-place compute_vacuum_response! matches wrapper" begin # The allocating wrapper is a thin caller of the in-place routine; verify the - # in-place entry populates caller-owned duck-typed (NamedTuple) storage identically. + # in-place entry populates caller-owned storage identically. Both a duck-typed + # NamedTuple and the real ForceFreeStates.VacuumData are exercised, so the + # struct is kept in sync with the documented field contract. for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.5)) inputs = _make_inputs() - wv, grri, grre, pp, wp = compute_vacuum_response(inputs, wall_settings) + wv, I_v, pp, wp = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) - vac = (wv=zeros(ComplexF64, num_modes, num_modes), grri=zeros(ComplexF64, 2 * numpoints, num_modes), grre=zeros(ComplexF64, 2 * numpoints, num_modes), + vac = (wv=zeros(ComplexF64, num_modes, num_modes), + I_v=zeros(ComplexF64, num_modes, num_modes), plasma_pts=zeros(numpoints, 3), wall_pts=zeros(numpoints, 3)) - compute_vacuum_response!(vac, inputs, wall_settings) + compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv=true) @test vac.wv ≈ wv - @test vac.grri ≈ grri - @test vac.grre ≈ grre + @test vac.I_v ≈ I_v @test vac.plasma_pts ≈ pp @test vac.wall_pts ≈ wp end @@ -581,17 +595,15 @@ @testset "compute_vacuum_response 3D nowall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(grri) == (2 * numpoints, num_modes) - @test size(grre) == (2 * numpoints, num_modes) - @test all(isfinite, grri) - @test all(isfinite, grre) + @test size(I_v) == (num_modes, num_modes) + @test all(iszero, I_v) @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -604,7 +616,7 @@ @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -612,10 +624,8 @@ @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(grri) == (2 * numpoints, num_modes) - @test size(grre) == (2 * numpoints, num_modes) - @test all(isfinite, grri) - @test all(isfinite, grre) + @test size(I_v) == (num_modes, num_modes) + @test all(iszero, I_v) @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -625,12 +635,12 @@ @testset "compute_vacuum_response 3D conformal wall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="conformal", a=0.3) - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test size(grri) == (2 * numpoints, num_modes) + @test size(I_v) == (num_modes, num_modes) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # Wall and plasma should differ (conformal wall offset from plasma) @@ -676,13 +686,13 @@ wall_settings = WallShapeSettings(shape="nowall") # Reduced (block-circulant) path - wv_red, _, _, plasma_pts_red, _ = compute_vacuum_response(inputs_red, wall_settings) + wv_red, _, plasma_pts_red, _ = compute_vacuum_response(inputs_red, wall_settings) # Full-torus reference: pre-expand so nfp=1 forces the dense path inputs_full = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(inputs_red) @test inputs_full.nfp == 1 @test inputs_full.nzeta == nzeta_p * nfp - wv_full, _, _, _, _ = compute_vacuum_response(inputs_full, wall_settings) + wv_full, _, _, _ = compute_vacuum_response(inputs_full, wall_settings) num_modes = length(m_modes) * length(n_modes) @test size(wv_red) == (num_modes, num_modes)