[rogeriojorge]

This updates the winding-surface optimization example and adds a resolved comparison with REGCOIL and QUADCOIL.

Corrections

• include winding-surface Jacobians and quadrature weights in the induction matrix
• use close=False to exclude duplicate periodic endpoints
• include the missing 4*pi^2 normalization
• prevent the entropy objective from improving by weakening the full operator
• add a nonlocal self-separation check alongside the area-element check
• use 64-bit JAX calculations

ESSOS peak-current Pareto objective

For every trial winding surface, ESSOS solves the regularized current-potential problem at three nearby peak-current limits:

$$ \Phi_i=\arg\min_{\Phi}\left[f_B(\Phi)+\lambda_i f_K(\Phi)\right], \qquad K_{\max}(\Phi_i)=K_i, $$

where

$$ f_B=\int_{S_p}(\mathbf B\cdot\mathbf n)^2,dA, \qquad f_K=\int_{S_c}|\mathbf K|^2,dA =\int_{S_c}|\nabla_s\Phi|^2,dA, $$

and $K_i=(0.95,1.00,1.05)K_{\mathrm{target}}$. The winding surface minimizes

$$ J_{\mathrm{Pareto}}(S_c)=\frac{1}{3}\sum_{i=1}^{3} \frac{f_B(S_c;K_i)}{f_B(S_{c,0};K_i)}+J_{\mathrm{geometry}}. $$

The three samples cover a short section of the trade-off between $f_B$ and $K_{\max}$ and reduce sensitivity to one regularization value. The geometry term preserves volume, controls the Fourier spectrum, retains at least 90% of the initial plasma clearance, and penalizes small Jacobian or nonlocal self-approach. The result is a local three-point approximation to the Pareto curve.

The optimization uses 32 points per angle per field period. Surface motion is capped at 2% of the plasma minor radius per active Fourier coefficient.

Results

For the included QA example at 7 MA/m, a separate 96x96 REGCOIL evaluation gives:

• sheet fB: 19.873 -> 2.712
• max |Bn|/B: 0.0569 -> 0.0267
• 40-filament fB: 20.331 -> 3.010
• dense Kmax: 6.989 MA/m

At 10 MA/m, the Pareto surface reduces resolved sheet fB by 27.6% relative to the normal offset for QA and by 88.2% for W7-X. W7-X filament fB is reduced by 88.1%. The QA filament result worsens despite its improved sheet result, which limits the coil-improvement result to W7-X in this comparison.

The current-complexity comparison uses the same $f_K$ as the inner regularization. Relative to the normal offset, ESSOS Pareto changes $f_K$ by +0.24% for QA and +22.2% for W7-X. W7-X RMS current increases from 2.208 to 2.432 MA/m while both surfaces satisfy the same 10 MA/m peak limit. The field-error improvement therefore has a visible current-complexity cost.

The driver, numerical tables, all figures, and detailed notes are under examples/winding_surface_comparison/. The W7-X QUADCOIL solves are labeled failed after exceeding the convergence limit.

Eight-coil study

An additive study under examples/winding_surface_comparison_8_coils/ repeats the comparison with eight coils per half field period and adds Landreman-Paul QH. At 96x96, the QH Pareto iterate reduces sheet fB by 75.1% and filament fB by 73.6% relative to the normal offset. The QA case shows the limitation clearly: its sheet result improves, but filament fB is 15.7% worse with the reduced coil count. The QH Pareto outer solve stopped on a line-search failure, so it is reported as a valid non-converged iterate.

[rogeriojorge]

@dseidita1919 Doing my own tests, to see if a couple things might help:

• Adding winding and plasma Jacobian/quadrature weights as Landreman & Boozer mention
• Explicit close=False, since ESSOS’s close=True includes a duplicate periodic endpoint (perhaps we should change this, as SIMSOPT does)
• Missing 4 pi^2 full-torus integration factor for minor radius and volume
• float64 calculations
• A tangent-point self-intersection check based on https://arxiv.org/abs/1102.3642 and https://arxiv.org/abs/2107.01664
• There seems to be a loophole where every singular value decreases while the spectrum flattens, giving higher entropy with a weaker operator

[rogeriojorge]

The resolved comparison is now included in the PR under examples/winding_surface_comparison.

Main findings:

• At 10 MA/m, the ESSOS Pareto surface reduces the 96x96 sheet fB by 27.6% for Landreman-Paul QA and 88.2% for W7-X relative to a normal-offset surface.
• The W7-X filament result improves by 88.1% (5.909 to 0.703), with max |Bn|/B reduced from 0.260 to 0.109.
• The QA sheet improves, but its filament fB becomes worse. This is an important limitation and is reported explicitly.
• The entropy proxy is worse than the normal-offset baseline in both resolved sheet comparisons.
• ESSOS and REGCOIL current solves agree when evaluated at the same resolution. QUADCOIL agrees for QA after refinement; its W7-X solves are labeled failed after exceeding the convergence limit.
• The standalone 7 MA/m example reduces validated sheet fB from 19.873 to 2.712. Its saved coefficients have dense Kmax = 6.989 MA/m.

Resolved sheet and filament metrics:

Resolved current potentials and coil contours:

Resolved winding surfaces and coils:

The remaining figures and CSV tables are in the same folder.

[rogeriojorge]

Added the requested Pareto definition and current-complexity comparison.

For each trial surface, the inner problems are

$$ \Phi_i=\arg\min_{\Phi}\left[f_B(\Phi)+\lambda_i f_K(\Phi)\right], \qquad K_{\max}(\Phi_i)=K_i, $$

with $K_i=(0.95,1.00,1.05)K_{\mathrm{target}}$ and

$$ f_K=\int_{S_c}|\mathbf K|^2,dA =\int_{S_c}|\nabla_s\Phi|^2,dA. $$

The outer objective minimizes the mean normalized $f_B$ at those three current limits, plus the explicit geometry controls. Sampling three nearby points improves a local part of the trade-off between $f_B$ and $K_{\max}$ and reduces sensitivity to the selected regularization value.

The new figure compares $f_K/f_{K,\mathrm{offset}}$ and $K_{\mathrm{rms}}=\sqrt{f_K/A_c}$ for all four winding-surface methods:

For QA, the Pareto surface has essentially the same complexity as the normal offset: $f_K$ changes by +0.24% and RMS current changes from 4.956 to 4.953 MA/m. For W7-X, $f_K$ increases by 22.2% and RMS current increases from 2.208 to 2.432 MA/m. Both still satisfy the same 10 MA/m peak-current limit. This makes the cost of the W7-X field-error improvement explicit.

The dimensional $f_K$ and RMS-current values are included in data/surface_validation_96.csv.

[rogeriojorge]

Added a second, additive study with eight coils per half field period. The original ten-coil results are unchanged. The new study also adds Landreman-Paul QH.

Main results from the independent 96x96 validation:

• QH: ESSOS Pareto reduces sheet fB by 75.1% and filament fB by 73.6% relative to the normal offset. Relative to the tested REGCOIL adjoint surface, the reductions are 65.1% and 63.1%.
• W7-X: filament fB remains 88.1% below the normal-offset result.
• QA: sheet fB improves by 27.6%, but filament fB is 15.7% worse. Going from ten to eight coils per half period raises QA filament error by about a factor of eight, showing the limitation of the present sheet objective for sparse coil sets.
• The QH Pareto run stopped after seven outer iterations on a line-search failure. Its saved surface is valid and feasible at 10 MA/m on the 96x96 grid; the optimizer status is non-converged. QUADCOIL exceeded six minutes and about 2.8 GB on QH, and those entries are labeled failed.

Resolved metrics:

Current potential and coil contours:

Winding surfaces and coils:

The reproducible driver, CSV tables, complexity plot, convergence plot, and remaining figures are in examples/winding_surface_comparison_8_coils/.

[rogeriojorge]

@dseidita1919 I tested a simpler entropy formulation on QA, QH and QI.

The full-torus induction matrix contains several Bloch sectors, but a field-periodic plasma and repeated coil set use only the symmetric sector:

$$ A_{\mathrm{sym}}=\sum_{k=0}^{N_{\mathrm{fp}}-1}A_k, \qquad p_i=\frac{\sigma_i(A_{\mathrm{sym}})}{\sum_j\sigma_j(A_{\mathrm{sym}})}, \qquad H_{\mathrm{sym}}=-\sum_i p_i\log p_i. $$

This changes the calculation from $N_{\mathrm{fp}}$ SVDs to one physically relevant SVD. It reproduced the full-entropy cut-coil results while making the surface optimization 2.3–4.7 times faster.

With eight coils per half field period, relative to the normal-offset surface:

Case	max $\lvert B_n\rvert/B$	filament $f_B$	mean $\int\kappa^2dl/L$
QA	−10.5%	−4.7%	+9.8%
QH	−13.5%	−25.4%	+3.7%
QI	−3.9%	+3.7%	−3.3%

A simple stage-1 use would be:

1. Build the control surface at a fixed normalized offset $d/a$ for every trial plasma.
2. Set a target for the required symmetric spectral dimension:

$$ J_{\mathrm{SVD}}=\left[\max(0,H_{\mathrm{target}}-H_{\mathrm{sym}})\right]^2. $$

Fixing $d/a$ removes the global operator-strength degree of freedom, so the separate SINGULAR_STRENGTH_WEIGHT can be removed. For a freely moving winding surface, a hard distance or $\sum_i\sigma_i$ constraint is still needed before removing it.

The main gain is a smaller and faster objective with consistent improvement in cut-coil maximum $\lvert B_n\rvert/B$. Filament $f_B$ and coil complexity improve inconsistently. The present entropy-plus-strength objective also gave a larger QH maximum-error reduction: 21.5% compared with 13.5%.

[rogeriojorge]

@dseidita1919 Four-coil follow-up. The winding surfaces and current-potential solutions are unchanged; only the extracted contours are reduced from eight to four coils per half field period. This isolates the sparse-coil realization.

Actual independently evaluated filament maximum $\lvert B_n\rvert/B$:

Case	Normal offset	Entropy + strength	Full-sector entropy	Symmetric entropy
QA	0.1834	0.1747	0.1741	0.1742
QH	0.1334	0.1165	0.1141	0.1138
QI	0.7002	0.5673	0.5755	0.5758

The symmetric objective changes maximum $\lvert B_n\rvert/B$ by −5.0%, −14.7%, and −17.8% for QA, QH, and QI. Filament $f_B$ changes by +0.2%, −9.8%, and +3.7%; mean $\int\kappa^2dl/L$ changes by +5.0%, +5.0%, and −6.7%.

Resolved current potential; white contours are the coil cuts:

Resulting coils:

Winding surfaces and coils:

Symmetric entropy remains almost identical to full-sector entropy while using one SVD. It reduces the worst pointwise error in all three cases; integrated $f_B$ improves only for QH. Relative to eight coils on the same symmetric surfaces, maximum $\lvert B_n\rvert/B$ is 6.88×, 2.04×, and 2.37× higher. Four coils expose the limitation of using the surface-current spectrum alone for very sparse contour cuts.

Full numerical data

[rogeriojorge]

@dseidita1919 Last idea from these tests: keep the SVD solve, but make the surface objective know how many coils will be cut. “Finite-coil SVD” refers to the combination below; the SVD acts on the current-sheet operator.

1. Current-normalized SVD solve. For current-potential coefficients $x$,

$$B_n=A(S)x+b,\qquad f_K=(x-x_0)^TC(S)(x-x_0)+f_{K,0}.$$

$A$ maps current potential to normal field and $C$ measures integrated surface current. Factor $C=LL^T$ and set $y=L^T(x-x_0)$. Then $y^Ty$ is the additional current cost. The operator in these coordinates and its SVD are

$$\widehat A=A L^{-T}=U,\mathrm{diag}(\sigma_i)V^T.$$

Let $t=Ax_0+b$. The regularized current solve is

$$ \min_y\left[|\widehat A y+t|_W^2+\sigma_{\ast}^2|y|^2\right]. $$

Substituting the SVD and defining $\beta_i=u_i^Tt$ gives one scalar problem per mode,

$$ \min_{y_i}\left[(\sigma_i y_i+\beta_i)^2+\sigma_{\ast}^2y_i^2\right]. $$

Setting its derivative to zero gives

$$y_i=-\frac{\sigma_i}{\sigma_i^2+\sigma_{\ast}^2}\beta_i.$$

The factor $\sigma_i/(\sigma_i^2+\sigma_{\ast}^2)$ is the filter produced by the regularized minimization. Modes with $\sigma_i\gg\sigma_{\ast}$ cancel field efficiently. Modes with $\sigma_i\ll\sigma_{\ast}$ receive small coefficients because they require excessive current. The example fixes $\sigma_{\ast}$ from the initial 10 MA/m current budget.

2. Error from cutting contours. The continuous potential gives $\mathbf K=\mathbf n\times\nabla_s\Phi$. Coils are level curves of $\Phi$. Replacing the smooth potential by equally spaced levels is a staircase $Q(\Phi)$. Define the contour error by

$$\delta\Phi_\alpha=Q_{\Delta\Phi,\alpha}(\Phi)-\Phi.$$

Its normal-field error is $A\delta\Phi_\alpha$. For $N$ coils per half field period, $M=2N$ levels per field period and $\Delta\Phi=I_{pol}/(N_{fp}M)$. The staircase error is periodic in $\Phi/\Delta\Phi$. Its leading Fourier harmonic is

$$ \delta\Phi_c=\frac{\Delta\Phi}{\pi}\sin\left(\frac{2\pi\Phi}{\Delta\Phi}\right),\qquad \delta\Phi_s=-\frac{\Delta\Phi}{\pi}\cos\left(\frac{2\pi\Phi}{\Delta\Phi}\right). $$

A shift of the first contour changes the phase $\alpha$. Any phase is a combination of the sine and cosine above. Averaging over this arbitrary phase gives the surface objective

$$ J_N=|A\Phi+b|_W^2+\frac12\left(|A\delta\Phi_c|_W^2+|A\delta\Phi_s|_W^2\right). $$

The first term is the continuous-sheet error. The second estimates the added error from replacing that sheet by $N$ coils.

$I_{pol}$, $N_{fp}$ and $M$ are fixed, so $\Delta\Phi$ is fixed. The objective still changes because the surface changes $A(S)$, $C(S)$ and the solved $\Phi(S)$. Hence the contour paths and $A(S)\delta\Phi(S)$ change during optimization. The same number of equal-current coils is being placed in a better geometry. Distance, Jacobian and self-intersection terms remain geometry safeguards.

All surfaces below were re-solved with REGCOIL at 96x96, with the same mpol=ntor=6 basis and 10 MA/m limit. Coils were cut at 96x96 and evaluated independently at 48x48. W7-X includes the virtual-casing field.

case	four-coil change from normal offset	max $\lvert B_n\rvert/B$	curvature change
QA	$f_B$ -10.1%, max -9.0%	0.1669	+3.4%
QH	$f_B$ -22.1%, max -14.0%	0.1147	+0.8%
W7-X	$f_B$ -70.5%, max -36.0%	0.1803	-34.1%

On the same QA surface, ESSOS, REGCOIL and QUADCOIL give max $\lvert B_n\rvert/B=0.16684$, $0.16684$ and $0.16691$. This agreement separates the surface result from the choice of current solver.

Advantages: one physical field-error objective, explicit coil-count dependence, differentiable, and usually 1-2 s for the surface solve. Limits: leading-harmonic and equal-current approximation; omitted coil clearance, forces and ports; and 21.5 s spent in QH L-BFGS-B line searches. The absolute four-coil errors are still 0.11-0.18. A subsequent ESSOS filament run can refine these coils, with the added cost of weights and several initial conditions.

The minimal runnable example is examples/winding_surface_opt_3.py; the full data and figures are in examples/winding_surface_sparse_coils/.

[dseidita1919]

Fourier-coefficient induction matrix: a new SVD/entropy surrogate variant

This adds an alternative to "ESSOS entropy": instead of representing the winding surface's current potential Φ as dipole strengths at discrete grid points, it represents Φ as a truncated Fourier series (mpol=ntor=6, matching the REGCOIL validation basis used everywhere else in this study) and builds the induction matrix from those coefficients directly. Same outer optimization grid, same geometry regularizers/weights, same optimizer settings and initialization as "ESSOS entropy" — only the induction-matrix construction differs, so the comparison below isolates that one change.

Not folded into SURFACE_METHODS here (that would mean touching the SURFACE_METHODS-keyed color maps/bar-chart offsets in winding_surface_comparison_8_coils.py, a separate change). It's a self-contained script, examples/winding_surface_fourier_comparison.py, that reuses this study's own helper functions (make_surface, baseline_surface, geometry_metrics, dipole_kernel, build_operators, and worker() itself for the 48×48 step) so every metric is defined identically to what's already tracked. Full writeup, reproduction steps, and raw CSVs: examples/winding_surface_fourier_comparison/README.md.

Result (96×96, sheet fB — dipole entropy vs. the new Fourier method, both freshly rerun on the same machine)

Case	normal offset	ESSOS entropy (dipole)	ESSOS Fourier entropy (new)	improvement
Landreman-Paul QA	0.000375	0.000677	0.000502	26%
Landreman-Paul QH	25.20	17.48	11.90	32%
W7-X	4.105	2.499	2.347	6%

The Fourier method beats the dipole entropy method it's meant to replace on all three cases. It costs more wall-time per optimization (2–9× dipole entropy: 2.4s/8.0s/19.4s vs. 1.0s/1.3s/2.1s for QA/QH/W7-X), but stays well under 20s absolute. Resolution stability (48/56/64) is comparable to or better than the dipole method's own (see figures).

⚠️ Found while reproducing baselines: W7-X's tracked data has a staleness gap

"normal offset" and "ESSOS entropy"'s W7-X rows in both winding_surface_comparison_8_coils/data/*.csv and winding_surface_comparison/data/*.csv disagree substantially (2–4×) with a fresh rerun of the current code on this machine — very likely "REGCOIL adjoint" too, though that one wasn't independently rerun here (needs the legacy Fortran binary). This isn't a cross-machine numerical difference: the fresh rerun matches, to machine precision, the untracked winding_surface_comparison/output/ cache (dated Sep 4), which postdates the tracked data/ promotion (dated Sep 3) and was apparently never re-promoted — the exact gap WINDING_SURFACE_AGENT_NOTES.md already flagged. "ESSOS Pareto"'s W7-X row happens to be unaffected, which is presumably why this wasn't caught earlier.

Practical effect: any current comparison against W7-X's promoted "normal offset" or "ESSOS entropy" numbers is comparing against stale data (e.g. a draft of this analysis, before I tracked this down, reported the Fourier method beating dipole entropy on W7-X by 74% — the real, same-machine figure is 6%, used above). Promoting winding_surface_comparison_8_coils/data/*.csv and winding_surface_comparison/data/*.csv from the current output/ cache looks like a separate, independent fix worth doing before trusting W7-X numbers elsewhere in this PR.

Also new: a memory-watchdog pattern for this machine's 12 GB ceiling

Every REGCOIL-touching step in winding_surface_fourier_comparison.py runs as a subprocess under RSS-based polling (kills + logs if system-available memory or process RSS crosses a threshold), not resource.setrlimit(RLIMIT_AS, ...) — that was tried first and rejected after it killed a JAX/XLA worker using only ~1.9 GB RSS, since XLA's CPU backend reserves large virtual address space independent of actual physical usage. One data point (W7-X resolution-64 in the resolution-convergence check) is missing for both the dipole and Fourier methods, watchdog-aborted at the genuine edge of available memory on this machine (1.48 GiB available vs. a 1.5 GiB floor) — not a bug, just this machine's headroom right now.

[dseidita1919]

Following up on the Fourier-coefficient induction matrix comment above: the first working version cost 2–9× the dipole method's optimization wall-time. This is a performance-only follow-up — same method, same accuracy, same resolutions throughout, no new modeling assumptions.

1.8–5.3× faster, verified bit-for-bit on the physics

Profiling showed the per-iteration induction-matrix construction dominated (not one-time JIT compilation). Two resolution-independent, machine-precision-exact changes to fourier_induction_matrix:

1. nfp-fold periodicity reduction — the Fourier basis functions sin(m·θ − n·nfp·φ) are exactly periodic with period 2π/nfp in φ (the mode number is always an integer multiple of nfp), so only one field period of plasma/winding points is needed, and the kernel contribution from all nfp copies of each winding point can be summed before contracting with the potential basis — an O(nfp²) → O(nfp) reduction in the dominant cost.
2. Matmul-restructured dipole kernel — algebraically identical to the existing dipole_kernel, but replaces a (P,Q,3) pairwise-difference tensor with four (P,3)@(3,Q) BLAS matrix multiplications via the |a−b|² = |a|²−2a·b+|b|² identity.

Verified end to end: identical L-BFGS-B iteration and evaluation counts on all three cases (the optimizer takes the exact same path, just faster), sheet fB matching the pre-optimization version to 13–14 significant figures at both 48×48 and 96×96.

Case	Before	After	Speedup
Landreman-Paul QA	2.37s	1.32s	1.80×
Landreman-Paul QH	8.04s	1.87s	4.30×
W7-X	19.37s	3.63s	5.33×

Two other candidates were tried and rejected: swapping the SVD for eigvalsh on the smaller Gram matrix gave a negligible 2–7% speedup and produced NaN gradients from the near-degenerate tail of the singular-value spectrum (a real correctness failure, not a tradeoff); float32 precision corrupted that same tail by ~9 orders of magnitude and degraded the full gradient to only ~1e-5 relative accuracy. Both land in exactly the numerically fragile territory WINDING_SURFACE_AGENT_NOTES.md already flags for this codebase's SVD-gradient behavior, for less benefit than the two changes above.

Updated "vs other methods" comparison (same numbers as before, now cheaper to produce)

Case	normal offset	ESSOS entropy (dipole)	ESSOS Fourier entropy (new)	improvement
Landreman-Paul QA	0.000375	0.000677	0.000502	26%
Landreman-Paul QH	25.20	17.48	11.90	32%
W7-X	4.105	2.499	2.347	6%

These are the same physical numbers as the comment above — the point of this follow-up is the runtime column, not new results. (Full writeup, including why two sets of figures now live side by side in the repo — the comment above embeds the original ones by URL, so they're frozen rather than overwritten — is in examples/winding_surface_fourier_comparison/README.md.)

Also: resolution-convergence data is more complete now

Freeing some background memory on this machine let 4 of the 5 previously-watchdog-aborted resolution-64 points complete (Landreman-Paul QH's Fourier-method point, and W7-X's dipole-entropy and Fourier-entropy points). W7-X's "normal offset" row still doesn't fit after four separate attempts (6.5–7.5 GB peak RSS each time, against available headroom) and is the one remaining gap in the resolution-convergence figure — a genuine, repeatable memory ceiling on this machine for that specific computation (in unmodified, shared build_operators code, not something either of these two comments' changes introduces), not a bug.