ForcingTerms Module

BoundaryGrid

Pre-computed plasma boundary grid for evaluating coil fields.

Fields

• mtheta, nzeta: grid dimensions
• R, Z: cylindrical coordinates [mtheta] (same for all ζ, axisymmetric)
• phi_grid: base toroidal angle grid [nzeta] in radians = -helicity × 2π × j/nzeta
• phi_offset: per-θ toroidal angle correction [mtheta] from SFL coordinates: ν(ψ, θ) scaled by -helicity, so the physical toroidal angle at (i, j) is phi_grid[j] + phi_offset[i]. Matches Fortran's phi = -helicity*(2π*ζ + dphi(ψ,θ)). Zero for axisymmetric equilibria on-axis; non-zero off-axis due to SFL coordinate transform (Hamada/SFL θ ≠ geometric θ introduces a toroidal offset).
• dR_dtheta, dZ_dtheta: poloidal derivatives w.r.t. unit-norm angle θnorm ∈ [0,1] `dRdtheta[i] = dR/dθnorm = 2π × dR/dθphys`

CoilConfig

Top-level coil configuration from the [ForcingTerms] TOML section when forcing_data_format = "coil".

Fields

• machine: machine name prefix for .dat files (e.g. "d3d")
• dat_dir: directory containing .dat files; defaults to bundled coil_geometries/
• mtheta_coil: poloidal grid resolution for boundary field evaluation (default: 480)
• nzeta_coil: toroidal grid resolution; 0 = auto (32 × n)
• coil_sets: parsed [[ForcingTerms.coil_set]] blocks

CoilConfig(ft_ctrl::ForcingTermsControl)

Construct a CoilConfig from a ForcingTermsControl populated from TOML.

CoilSet

Geometry and current data for a single coil set (one .dat file).

Fields

• name: coil set identifier (e.g. "il", "iu")
• ncoil: number of independent conductors
• s: sub-systems (strand groups) per conductor
• nw: winding multiplier (turns per conductor element)
• nsec: number of Cartesian points per strand
• x, y, z: conductor coordinates [ncoil, s, nsec] in meters (Cartesian)
• currents: current per conductor [ncoil] in Amperes

CoilSetConfig

Per-coil-set configuration from a [[ForcingTerms.coil_set]] TOML block.

A coil set comes from one of several sources, selected by source:

• "" or "file": read a geometry file — .dat (legacy ASCII) or .h5 (modern), chosen by extension. The file path comes from dat_file/h5_file or the machine+name convention.
• "pf_hoop": an analytic horizontal circular loop (radius, height).
• "window_pane": an analytic rectangular picture-frame array, placed EITHER by explicit corners (rz_corners) OR by physical standoff from the plasma surface (standoff, poloidal_angle, poloidal_length, poloidal_tilt); shared toroidal layout via ncoil_gen, gap_fraction, phi0.
• "helical": an analytic helical array on a circular cross-section torus (R0, a, m_hel, n_hel, theta_lo, theta_hi, n_coils, gap_fraction, phi0, offset).

Shifts/tilts apply to every source (the analytic geometry is built first, then transformed).

Fields

• name: coil set name; resolves to {dat_dir}/{machine}_{name}.dat, or selects a set by name from a multi-set .h5 file
• source: coil source discriminator (see above); empty/"file" reads a geometry file
• dat_file: explicit path to .dat file (overrides machine+name convention)
• h5_file: explicit path to .h5 file (multi-set; name selects which set, else the sole set)
• currents: current [A] per conductor [ncoil]; shorter arrays pad with zeros
• shiftx, shifty, shiftz: per-conductor translation [m] [ncoil]
• tiltx, tilty, tiltz: per-conductor tilt in degrees (or meters if tilt_in_meters)
• xnom, ynom, znom: explicit rotation center [m]; defaults to arc-length-weighted center of mass
• n_tilt: toroidal mode number for tilt/shift modulation; -1 means inherit run's n
• tilt_in_meters: interpret tilt as displacement [m] instead of angle [degrees]

Analytic-source parameters (used only for the matching source):

• radius, height: PF hoop circle radius and Z height [m]
• ncoil_gen: number of generated coils (window-pane / helical arrays)
• rz_corners: window-pane cross-section as two opposite [R, Z] corners in m (corner mode)
• gap_fraction: toroidal gap between adjacent coils as a fraction of their spacing
• phi0: toroidal angle offset of the first generated coil [rad]
• R0, a: helical torus major radius and minor (circle) radius [m]
• m_hel, n_hel: helical poloidal/toroidal mode numbers setting the pitch
• theta_lo, theta_hi: helical poloidal extent in degrees; a full turn when theta_hi-theta_lo == 360
• n_coils: number of helical coils in the array
• offset: helical coil radial offset from the torus surface [m]
• nsec_gen: optional override of the generated point count per coil

Window-pane standoff placement (alternative to rz_corners; needs an equilibrium):

• standoff: outward-normal distance from the control surface to the frame center [m]
• poloidal_angle: cylindrical poloidal angle of the frame center, in degrees
• poloidal_length: poloidal leg span (corner-to-corner) along the surface tangent [m]
• poloidal_tilt: extra tilt of the legs in degrees; 0 means tangential to the surface

ForcingMode

Data structure for a single forcing mode.

Fields

• n::Int - Toroidal mode number

• m::Int - Poloidal mode number

• amplitude::ComplexF64 - Complex amplitude in unit-norm convention (= Fortran Phi_x, T·m² per unit-norm cell) after loading. File inputs are tagged by their input convention:

  • "normal_field_T": B·n̂ in Tesla (2π-angle); converted to unit-norm on load
  • "sfl_flux_Wb": SFL flux in 2π-angle convention; multiplied by (2π)² on load

ForcingTermsControl

User-facing control parameters from TOML [ForcingTerms] section.

Fields

Forcing Data:

• forcing_data_file::String - Path to forcing data file (n, m, complex amplitude)
• forcing_data_format::String - Format: "ascii", "hdf5", or "coil"

Coil settings (used when forcing_data_format = "coil"):

• machine::String - Machine name prefix for .dat files (e.g. "d3d")
• dat_dir::String - Directory containing .dat files; defaults to bundled coil_geometries/
• mtheta_coil::Int - Poloidal grid resolution for boundary field evaluation (default: 480)
• nzeta_coil::Int - Toroidal grid resolution; 0 = auto (32 × n)
• coil_sets_raw::Vector{Dict{String,Any}} - Parsed [[ForcingTerms.coil_set]] TOML blocks

PrecomputedSegments

Pre-computed segment vectors and midpoints for one conductor strand, stored in contiguous arrays for cache-friendly access in the inner Biot-Savart loop.

_accumulate_precomputed(ox, oy, oz, seg) -> (dbx, dby, dbz)

Inner kernel: accumulate Biot-Savart field from pre-computed segments. Returns (dbx, dby, dbz) increments. Matches Fortran field_bs_psi midpoint rule.

Operates on PrecomputedSegments with contiguous arrays for maximum cache efficiency.

_arc_length_center(x, y, z) -> (x0, y0, z0)

Compute the arc-length-weighted centroid of a strand (1D array of points). Matches Fortran center-of-mass convention using midpoints weighted by segment length.

_build_raw_coil_set(cfg::CoilConfig, sc::CoilSetConfig) -> CoilSet

Produce the un-transformed coil set for a config block, dispatching on sc.source: analytic generators (pf_hoop/window_pane/helical) or a geometry file (.h5 → read_coil_h5, else .dat → read_coil_dat). Currents and shifts/tilts are applied by the caller, so all sources share the same downstream path.

_resolve_coil_path(cfg::CoilConfig, set_cfg::CoilSetConfig) -> String

Resolve the path to a coil geometry file. Prefers set_cfg.h5_file, then set_cfg.dat_file, otherwise constructs {dat_dir}/{machine}_{name}.dat.

accumulate_strand_field!(Bx, By, Bz, i, obs_x, obs_y, obs_z, xs, ys, zs, current_nw)

Accumulate the Biot-Savart contribution from one conductor strand onto observation point i.

Uses midpoint-rule integration: each segment from point l to l+1 contributes: dB = μ₀/(4π) × current_nw × (dl⃗ × r⃗) / |r|³ where r⃗ is the displacement from the segment midpoint to the observer.

Arguments

• Bx, By, Bz: pre-allocated output arrays (Cartesian field, Tesla); updated in-place
• i: observation point index
• obs_x/y/z: observer Cartesian coordinates [m]
• xs, ys, zs: strand point coordinates [nsec], Cartesian [m]
• current_nw: effective current = I × nw [A⋅turns]

apply_transforms(cs::CoilSet, cfg::CoilSetConfig; n_tilt::Int=1) -> CoilSet

Apply per-conductor shifts and tilts to a coil set, returning a modified copy.

Replicates the Fortran coil_read shift/tilt logic (coil.F lines 240–340):

• Tilts are rotations around the arc-length-weighted center of mass (unless xnom/ynom/znom specified)
• n_tilt controls the toroidal periodicity of tilt/shift modulation
• n_tilt = 0: rigid shift only (no tilts applied)
• n_tilt ≥ 1: n-fold modulated perturbations

compute_biot_savart_boundary!(B_R, B_phi, B_Z, obs_R, obs_phi, obs_Z, coil_sets)

Compute the total Biot-Savart field from all coil sets at every observation point.

Observation points are specified in cylindrical coordinates (R, φ, Z) and the returned field components are also cylindrical.

Uses Threads.@threads over observation points – each thread writes to a unique index so no synchronization is required.

Performance: segment vectors and midpoints are pre-computed once outside the observation-point loop, then accessed from contiguous arrays in the inner kernel.

Arguments

• B_R, B_phi, B_Z: output field arrays [Tesla], length nobs; overwritten
• obs_R, obs_phi, obs_Z: observation point cylindrical coordinates [m, rad, m]
• coil_sets: conductor geometry and current data

compute_coil_forcing_modes!(forcing_modes, coil_sets, equil, cfg, n, m_low, m_high; verbose)

Top-level entry point: compute Fourier mode amplitudes of the normal magnetic flux from all coil sets on the plasma boundary.

Pipeline

• Build boundary grid at (mtheta × nzeta) resolution, with helicity from equil.params
• Lay out observation points in (R, φ, Z) for all (θ, ζ) combinations
• Run threaded Biot-Savart summation over all coil sets (matches Fortran field_bs_psi)
• Project B field onto plasma boundary normal flux (project_normal_flux!)
• 2D Fourier decompose to get bmn amplitudes for mode range

Output amplitudes are in unit-norm convention (= Fortran Phi_x). No normalization conversion is needed when using these modes with compute_plasma_response!.

Result is appended to forcing_modes (existing content is cleared first).

convert_coil_dat_to_h5(dat_path, h5_path; group="coils") -> String

Convert a legacy ASCII .dat coil file to the modern HDF5 format (single set).

convert_coil_h5_to_dat(h5_path, dat_path; group="coils", set=nothing) -> String

Convert one coil set from an HDF5 file to a legacy ASCII .dat file. If the file holds multiple sets, set selects which one by name (required when more than one is present).

convert_forcing_normalization!(modes, from_tag, equil, n, m_low, m_high; psi, mtheta, nzeta)

Convert ForcingMode amplitudes from from_tag normalization to the internal unit-norm convention (= Fortran Phi_x), expected by compute_plasma_response!.

If from_tag == "sfl_flux_Wb", the user has provided amplitudes in the 2π-angle SFL-flux convention. These are scaled by (2π)² to reach unit-norm.

If from_tag == "normal_field_T", the amplitudes represent Fourier modes of B·n̂ in Tesla (2π-angle convention). Conversion to unit-norm Phi_x:

• Inverse-Fourier reconstruct B·n̂(θ, ζ) from the input modes
• Multiply pointwise by 2π × R(θ) × |dr/dθ_norm(θ)|
• Re-Fourier transform to get unit-norm mode amplitudes

The 2π factor comes from the toroidal Jacobian ∂r/∂ζnorm = 2π·R·êφ. This conversion is a mode-mixing operation because R and |dr/dθ| vary poloidally.

Arguments

• modes: Vector{ForcingMode} with amplitudes to convert (modified in place)
• from_tag: source normalization — "normal_field_T" or "sfl_flux_Wb"
• equil: PlasmaEquilibrium providing boundary geometry
• n: toroidal mode number
• m_low, m_high: poloidal mode range (must cover all modes in modes)
• psi: flux surface for boundary geometry (default 1.0)
• mtheta, nzeta: grid resolution for the conversion (defaults: 256, 64)

fourier_decompose_bn(bn, grid, n, m_low, m_high) -> Vector{ForcingMode}

2D Fourier decompose bn[mtheta, nzeta] to extract mode amplitudes for toroidal mode number n and poloidal range m_low:m_high.

Coefficients are computed as: bmn = (2 / (mtheta × nzeta)) × Σ_{i,j} bn[i,j] × exp(-i(m×θᵢ - n×ζⱼ))

The factor 2 matches the GPEC/DCON convention for real signals where positive and negative m modes are related by conjugation.

When called after project_normal_flux!, the returned amplitudes are in unit-norm convention equal to Fortran Phi_x (T·m² per unit-norm cell).

Uses compute_fourier_coefficients from Utilities.FourierTransforms with the 3D (mtheta×nzeta, mpert) basis matrix (npert=1, nlow=n).

load_coil_sets(cfg::CoilConfig, n_tilt::Int; equil=nothing) -> Vector{CoilSet}

Build all coil sets from config (file or analytic source), apply currents and shifts/tilts, and return ready-to-use coil data. n_tilt is the toroidal mode number for modulated perturbations (typically the run's n). equil is required only when a window_pane set uses standoff placement (it locates the surface point and normal).

load_coils_from_h5_group!(coil_sets::Vector{CoilSet}, group) -> Vector{CoilSet}

Populate coil_sets from an already-open HDF5 group, reading every subgroup that holds x/y/z datasets as a CoilSet. Any other datasets or subgroups under group are silently ignored, so the group may live inside a larger file (e.g. gpec.h5). Mirror of save_coils_to_h5.

load_forcing_ascii!(forcing_modes, filepath, verbose) -> String

Load forcing data from ASCII file. Returns normalization tag ("normal_field_T" by default).

Optional first # normalization: <tag> header line sets the normalization. Remaining #-prefixed lines and blank lines are ignored. Data columns: n, m, amplitudereal, amplitudeimag.

load_forcing_data!(forcing_modes, dir_path, forcing_data_file, forcing_data_format, verbose) -> String

Load forcing data from ASCII or HDF5 file. Returns the normalization tag read from the file (or "normal_field_T" by default if no tag is present).

Normalization tags (set in the file, not in the TOML):

• "normal_field_T" (default): Fourier modes of B·n̂ [Tesla], 2π-angle convention. Most intuitive for users; converted to unit-norm (Phi_x) on load.
• "sfl_flux_Wb": SFL flux R×(BR·∂Z/∂θ - BZ·∂R/∂θ) [T·m²], 2π-angle convention. Multiplied by (2π)² on load to reach unit-norm (Phi_x).

ASCII format — optional # normalization: <tag> header line, then data rows:

# normalization: normal_field_T
1  2  0.5  0.1
1  3  0.3  -0.2

Columns: n, m, amplitudereal, amplitudeimag

HDF5 format — optional root attribute normalization (string), plus datasets:

• "n": Integer array of toroidal mode numbers
• "m": Integer array of poloidal mode numbers
• "amplitude_real": Real parts of amplitudes
• "amplitude_imag": Imaginary parts of amplitudes

load_forcing_from_h5_group!(forcing_modes::Vector{ForcingMode}, group)

Populate forcing_modes from an already-open HDF5 group with datasets n, m, amplitude_real, amplitude_imag. Mirror of save_forcing_to_h5 for the rerun path — deliberately accepts an open group rather than a file path so the snapshot data can live inside gpec.h5/input/raw_inputs/forcing_terms/.

load_forcing_hdf5!(forcing_modes, filepath, verbose) -> String

Load forcing data from HDF5 file. Returns normalization tag ("normal_field_T" by default).

Optional root attribute normalization (string) sets the normalization. Required datasets: n, m, amplitude_real, amplitude_imag.

make_helical(; R0, a, m, n, theta_lo=0.0, theta_hi=360.0, n_coils,
               gap_fraction=0.0, phi0=0.0, offset=0.0, npts=720,
               name="helical") -> CoilSet

Build an array of n_coils analytic helical coils wound on a circular cross-section torus of major radius R0 and minor radius a [m], offset radially outward by offset [m]. The helical pitch is set by the poloidal/toroidal mode numbers m/n via ζ = (m/n)·θ + φ.

The poloidal extent theta_lo/theta_hi is in degrees. For a full poloidal turn (theta_hi - theta_lo == 360) each coil is a single helical strand. For a partial poloidal range the coil is closed into a loop with two helical legs joined by toroidal arcs (matching OMFIT form_helical_array.py). npts sets the strand resolution. Currents are zeroed; load_coil_sets fills them.

make_pf_hoop(; radius, height, nsec=361, name="pf_hoop") -> CoilSet

Build a single horizontal circular loop (PF hoop) of radius [m] centered on the machine axis at vertical position height [m], discretized into nsec points.

The loop is closed (point nsec coincides with point 1). A bare hoop is axisymmetric and produces only an n=0 field; tilting or shifting it (via apply_transforms) breaks axisymmetry to produce n≠0 forcing. Currents are zeroed; load_coil_sets fills them.

make_window_pane(; ncoil, rz_corners, gap_fraction=0.0, phi0=0.0,
                   np_pol=20, np_tor=40, name="window_pane") -> CoilSet

Build an array of ncoil rectangular "picture-frame" (window-pane) saddle coils, evenly spaced toroidally. The two [R, Z] points in rz_corners are the endpoints of the poloidal legs; the legs run along the straight segment between them at the two toroidal edges of the coil, joined by toroidal arcs at each corner. For a vertical frame give the two corners the same R (legs run in Z); for a horizontal frame give them the same Z. Adjacent coils are separated by a toroidal gap of gap_fraction times their angular spacing; the first coil is centered at toroidal angle phi0 [rad].

Each coil is a closed loop with np_pol points per poloidal leg and np_tor per toroidal arc. Currents are zeroed; load_coil_sets fills them.

make_window_pane_standoff(equil; standoff, poloidal_angle, poloidal_length,
                          poloidal_tilt=0.0, ncoil, gap_fraction=0.0, phi0=0.0,
                          name="window_pane") -> CoilSet

Build a window-pane array placed by physical standoff instead of explicit corners. The frame center sits standoff [m] along the outward surface normal at cylindrical poloidal_angle (degrees); its poloidal legs span poloidal_length [m] along the surface tangent, rotated by poloidal_tilt (degrees, 0 = tangential). The control surface is the boundary equil.rzphi_xs[end]. Toroidal layout (ncoil, gap_fraction, phi0) and the leg/arc construction are delegated to make_window_pane.

The poloidal legs wind in the same sense as make_window_pane's c1 -> c2 ordering (lower-Z to higher-Z at the outboard midplane), so positive current produces positive b_n_x near the coil, matching the corner-built window pane.

precompute_segments(cs, j, k, current_nw) → PrecomputedSegments

Build contiguous segment-vector and midpoint arrays for strand (j,k) of coil set cs.

project_normal_flux!(bn, B_R, B_Z, grid)

Project the cylindrical magnetic field (BR, BZ) onto the plasma boundary normal direction ∇ψ and store in bn[mtheta, nzeta].

Computes the unit-norm flux element Phix per (θnorm, ζnorm) cell [T·m²]: bn(θnorm, ζnorm) = 2π × R(θnorm) × (BR × ∂Z/∂θnorm - BZ × ∂R/∂θnorm)

The 2π factor comes from the toroidal Jacobian ∂r/∂ζnorm = 2π·R·êφ in the cross-product ∂r/∂θnorm × ∂r/∂ζnorm. The derivatives dR/dθnorm and dZ/dθnorm are stored in grid.dR_dtheta and grid.dZ_dtheta (unit-norm convention from sample_boundary_grid).

The output matches Fortran GPEC's Phi_x convention directly (no extra factor needed).

Arguments

• bn: output array [mtheta, nzeta]; overwritten in-place
• B_R, B_Z: cylindrical field components, length mtheta × nzeta (flat, θ-major)
• grid: pre-computed boundary geometry from sample_boundary_grid

read_coil_dat(filepath::String) -> CoilSet

Parse an ASCII coil geometry file.

Format:

ncoil s nsec nw
X₁ Y₁ Z₁
X₂ Y₂ Z₂
...

Data rows are ordered: for each conductor 1:ncoil, for each sub-system 1:s, read nsec points. Coordinates are Cartesian in meters (lab frame).

read_coil_h5(filepath; group="coils", set=nothing) -> Vector{CoilSet}

Read coil sets from an HDF5 file. Coil data is expected under the group path; all coil-set subgroups there are returned (or, if set names one, just that set). Unrelated content elsewhere in the file is ignored. See save_coils_to_h5 for the layout.

sample_boundary_grid(equil, mtheta, nzeta; psi=1.0) -> BoundaryGrid

Evaluate plasma geometry at a uniform (mtheta × nzeta) grid on the flux surface psi.

Uses equil.rzphi_rsquared and equil.rzphi_offset splines at the given psi. Defaults to psi=1.0 (true plasma boundary). Use a smaller value (e.g. the Fortran psilim) to evaluate on an interior truncation surface. Poloidal derivatives dR/dθnorm and dZ/dθnorm (unit-norm θnorm ∈ [0,1]) are computed via periodic cubic splines on the resulting R(θnorm), Z(θ_norm) data.

The toroidal grid direction follows the Fortran GPEC convention: phi_j = -helicity × 2π × j/nzeta where helicity = sign(Bt) × sign(Ip). This is derived from equil.params.bt_sign and equil.params.crnt. For DIII-D (Bt < 0, Ip > 0 → helicity = -1): phi increases with j (standard direction). For positive-helicity machines (Bt > 0, Ip > 0 → helicity = +1): phi decreases with j.

save_coils_to_h5(coil_sets::Vector{CoilSet}, group)

Write each coil set as a subgroup (named by CoilSet.name) under an already-open HDF5 group. Inverse of load_coils_from_h5_group!. Used by the gpec.h5 snapshot writer so a coil run can be replayed from the output file alone.

save_forcing_to_h5(forcing_modes::Vector{ForcingMode}, group)

Write forcing modes to an open HDF5 group using the same layout that load_forcing_hdf5! consumes: datasets n, m, amplitude_real, amplitude_imag. Used by the gpec.h5 snapshot writer so the replay path can re-read forcing data from the output file with no reference to the original ASCII/HDF5 ingest file.

surface_point_and_normal(equil, theta_cyl; psi=equil.rzphi_xs[end]) -> (R, Z, nR, nZ)

Locate the control-surface point at cylindrical poloidal angle theta_cyl (radians) and return it with the outward unit normal (nR, nZ) in the R-Z plane.

The flux-surface splines are parametrized by the SFL angle theta_sfl, related to the cylindrical angle by theta_cyl = 2π*(theta_sfl + rzphi_offset(psi, theta_sfl)). Since the offset is a small smooth perturbation, theta_sfl is recovered by a damped fixed-point iteration. The tangent is taken by central-differencing R(theta_sfl), Z(theta_sfl), and the normal (dZ/dθ, -dR/dθ) is normalized and oriented outward (positive dot with (R-ro, Z-zo)), matching the normal-flux convention of project_normal_flux!.

write_coil_dat(filepath, cs::CoilSet)

Write a single coil set to an ASCII .dat file in the legacy Fortran GPEC format (ncoil s nsec nw header, then ncoil×s×nsec rows of X Y Z).

write_coil_h5(filepath, coil_sets; group="coils")

Write coil_sets to a new HDF5 file under the group path, one subgroup per set. Convenience file-level wrapper around save_coils_to_h5 for export and the .dat→.h5 converter.