examples: adjoint boundary sensitivities + SIMSOPT analytic gradient - #11
Draft
krystophny wants to merge 47 commits into
Draft
examples: adjoint boundary sensitivities + SIMSOPT analytic gradient#11krystophny wants to merge 47 commits into
krystophny wants to merge 47 commits into
Conversation
krystophny
force-pushed
the
simsopt-adjoint
branch
from
June 14, 2026 08:40
a86496f to
50ec5ab
Compare
This was referenced Jun 14, 2026
…proximafusion#580) * build: bump CMake abseil pin to 20260107.1 for Clang >= 21 The CMake FetchContent abseil pin (2024-08) fails to compile under Clang >= 21: absl::Nonnull SFINAE in absl/strings/ascii.cc and the numbers.cc nullability annotations are rejected by the newer frontend. Bump to the 20260107.1 LTS, which compiles cleanly under Clang 21.1.8 and GCC. Clang is the compiler required for the Enzyme autodiff build. The Bazel build keeps its own (BCR) abseil pin and is unaffected. * enzyme: opt-in Clang/Enzyme build option and AD smoke test Add VMECPP_ENABLE_ENZYME (OFF by default), which requires a Clang compiler and a ClangEnzyme plugin path and builds a self-contained autodiff smoke test. The test differentiates a scalar objective written over Eigen::Map'd caller buffers and checks reverse- and forward-mode Enzyme gradients against the closed form and central finite differences. enzyme.h documents the intrinsic ABI and the allocation constraint that shapes the differentiable kernels: Enzyme cannot track Eigen's aligned allocator, so differentiable paths use Eigen::Map over caller-owned buffers and avoid heap expression temporaries. With the option off the build is unchanged. * pybind: expose the unpreconditioned internal-basis gradient Add a precondition flag to VmecModel.evaluate (default true, unchanged behaviour). With precondition=false the forward model returns at the INVARIANT_RESIDUALS checkpoint, so get_forces() yields the raw, unpreconditioned force: the gradient of VMEC's augmented functional (MHD energy plus the spectral-condensation and lambda constraints) with respect to the decomposed internal-basis state. This is the consistent state/gradient pair an external optimizer needs to minimise in VMEC's own basis. The native solver's preconditioned search direction (precondition=true) is a different vector; the raw gradient is the equilibrium residual and vanishes at convergence. Tests: raw force is finite and differs in direction from the preconditioned force, and drops by >1e6 from the initial guess to the converged equilibrium. * examples: drive VMEC++ from external optimizers in the internal basis Treat the equilibrium as the root problem F(x) = 0, where F is the raw internal-basis force (gradient of VMEC's augmented functional) exposed by evaluate(precondition=False). Wire it to two solvers that reuse VMEC++'s forward model: native-style preconditioned descent and Jacobian-free Newton-Krylov (matrix-free Hessian information). Both reach the native solver's equilibrium. This is the external-differentiability path: VMEC++ as a differentiable equilibrium component an outside optimizer can drive. Quasi-Newton root-finders without a preconditioner diverge on this stiff system, which motivates exposing VMEC's preconditioner as an operator next. Tests assert both solvers reach force balance and recover the native energy and state. * pybind: expose VMEC preconditioner as an operator; preconditioned JFNK Add VmecModel.apply_preconditioner(v): applies VMEC's preconditioner M^-1 (m=1, radial, lambda steps) to a vector in the decomposed basis. M^-1 is VMEC's hand-built approximate inverse Hessian; this exposes it as a reusable linear operator for preconditioned Krylov / quasi-Newton and for the Hessian solve in adjoint sensitivities. It requires a prior evaluate(precondition=true), which assembles the radial preconditioner. Validated exactly: apply_preconditioner(raw force) equals the native preconditioned search direction; the operator is linear and, once assembled, state-invariant. Use it as the inner Krylov preconditioner in Newton-Krylov: on solovev (ns=11) this cuts force evaluations from 2242 to 505 (4.4x) versus unpreconditioned JFNK, converging to the same equilibrium. * pybind: Hessian-vector product inside VMEC++; internal Newton-Krylov Add VmecModel.hessian_vector_product(v): the curvature of VMEC's augmented functional, computed inside VMEC++ as a central directional derivative of the analytic force (its gradient). The force is exact; only the directional step is finite-differenced. Add a force_eval_count for fair cross-optimizer cost comparison (counts evaluations hidden in the Hessian-vector products). Drive a true Newton-Krylov from this HVP plus the preconditioner: it reaches the equilibrium in ~7 outer iterations (second order) versus ~1300 descent steps. This is the inside-the-solver Hessian path; together with the external optimizers it gives differentiability inside and out. Benchmark (solovev, ns=11, force evals counted in VMEC++): preconditioned descent 2606 evals 1302 iters Newton-Krylov (JFNK) 2243 evals Newton-Krylov (preconditioned) 507 evals Newton (VMEC++ HVP + M^-1) 9194 evals 7 iters The HVP-Newton's higher force-eval count (two evals per finite-difference HVP) is what the exact Enzyme Hessian will remove. * ideal_mhd_model: make computeMHDForces allocation-free The force kernel allocated 17 dynamic Eigen vectors per radial surface (the _o half-grid quantities and the avg/wavg surface averages). Move them to preallocated per-thread ThreadLocalStorage scratch and assign in place, so the radial loop allocates nothing. Two benefits: it removes per-surface heap churn from the hot force loop, and it makes the kernel differentiable by Enzyme, which cannot trace dynamic Eigen temporaries (forward and reverse mode both abort on them). This is the allocation-free prerequisite for an exact autodiff Hessian. Pure refactor, identical arithmetic. Verified bit-for-bit: vmec_standalone MHD energy unchanged on solovev (2.548352e+00) and cth_like_fixed_bdy (5.057191e-02). * examples: globalize HVP Newton with a backtracking line search The full Newton step overshoots on stiff 3D equilibria (cth_like stalled at the iteration cap with ||F|| ~ 5e-2). Add a backtracking line search on ||F|| so each step is damped to a decrease. With it the HVP-Newton converges on cth_like in 9 outer iterations (||F|| = 1.8e-10) and still converges solovev in 8. * dft_toroidal: make ForcesToFourier allocation-free The forces transform materialized two per-(surface,m,zeta) Eigen temporaries (tempR_seg, tempZ_seg) inside the inner loop. Reuse per-thread scratch instead, so the whole FFTX-off force path (geometryFromFourier, computeJacobian/Metric/BContra/BCo, pressureAndEnergies, computeMHDForces, forcesToFourier) is now allocation-free end to end. Same arithmetic as the previous .eval(); verified bit-for-bit: solovev 2.548352e+00, cth_like_fixed_bdy 5.057191e-02. * enzyme: exact autodiff of the VMEC Jacobian kernel (forward vs reverse) Demonstrate exact automatic differentiation of a real VMEC nonlinear kernel. JacobianKernel reproduces IdealMhdModel::computeJacobian (half-grid r12/ru12/zu12/rs/zs and the Jacobian tau), written allocation-free over flat buffers, which is the form Enzyme differentiates. For L = 0.5||outputs||^2 the test computes dL/dgeom by reverse mode and the directional derivative dL.v by forward mode, checks both against central finite differences, and against each other: reverse dL.v vs FD : 1.9e-9 forward dL.v vs FD : 1.9e-9 forward vs reverse : 2.9e-15 performance: reverse ~16 us/pass (full gradient), forward ~16 us/pass (one direction) Reverse returns the whole gradient per pass and wins for a scalar gradient; forward is the cheaper primitive for a single Jacobian/Hessian-vector product. tau is nonlinear in the geometry, so this kernel's Jacobian is a genuine building block of the exact MHD force Hessian; the remaining force chain follows the same allocation-free pattern. * ideal_mhd_model: share the Jacobian kernel between solver and autodiff Move the half-grid Jacobian arithmetic into jacobian_kernel.h (ComputeHalfGridJacobian), allocation-free over flat buffers. Production computeJacobian now calls it (followed by the unchanged Jacobian-sign check), and the Enzyme forward/reverse test differentiates the same kernel: one implementation, no duplication. Bit-exact: vmec_standalone MHD energy unchanged on solovev (2.548352e+00) and cth_like_fixed_bdy (5.057191e-02). Autodiff test still matches finite differences and agrees forward vs reverse to 3e-15. * ideal_mhd_model: share the metric kernel (gsqrt, guu, guv, gvv) Extract computeMetricElements into the shared, allocation-free kernel ComputeMetricElements (metric_kernel.h), over flat buffers, and call it from the solver. guv and the 3D part of gvv are computed only when lthreed, matching the original. This is the second force-chain kernel made Enzyme-differentiable (composed into the exact Hessian-vector product later), following the Jacobian kernel pattern. Bit-exact: vmec_standalone MHD energy unchanged on solovev (2.548352e+00, 2D) and cth_like_fixed_bdy (5.057191e-02, 3D path with guv/gvv). * ideal_mhd_model: share the contravariant-field kernel (bsupu, bsupv) Factor the bsupu/bsupv arithmetic out of computeBContra into the shared, allocation-free kernel ComputeBsupContra (bcontra_kernel.h). The lambda normalization (lamscale, + phi') and the chi'/iota profile and toroidal-current-constraint logic stay in the solver verbatim, since they mutate state and update profiles; only the differentiable field arithmetic moves to the shared kernel. Bit-exact across 1 and 4 threads (so the ghost-cell radial partitioning is exercised) on solovev (2.548352e+00, 2D) and cth_like_fixed_bdy (5.057191e-02, 3D). * ideal_mhd_model: share the covariant-field kernel (bsubu, bsubv) Extract the metric index-lowering (bsubu = guu B^u + guv B^v, bsubv = guv B^u + gvv B^v; guv absent in 2D) from computeBCo into the shared, allocation-free kernel ComputeBCo (bco_kernel.h). Bit-exact across 1 and 4 threads on solovev (2.548352e+00) and cth_like_fixed_bdy (5.057191e-02). * ideal_mhd_model: share the magnetic-pressure kernel Extract the field-dependent magnetic pressure |B|^2/2 = 0.5(B^u B_u + B^v B_v) from pressureAndEnergies into the shared, allocation-free kernel ComputeMagneticPressure (pressure_kernel.h). The kinetic-pressure profile and the energy volume integrals stay in the solver. Bit-exact across 1 and 4 threads on solovev (2.548352e+00) and cth_like_fixed_bdy (5.057191e-02). Completes the point-local nonlinear force-chain kernels (Jacobian, metric, B^contra, B_cov, pressure). * ideal_mhd_model: share the MHD force-density kernel Extract computeMHDForces' real-space force-density assembly (armn/azmn/ brmn/bzmn, and crmn/czmn in 3D, even+odd) into the shared, allocation-free kernel ComputeMHDForceDensity (mhdforce_kernel.h). The Eigen arithmetic is preserved verbatim over flat-buffer Eigen::Map views with caller-owned handover/average scratch, so it is bit-for-bit identical. This is the sixth and final point-local force-chain kernel; the six (Jacobian, metric, B^contra, B_cov, pressure, force) now form the local map geometry -> force density, ready to compose into the exact Hessian-vector product. (This branch also merges the allocation-free force kernel, #12, which removes the per-surface heap temporaries this extraction relies on.) Bit-exact across 1 and 4 threads on solovev (2.548352e+00) and cth_like_fixed_bdy (5.057191e-02). * enzyme: exact Hessian of the composed local force map Compose the six shared force-chain kernels (Jacobian, metric, B^contra, B_cov, magnetic pressure, MHD force density) into the single local map g: real-space geometry -> real-space force density, the nonlinear core of VMEC's force. The full MHD force is T^T . g . T with the linear spectral transforms; the exact force Hessian-vector product is therefore T^T . J_g . T . v, and this provides J_g by autodiff. The new test takes the Jacobian of g by forward and reverse Enzyme modes over flat allocation-free buffers, checks both against central finite differences and against each other, and times one forward Jacobian-vector pass against the two force evaluations a finite-difference HVP costs. * ideal_mhd_model: share the hybrid lambda-force kernel Extract hybridLambdaForce's full-grid lambda force (blmn, and clmn in 3D) into lambda_force_kernel.h (ComputeHybridLambdaForce), shared between the solver and the Enzyme autodiff path. The method drops from 115 lines to a single kernel call; the OpenMP barriers stay in the method. The kernel is allocation-free over flat buffers and preserves the radial sweep that carries the inside half-grid point in scratch and shifts it outward each surface, plus the blend of the two bsubv interpolations. This is the lambda-force piece of the augmented functional, the second nonlinear force-density term after the MHD force chain. * ideal_mhd_model: share the constraint-force kernels Extract the two local (non-transform) pieces of the spectral-condensation constraint force into constraint_force_kernel.h, shared between the solver and the Enzyme autodiff path: - ComputeEffectiveConstraintForce: gConEff = (rCon-rCon0) ru + (zCon-zCon0) zu (effectiveConstraintForce), skipping the axis surface. - AddConstraintForces: add the bandpass-filtered gCon back into the MHD R/Z forces and write frcon/fzcon (the constraint part of assembleTotalForces). The Fourier-space bandpass between them stays the shared free function deAliasConstraintForce; the free-boundary rBSq contribution stays in assembleTotalForces. Allocation-free over flat buffers. This completes the local force-density terms of the augmented functional (MHD + lambda + constraint), the nonlinear core of the exact Hessian. * enzyme: extend the composed-force Hessian test with the lambda force Add the hybrid lambda force (lambda_force_kernel.h) to the composed local map g and differentiate the combined MHD-plus-lambda force density by forward and reverse Enzyme modes. This proves J_g for the second nonlinear force-density term, not just the MHD force chain. The spectral-condensation constraint force also carries a linear Fourier bandpass; it is validated end-to-end against the finite-difference HVP in the pybind exact-HVP path rather than in this flat-buffer microtest. * apply pre-commit formatting (ruff, docformatter, clang-format) * apply pre-commit formatting (ruff, docformatter, clang-format) * apply pre-commit formatting (ruff, docformatter, clang-format) * apply pre-commit formatting (ruff, docformatter, clang-format) * apply pre-commit formatting (ruff, docformatter, clang-format) * apply pre-commit formatting (ruff, docformatter, clang-format) * apply pre-commit formatting (ruff, docformatter, clang-format) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * bazel: declare force-chain kernel headers in ideal_mhd_model (sandbox fix) * test: docformatter-format test_internal_gradient docstrings Satisfies the docformatter pre-commit hook (was failing CI). * test: docformatter-format external/internal optimizer test docstrings Satisfies the docformatter pre-commit hook (was failing CI). * ci: re-trigger (transient apt-403 on packages.microsoft.com) * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * ci: skip benchmark result upload on fork PRs (token is read-only) The 'Compare benchmark result' step uses github-action-benchmark with comment-on-alert and the GITHUB_TOKEN, which is read-only for pull requests from forks -> 'Resource not accessible by integration'. Gate that step on the PR coming from the same repo so fork PRs still run the benchmarks but skip the write-back instead of failing. * ci: build VMEC2000 from source so the compat test runs on numpy 2 The pinned vmec-0.0.6 cp310 wheel was f90wrapped against numpy 1.x. Under the numpy 2.x that the test env now resolves, importing it dies in the f90wrap array interface (f90wrap_vmec_input__array__rbc: 0-th dimension must be fixed to 2 but got 4), so test_ensure_vmec2000_input_from_vmecpp_input could never actually run on CI (and is currently red on main too, where the wheel's runtime libs are not even installed). Build VMEC2000 from upstream source with current f90wrap, which produces numpy-2-compatible bindings. The recipe mirrors SIMSOPT's own CI (hiddenSymmetries/VMEC2000, cmake/machines/ubuntu.json). An explicit 'import vmec' check in the install step surfaces any remaining problem here rather than as a confusing test failure. * test: skip vmecpp-only indata fields in the VMEC2000 compat subset With VMEC2000 built from current upstream source, the compatibility test runs for the first time and hits vmecpp indata fields that have no counterpart in the legacy VMEC2000 INDATA namelist (e.g. free_boundary_method), which raised AttributeError. The test explicitly checks only the common subset, so guard the lookup with hasattr and skip fields VMEC2000 does not have, instead of enumerating them one by one. * build: pin abseil to the 20260107.1 commit hash Pin the FetchContent abseil dependency to commit 255c84d (the exact commit behind the 20260107.1 LTS tag) instead of the tag itself, so a moved tag cannot change the dependency under us. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash, not the tag. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash for Clang >= 21. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash for Clang >= 21. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash for Clang >= 21. * ci: sync VMEC2000-from-source build, benchmark fork guard, abseil commit pin Bring this stack branch up to the corrected CI baseline (from #583/#564): - tests.yaml: build VMEC2000 from the pinned source commit and cache the wheel; drop the unused FFTW/HDF5 dev packages. - benchmarks.yaml: skip the result upload on fork PRs (read-only token). - test_simsopt_compat.py: skip vmecpp-only INDATA fields. - CMakeLists: pin abseil to the 20260107.1 commit hash for Clang >= 21. * ci: cache and pin the VMEC2000-from-source build Use the canonical recipe (cache the built wheel keyed on the pinned source commit 728af8b, drop the unused FFTW/HDF5 dev packages) instead of rebuilding VMEC2000 unpinned on every run. * ideal_mhd_model: mark Jacobian kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: hoist ForcesToFourier scratch out of the inner loop The allocation-free rewrite placed tempR_seg/tempZ_seg in a block-scope thread_local inside the (jF, m, zeta) inner loop, which emits a __tls_get_addr call and an init-guard branch every iteration. Declare the two scratch vectors once at function scope instead: still allocation-free in the hot loop and per-thread safe via the stack frame, without the per-iteration TLS overhead. Same arithmetic; cma and w7x wout are bit-for-bit unchanged. * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: hoist ForcesToFourier scratch out of the inner loop The allocation-free rewrite placed tempR_seg/tempZ_seg in a block-scope thread_local inside the (jF, m, zeta) inner loop, which emits a __tls_get_addr call and an init-guard branch every iteration. Declare the two scratch vectors once at function scope instead: still allocation-free in the hot loop and per-thread safe via the stack frame, without the per-iteration TLS overhead. Same arithmetic; cma and w7x wout are bit-for-bit unchanged. * ideal_mhd_model: hoist ForcesToFourier scratch out of the inner loop The allocation-free rewrite placed tempR_seg/tempZ_seg in a block-scope thread_local inside the (jF, m, zeta) inner loop, which emits a __tls_get_addr call and an init-guard branch every iteration. Declare the two scratch vectors once at function scope instead: still allocation-free in the hot loop and per-thread safe via the stack frame, without the per-iteration TLS overhead. Same arithmetic; cma and w7x wout are bit-for-bit unchanged. * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: hoist ForcesToFourier scratch out of the inner loop The allocation-free rewrite placed tempR_seg/tempZ_seg in a block-scope thread_local inside the (jF, m, zeta) inner loop, which emits a __tls_get_addr call and an init-guard branch every iteration. Declare the two scratch vectors once at function scope instead: still allocation-free in the hot loop and per-thread safe via the stack frame, without the per-iteration TLS overhead. Same arithmetic; cma and w7x wout are bit-for-bit unchanged. * ideal_mhd_model: hoist ForcesToFourier scratch out of the inner loop The allocation-free rewrite placed tempR_seg/tempZ_seg in a block-scope thread_local inside the (jF, m, zeta) inner loop, which emits a __tls_get_addr call and an init-guard branch every iteration. Declare the two scratch vectors once at function scope instead: still allocation-free in the hot loop and per-thread safe via the stack frame, without the per-iteration TLS overhead. Same arithmetic; cma and w7x wout are bit-for-bit unchanged. * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * ideal_mhd_model: mark Jacobian metric kernel buffers __restrict Raw double* kernel params over the same flat layout prevent the compiler from vectorizing the pointwise loop (assumed aliasing), so on w7x these kernels ran ~2x slower than the Eigen-expression code they replaced. The buffers never overlap; mark them __restrict to restore SIMD. Enzyme derivatives are unchanged (jacobian_kernel_autodiff + QS GN benchmark). * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. (cherry picked from commit 27d36d21e1dd8ea6f73127b95bdc81d529f81672) * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. (cherry picked from commit 27d36d21e1dd8ea6f73127b95bdc81d529f81672) * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. (cherry picked from commit 27d36d21e1dd8ea6f73127b95bdc81d529f81672) * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. (cherry picked from commit 27d36d21e1dd8ea6f73127b95bdc81d529f81672) * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. (cherry picked from commit 27d36d21e1dd8ea6f73127b95bdc81d529f81672) * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. (cherry picked from commit 27d36d21e1dd8ea6f73127b95bdc81d529f81672) * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance, so every other caller is unchanged) and have the two vmec_in_memory_mgrid_test comparisons pass 2e-7 for jcuru/jcurv only, keeping 1e-7 for all profiles and geometry. (cherry picked from commit 27d36d21e1dd8ea6f73127b95bdc81d529f81672) * output_quantities: compare jcuru/jcurv at a looser opt-in tolerance The free-boundary in-memory-vs-disk mgrid golden compares two independent solves. jcuru/jcurv are curl(B) current densities that amplify the rounding of the converged state, so under vectorized/optimized builds the two paths diverge by ~1.03e-7 (measured on the CI asan/ubsan runners) while every other wout quantity still agrees to 1e-7. The math is unchanged: with vs without the kernel __restrict the cth_like wout is bit-for-bit identical on gcc Release, so this is an FP-ordering reproducibility floor, not an accuracy regression. Add an opt-in current_density_tolerance to CompareWOut (default 0 = use the main tolerance…
…sion#619) moved vmecpp.run_continuation() logic directly into vmecpp.run() Co-authored-by: Philipp Jurašić <166746189+jurasic-pf@users.noreply.github.com>
…n#621) The Eigen3 migration (proximafusion#410) and hot-loop rework (proximafusion#454) replaced the fused scalar poloidal accumulation in the toroidal transforms with per-quantity Eigen .dot() calls, and the FFT path additionally allocates two Eigen vectors per innermost (m,k) iteration via .eval(). On the short theta axis (nThetaReduced ~9-16) and small ntor+1 this is a pessimization: benchmark-runs history shows ToroidalForcesToFourier regressed ~2x from the pre-proximafusion#410 baseline and never recovered, including at the flagship 12x12 FFT size. Two changes: - dft_toroidal.cc: restore the pre-proximafusion#410 fused-scalar-loop DFT code verbatim (single pass over theta reading each basis value once; the original "auto-vectorize was a pessimization" note is retained). Fixes the DFT-fallback resolutions and the fftx-disabled build. - fft_toroidal.cc: the FFT path only replaces the toroidal direction; its poloidal fill kept the .dot()+.eval() pattern. Fuse it into one allocation-free scalar pass. Measured (same-machine A/B, --config=opt, OMP=1): 12x12 FFT forces 1.57x faster (1.00e-3 -> 6.35e-4), 6x8 neutral. The c2r FourierToReal output scatter already uses plain segment += and is left unchanged. fft_toroidal_test and vmec_test pass. CI benchmarks will confirm the recovery and inform whether any FFT shapes should still fall back to DFT. Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Fix typos in docs
Run every solver once on explicit Solov'ev and CTH-like inputs. Keep the strict internal-state comparison on the reproducible 2D case, and verify force balance and energy in 3D without optional diagnostics or repeated evaluations. Co-authored-by: Philipp Jurašić <166746189+jurasic-pf@users.noreply.github.com>
krystophny
force-pushed
the
simsopt-adjoint
branch
from
July 9, 2026 20:30
839a8e6 to
eb39c01
Compare
Add the implicit-function adjoint that turns VMEC++ into a gradient-providing equilibrium component for SIMSOPT, the original goal. vmecpp_adjoint.py: for a converged fixed-boundary equilibrium F_I(x)=0, the boundary sensitivity of a scalar objective J follows from H_II lambda = dJ/dx_I, dJ/dx_B = dJ/dx_B - (dF_I/dx_B)^T lambda, with H the symmetric Hessian of the augmented functional. It is matrix-free via hessian_vector_product and apply_preconditioner (the SPD interior system is solved with preconditioned CG). One Hessian solve gives the whole boundary gradient, versus one equilibrium re-solve per boundary DOF for finite differences. simsopt_vmec_gradient.py: VmecEnergy wraps this as a SIMSOPT Optimizable whose dJ is the adjoint gradient, plus a gradient-cost benchmark. Verified: the adjoint gradient matches brute-force re-solve finite differences (rel 2.4e-4) and the SIMSOPT Optimizable's dJ matches finite differences of J (rel ~1e-6). On solovev (ns=11, 18 boundary DOFs) the adjoint boundary gradient costs 762 force evaluations versus 9112 for finite differences (12x), and the gap grows with the boundary DOF count.
Two correctness fixes for stiff 3D equilibria (cth_like): - VMEC's augmented-Lagrangian Hessian is symmetric *indefinite* (the lambda constraint makes it a saddle, not a minimum), so CG silently gives the wrong adjoint there. Use GMRES, which handles indefinite systems, for the H_II solve and the interior Newton solve. With a loose, restarted tolerance the adjoint solve stays cheap. - Add a backtracking line search to solve_interior so the interior re-solve (used by the SIMSOPT wrapper and the finite-difference reference) converges on 3D instead of overshooting. Verified with a directional-derivative check against a re-converged finite-difference reference: solovev 1.5e-4, cth_like 2.2e-2 relative; both previously agreed only in 2D. Boundary-gradient cost on solovev: 626 force evaluations (analytic adjoint) versus 10460 (finite differences).
krystophny
force-pushed
the
simsopt-adjoint
branch
from
July 9, 2026 20:38
eb39c01 to
2f258ed
Compare
Fix formatting in docstring for clarity.
* Make raw VMEC forces history independent Separate the legacy previous-residual m=1 projection used by native iteration from the exact constrained projection used by external force evaluations. Cover both directions across the residual threshold. * Use compact storage for constraint policy * Name the m=1 gauge policy explicitly
* Make raw VMEC forces history independent Separate the legacy previous-residual m=1 projection used by native iteration from the exact constrained projection used by external force evaluations. Cover both directions across the residual threshold. * Use compact storage for constraint policy * Report comparable external optimizer metrics Use the low-level VMEC model for the native benchmark row so every solver reports the same force norm and MHD functional. Count SciPy Newton-Krylov outer iterations through its callback and correct the preconditioner description. * Count native force evaluations directly * Satisfy const wrapper contract * Name the m=1 gauge policy explicitly * Keep optimizer metric checks with their example tests * pybind: count force evaluations at model boundary --------- Co-authored-by: Philipp Jurašić <166746189+jurasic-pf@users.noreply.github.com>
…roximafusion#631) docs: fix README 'As a Python package' example to reference vmec_output The prose after the snippet referred to an `output` object, but the snippet defines the run result as `vmec_output`. Following the prose verbatim (`output.threed1_volumetrics`) raises a NameError. Rename the reference to match the variable defined in the same snippet. Fixes proximafusion#630
…rong `myst` dep (proximafusion#632) * Initial plan * Initial plan Co-authored-by: jurasic-pf <166746189+jurasic-pf@users.noreply.github.com> * fix: prevent duplicate github-pages artifacts on job reruns and remove incorrect myst dependency - docs.yaml: use `github-pages-${{ github.run_attempt }}` as the artifact name in both upload-pages-artifact and deploy-pages steps. When the build-docs job is re-run after a transient failure, the previous attempt's artifact is still present in the same workflow run. The deploy-pages action uses `listArtifacts()` which returns all attempts' artifacts, causing it to find 2 artifacts named "github-pages" and throw "Multiple artifacts named 'github-pages' were unexpectedly found". Using a per-attempt unique name avoids this conflict. - pyproject.toml: remove the `myst` package from the docs optional dependencies. The `myst` package (v1.0.4) is a Google authentication library unrelated to Sphinx/MyST markdown. The correct package `myst_parser` is already listed. Having `myst` installed pulls in google-auth, google-auth-oauthlib, and downgrades urllib3 to <2.0, introducing unnecessary dependency conflicts. Co-authored-by: jurasic-pf <166746189+jurasic-pf@users.noreply.github.com> --------- Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com> Co-authored-by: jurasic-pf <166746189+jurasic-pf@users.noreply.github.com>
…nd include/eigen3/ src files (proximafusion#642)
…roximafusion#641) The existing MultiGridFreeBoundary test only checked that the run converged, so it didn't catch the proximafusion#330/proximafusion#640 vacuum_pressure_state_ carryover regression: the run still converged, just via a different (longer) path once the vacuum term switched on a stage too early.
…ximafusion#645) * Decouple NESTOR vacuum solve thread count from the radial solver The whole equilibrium solve runs inside one persistent OpenMP parallel region whose team size is capped at ns/2 (>=2 flux surfaces per thread). The free-boundary NESTOR vacuum solve ran inside that region and so inherited the radial thread count: at the first multigrid step (ns=5) only 2 threads, even with -t 16. But NESTOR is parallelized over the tangential boundary grid (nZnT ~ thousands of points) and could use the full thread budget. At coarse grids the vacuum solve was ~8x thread-starved (Fortran PARVMEC avoids this with a separate vacuum communicator, VNRANKS). Give the vacuum solve its own thread count, decoupled from the radial one: - vmec_adjust_vacuum_num_threads(max_threads, nZnT) = min(max_threads, nZnT); computed once, ns-independent. - Build the vacuum solvers (fb_vac_/tp_vac_) exactly once in Vmec::SetupVacuumSolvers(), sized to vac_num_threads_, independent of the per-radial-thread setup loop. - IdealMhdModel::update drives the solve from a nested parallel region: a single radial thread spawns a team of vac_num_threads_ threads, each running NESTOR on its tangential slice. omp_set_max_active_levels(2) enables nesting; a CHECK_EQ guards against an under-provisioned team (which would silently under-cover the tangential grid). - The checkpoint early-exit result is broadcast to the radial team via a shared HandoverStorage::vacuum_reached_checkpoint flag published by the omp single barrier. - Retarget the free-boundary component tests from fb_/tp_/num_threads_ to fb_vac_/tp_vac_/vac_num_threads_. Numerics change at ULP level due to the changed tangential reduction order; the free-boundary multigrid regression (cth_like_free_bdy) still passes without re-baselining. Verified directly: at ns=5 with -t 16 the radial team is 2 threads while the vacuum team is 16. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Update vmec.h --------- Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
) * Tolerate corrupted string variables in from_wout_file Some Fortran VMEC codes (e.g. PARVMEC) can leave cosmetic string fields such as curlabel filled with uninitialized memory instead of valid text, when they fail to read a corresponding mgrid attribute (observed with a single-coil-group mgrid file). from_wout_file previously ASCII-decoded every string/char netCDF variable unconditionally, so one such garbage field aborted loading an otherwise perfectly valid wout file. Catch UnicodeDecodeError per variable and fall back to an empty string with a logged warning instead, since these fields are never used in the physics. * Update __init__.py * Update test_init.py --------- Co-authored-by: Philipp Jurašić <jurasic-pf@users.noreply.github.com>
* Decouple NESTOR vacuum solve thread count from the radial solver The whole equilibrium solve runs inside one persistent OpenMP parallel region whose team size is capped at ns/2 (>=2 flux surfaces per thread). The free-boundary NESTOR vacuum solve ran inside that region and so inherited the radial thread count: at the first multigrid step (ns=5) only 2 threads, even with -t 16. But NESTOR is parallelized over the tangential boundary grid (nZnT ~ thousands of points) and could use the full thread budget. At coarse grids the vacuum solve was ~8x thread-starved (Fortran PARVMEC avoids this with a separate vacuum communicator, VNRANKS). Give the vacuum solve its own thread count, decoupled from the radial one: - vmec_adjust_vacuum_num_threads(max_threads, nZnT) = min(max_threads, nZnT); computed once, ns-independent. - Build the vacuum solvers (fb_vac_/tp_vac_) exactly once in Vmec::SetupVacuumSolvers(), sized to vac_num_threads_, independent of the per-radial-thread setup loop. - IdealMhdModel::update drives the solve from a nested parallel region: a single radial thread spawns a team of vac_num_threads_ threads, each running NESTOR on its tangential slice. omp_set_max_active_levels(2) enables nesting; a CHECK_EQ guards against an under-provisioned team (which would silently under-cover the tangential grid). - The checkpoint early-exit result is broadcast to the radial team via a shared HandoverStorage::vacuum_reached_checkpoint flag published by the omp single barrier. - Retarget the free-boundary component tests from fb_/tp_/num_threads_ to fb_vac_/tp_vac_/vac_num_threads_. Numerics change at ULP level due to the changed tangential reduction order; the free-boundary multigrid regression (cth_like_free_bdy) still passes without re-baselining. Verified directly: at ns=5 with -t 16 the radial team is 2 threads while the vacuum team is 16. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Update vmec.h * Fixed OMP nested behavior on older glibc versions --------- Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
…imafusion#650) Fix vmecpp CLI: restore original signature, keep convert as its own subcommand better-errors and convert-cli (proximafusion#648, proximafusion#649) auto-merged mid-review with a --help-visible "run" subcommand that changed the documented default invocation from `python -m vmecpp <input_file>` to requiring `python -m vmecpp run <input_file>`. That was never requested; the review only asked for `convert` to become its own subcommand. This dispatches on `sys.argv[0] == "convert"` instead of using argparse subparsers for the default path, so the top-level parser and its --help output are unchanged from before, while `python -m vmecpp convert <input_file>` gets its own dedicated argument parser.
…oximafusion#654) Fortran VMEC 8.52 stores the rollback backup (xstore <- xc) BEFORE interp() overwrites the state with the interpolated coarse-grid solution, so the first restart of a continuation stage silently discards the interpolated seed and the fine stage effectively re-solves from a cold start (measured: identical iteration counts with and without the multigrid ladder on w7x/cth). PARVMEC/VMEC2000 fixed this on 2017-01-24 ('SPH 012417: move this AFTER interpolation call' in initialize_radial.f). Adopt the same ordering as the VMEC++ default: the first rollback target of a continuation stage is now the interpolated coarse-grid solution. This is a deliberate deviation from VMEC 8.52 (documented in README.md and the solver AGENTS.md). CHANGELOG.md catalogs all iteration flow-control differences between VMEC 8.52 and PARVMEC/VMEC2000 (~2017) and where VMEC++ stands on each, so further adoptions can be made deliberately; notably the PARVMEC restart criteria (iteration_style='parvmec') must not become the default before their companion changes (pre-step store, rollback iteration counting) are implemented -- see CHANGELOG.md items 2-4.
agonal solver package (squashed) Standalone block-tridiagonal factorization/solve and FD block-Jacobian assembly from proximafusion#616, needed by the experimental 2D preconditioner in this stack. Will rebase away once proximafusion#616 lands on main. Convergence improvements: multigrid transition handling and experimental 2D preconditioner All env-gated; default behavior unchanged (verified: identical iteration counts and physics on cth/w7x; full pytest suite passes). - VMECPP_MULTIGRID_INTERP=cubic|cubic_rho: 4-point Lagrange radial interpolation of the scaled spectral coefficients at multigrid transitions (in s or rho=sqrt(s)), replacing 2-point linear; reduces the re-injected force imbalance ~50x. - VMECPP_DELT_START=<frac>: continuation stages enter the force iteration at frac * delt, avoiding the stage-entry instability at full delt that otherwise destroys the interpolated seed. - VMECPP_DELT_RECOVERY: delt grows back (2 percent per accepted step) after restarts, bounded by a per-stage stability ceiling learned from failures (downward-only ratchet), with a stagnation guard (rollback to best state when no new residual minimum for 100 iterations) that catches weakly unstable marginal delt within ~100 iterations. - VMECPP_PRESERVE_INTERP_SEED: the bad-Jacobian axis-recovery path restores the interpolated backup at reduced delt instead of discarding an interpolation-seeded stage after a clean first force evaluation. - VMECPP_PREC2D=coupled: poloidally coupled 2D preconditioner on top of the proximafusion#616 block-tridiagonal infrastructure (zeta-averaged angle-resolved coupling kernels, cached factorizations). Verified correct; documented negative performance result on strongly shaped 3D configurations -- see docs/convergence_study.md Finding 5. Together with the backup-ordering fix underneath, the bundle gives 1.6-4.3x end-to-end speedups on w7x (incl. mpol=16 and 3x-pressure variants to ns=499), cth and cma (to ns=801) at identical converged physics; see docs/convergence_study.md. Convergence study: mechanism analysis, validation matrix, and figures Findings 1-8: why the multigrid ladder saved nothing (state discard at the first restart + stage-entry delt instability + permanent delt loss), the fix bundle and its validation (1.6-4.3x, identical physics), the 2D preconditioner phase-1 negative result, Fortran lineage provenance, and the truncation-error analysis of the transition residual floor with Richardson extrapolation as the follow-up. Remove the experimental coupled 2D preconditioner again The zeta-averaged poloidally-coupled 2D preconditioner (Finding 5 of docs/convergence_study.md) was evaluated and came out negative: correct (identical converged physics, bit-compatible in diagonal-block mode) but slower on strongly 3D-shaped configurations, because the zeta-averaged kernels discard exactly the toroidal mode structure that makes w7x-class cases stiff. All measured speedups of this branch come from the multigrid transition bundle, none from the preconditioner. This removes the implementation from ideal_mhd_model / handover_storage and drops the imported block-tridiagonal package (squashed from PR proximafusion#616) that Fix InterpTest after the multigrid interpolation rewrite The rewritten InterpolateToNextMultigridStep no longer stores the interpolation scratch arrays (sj, js1, js2, s1, xint) on the Vmec object; drop the checks of those internals and keep the behavioral checks of the interpolated coefficients (xold/xnew) against the golden reference data, which still pass since the default linear interpolation is bit-compatible. Port the delt-recovery time-step control to the Python iteration loop Adds "delt_recovery" as a fourth iteration style in vmecpp._iteration: the VMEC 8.52 control plus the time-step recovery scheme prototyped in the C++ core on this branch (VMECPP_DELT_RECOVERY / VMECPP_DELT_START): - every revert ratchets a one-directional stability ceiling to 0.95x the step that just proved unstable, - every store grows the step 2% back towards min(user delt, ceiling), - a stagnation guard reverts gently when no new residual minimum appears for 100 iterations (weakly unstable step below the 100x leash), - solve_multigrid enters continuation stages at 0.5x delt (the ramp-in), via a new delt_start_fraction parameter of solve_equilibrium. On a run without reverts the step already sits at the user delt and the scheme is bit-identical to vmec_8_52 (asserted in test_styles_agree_when_no_restart). IterationState gains a delt_ceiling field so the ratchet can be traced. examples/iteration_dynamics.py now overlays all four styles (plus the learned ceiling in the time-step panel) on the cma ns=72 case. delt recovery: no regression on cold starts; compare styles over the multigrid ladder Gate the delt-recovery scheme (C++ VMECPP_DELT_RECOVERY and the Python "delt_recovery" style) to interpolation-seeded continuation stages: on a cold start there is no seed to protect and no ramp-in deficit to recover, and re-probing near the stability margin after genuine-marginality restarts costs iterations. Cold stages now run the unmodified 8.52 control, bit-identical trajectories (asserted in the tests). Measured, w7x ns=[51, 99] (Python loop): vmec_8_52 stages (1551, 9 restarts), (1817, 6) total 3368 delt_recovery stages (1551, 9), (1108, 1) total 2659 (-21%) Stage 1 is now identical to the baseline (was +15% before the gate); the continuation stage keeps the full win (-39%). C++ core confirms: final stage 1813 -> 1109 iterations with the env bundle. On easy ladders (cma [51, 99]) the cost is +4 iterations on the continuation stage, zero on the cold stage. examples/iteration_dynamics.py now compares the styles over a cma [51, 99] multigrid ladder via solve_multigrid (per-stage counts, stage-boundary markers, global iteration axis), which is where the styles actually diverge. solve_equilibrium is no longer exported as public API -- the supported entry points are solve_multigrid and iterate; tests import the inner solve from vmecpp._iteration directly. Expose the multigrid interpolation scheme in the API Promote MultigridInterpolationScheme (kLinear / kCubic / kCubicRho) from an anonymous-namespace detail of vmec.cc to a public enum in vmec.h, thread it through InitializeRadial / InterpolateToNextMultigridStep as an optional parameter (default: the VMECPP_MULTIGRID_INTERP environment variable, i.e. linear), and expose it via pybind (MultigridInterpolationScheme enum, VmecModel.refine_to(ns, interpolation=...)). vmecpp.solve_multigrid re-exposes it as a string literal: interpolation='linear' | 'cubic' | 'cubic_rho'. The 4-point Lagrange interpolants reduce the interpolation-added error of a continuation-stage seed down to the coarse grid's own discretization error (measured on cma [51, 99]: entry residual 2.7e-2 -> 2.2e-4, stage-2 iterations -11%); asserted in test_multigrid_interpolation_scheme_improves_seed. examples/iteration_dynamics.py gains an INTERPOLATION knob (default cubic). Port the interp-seed preservation guard to the delt_recovery style The C++ VMECPP_PRESERVE_INTERP_SEED guard: when a continuation stage seeded by multigrid interpolation goes unstable during its first iterations (at least one clean force evaluation happened, so the seed is usable), restore the backup -- the interpolated seed -- at a reduced time step instead of discarding it via the axis-recovery cold reset. In the Python loop most of this guard's C++ role is already structural: a mid-run bad Jacobian on a seeded stage goes through the 8.52 leash branch, which restores the backup, and since the backup ordering fix the backup IS the seed. The remaining gap was a non-finite/bad-Jacobian evaluation at iter2 > 1 with ijacob == 0, which previously fell into the axis-recovery cold reset; it now restores the seed (armed together with the rest of the delt_recovery scheme, i.e. only on seeded continuation stages). With this, all four env-gated C++ levers have Python counterparts: VMECPP_MULTIGRID_INTERP (solve_multigrid interpolation=...), VMECPP_DELT_START + VMECPP_DELT_RECOVERY + VMECPP_PRESERVE_INTERP_SEED (the delt_recovery style). Finding 9: stage-entry stability limit is equilibrium-dominated, weakly jump-dependent Numerical scan (4 equilibria x 4 jump ratios x 2 source resolutions x 6 entry fractions, identical-seed probes via the Python loop, script in examples/delt_start_study.py): cth and cma enter stably at the full user delt at any jump up to 4x; w7x and li383 are entry-unstable already at a 1.5x jump, and quadrupling the jump lowers their stability limit by only ~25%. Scheduling VMECPP_DELT_START on the jump ratio is therefore not the right lever; the fixed 0.5 entry is stable for every tested transition and delt recovery erases its cost. Follow-up: seed the next stage's entry delt from the previous stage's learned stability ceiling. Finding 10: tail anatomy; raise delt input bound to ]0,10] - examples/delt_tail_study.py: snapshot probes of the tail stability boundary (static across the tail; learned ceiling 20-40% stale; cth/cma bounded by the user delt itself) and force-residual decomposition (the preconditioned residual is 78-99% lambda through the 1e-6..1e-9 bulk; the highest resolved m at mid/edge radius dominates at 1e-12). - Raise the delt input bound from ]0,1] to ]0,10]: no formula requires delt <= 1 (the damping term is scale-invariant in delt) and the tail boundary of cth/cma sits at 1.5-2x their conventional input value; raising the user delt measured -14% (cth at 0.9) and -18% (cma at 0.7) total iterations on [25,99] ladders. - Ceiling forgetting (letting the learned stability ceiling drift back up in quiet stretches) is a confirmed negative result: structural probe/fail limit cycle, every case regresses (w7x stage 2: 1646 -> 4354 iterations); documented in Finding 10 and at the recovery site. Finding 11: lambda preconditioner hyperparameters and boost study - Consolidate the lambda preconditioner's inherited magic scalings into vmec_algorithm_constants.h as documented named constants (kLambdaPreconditionerDampingFactor, kLambdaPreconditionerZeroGuard, kLambdaHighMDampingReferenceM, kLambdaHighMDampingMaxPower), replacing the vague kEigenvalueAvoidanceFactor / kModeDampingLarge / -Small and the private dampingFactor; document that faclam is the diagonal second variation of the magnetic energy with respect to lambda_mn. - Add experimental env knobs VMECPP_LAMBDA_PRECOND_SCALE and
* Name the lambda preconditioner constants Replace the magic numbers in IdealMhdModel::updateLambdaPreconditioner with named constants in vmec_algorithm_constants.h, documented as the tuning hyperparameters of the lambda preconditioner: - kLambdaPreconditionerDampingFactor (was the undocumented private dampingFactor = 2.0 in ideal_mhd_model.h) - kLambdaPreconditionerZeroGuard (was a literal -1.0e-10; replaces the misnamed kEigenvalueAvoidanceFactor, which was never referenced) - kLambdaHighMDampingReferenceM and kLambdaHighMDampingMaxPower (were the literals 16.0 * 16.0 and 8.0; replace the unreferenced kModeDampingLarge / kModeDampingSmall) Also documents the faclam stiffness diagonal and the pFactor scaling, replacing the TODO(jons) placeholders. No behavior change. * Apply suggestion from @jurasic-pf
…n#658) Add an optional MultigridInterpolationScheme enum (linear, cubic, cubic_rho) for Vmec::InterpolateToNextMultigridStep and rewrite the transfer loop: - linear keeps the VMEC 8.52 2-point scheme (default, and forced when the coarse grid has fewer than 4 surfaces); cubic uses a 4-point Lagrange stencil in s; cubic_rho interpolates in rho = sqrt(s), the natural radial variable near the magnetic axis. The 4-point stencil is shifted inward at the radial-domain ends and reproduces old-grid values exactly on coinciding grid points. - The higher-order interpolants reduce the interpolated-seed force residual of a continuation stage down to the coarse solution's own truncation error
…ximafusion#663) The first iteration of a free-boundary continuation stage skips the whole vacuum block (the iter2 > 1 gate in IdealMhdModel::update, inherited from Fortran VMEC funct3d) while assembleTotalForces still applies the edge term with the freshly zeroed rBSq. The LCFS row therefore takes one step under the raw, unbalanced plasma pressure: the stage-entry residual is FSQR ~ 9 regardless of the interpolation scheme, the MHD energy drops 12 percent in a single step, the LCFS pressure mismatch DELBSQ blows up to ~400x its converged value, and the following ~50-100 iterations are spent on a NESTOR ring-down chasing the kicked boundary. The vacuum solution of the converged coarser stage is still exactly valid at stage entry: the radial interpolation changes neither the angular grid nor the LCFS geometry. A continuation stage is a hot restart in this respect, and the hot-restart path already sets vacuum_pressure_state_ = kInitialized for exactly this purpose. Doing the same in InitializeRadial on free-boundary continuation stages makes iteration 1 run a full NESTOR solve on the preserved boundary and apply a force-balanced edge term: stage-entry FSQR drops from 9.19 to 2.6e-5 and the transition ringing disappears (W_MHD and DELBSQ stay at their converged values through the transition). Measured (identical converged physics to all printed digits): - cth_like_free_bdy_multigrid [15,25]: second stage 343 -> 320 iterations (niter regression guard updated 344 -> 321) - solovev_free_bdy [16,32]: second stage 828 -> 630 iterations (-24%) Single-stage and fixed-boundary runs are unaffected (ns_old == 0 or lfreeb == false); the state transition only fires when the vacuum contribution was already fully active on the previous stage.
…ig recompute (proximafusion#666) Speed up output quantities computation
…rs to aid debugging (proximafusion#665)
…lently instead of raising) (proximafusion#672)
* Near axis test * Fix CI: add test data file referenced by near-axis iota test The test references near_axis_iota_nfp4.json but the previous commit added the orphaned iota_05_near_axis.json instead, causing a FileNotFoundError. Commit the file the test actually uses and drop the unreferenced one.
* Add QUASR free-boundary integration tests Use QUASR SIMSOPT configurations (surfaces + coils) as a free-boundary convergence test bed. For each configuration three physics regimes are exercised (vacuum, ~1% beta, ~2% beta with net toroidal current) at ns=[8,24,71], mpol=ntor=10. - boundary from the outermost QUASR flux surface, - external field via VMEC++'s own mgrid response table (built once per config and reused across profiles); the makegrid coils file is written directly so the field computation stays inside VMEC++, - phiedge from the enclosed vacuum toroidal flux (Biot-Savart). Non-convergence surfaces as a test failure on purpose: the suite doubles as a convergence diagnostic. Marked 'slow'; resolution and config IDs are overridable via environment variables for a cheaper tier. * Check in QUASR configs via Git LFS for network-free CI Commit the 13 QUASR SIMSOPT serial files under tests/data/quasr (tracked with Git LFS) so the free-boundary suite runs without network access. The loader prefers the checked-in copy and falls back to downloading from the QUASR database only if a requested ID is missing locally. The default ID set now covers all checked-in configurations. Exclude the data directory from the whitespace/EOL pre-commit fixers so the serial files stay byte-exact. * Add physics-correctness assertions to QUASR free-boundary tests For converged runs, validate the physics rather than convergence alone: - vacuum: the free-boundary LCFS reproduces the QUASR boundary, compared by enclosed volume (a parametrisation-invariant geometric measure; absolute shape agreement is mgrid-resolution limited), - finite beta: the magnetic axis shifts outboard (Shafranov shift) relative to the vacuum axis, - net current: the prescribed curtor appears in the equilibrium, - new cross-check test: the vacuum free-boundary equilibrium agrees with an independent fixed-boundary equilibrium on the same boundary (magnetic axis and volume). Solves are memoised in a module-scoped cache so each (config, profile) runs at most once even when reused across assertions (e.g. the Shafranov comparison against the vacuum axis). * Fix CI: deselect slow QUASR tests by default; exact phiedge; hardcode defaults - Deselect the long-running 'slow' suite by default (pytest addopts -m 'not slow'); it was running in CI and timing out. Run explicitly with 'pytest -m slow'. - Compute phiedge exactly as the coil vector-potential line integral (oint A.dl, Stokes) instead of a masked grid integral of B_phi, which was biased ~5% low and shrank the free-boundary plasma; the vacuum LCFS now reproduces the QUASR boundary to <1% in volume (verified NESTOR and only_coils agree, so the residual is not a NESTOR artifact). This also drops the matplotlib dependency. - Import SIMSOPT unconditionally (fail loudly if missing) rather than importorskip; the slow marker is the guard. - Remove the VMECPP_QUASR_* environment overrides now that good defaults are fixed in the module; widen the mgrid margin so current-carrying LCFS stay inside the grid. * Run QUASR free-boundary suite in CI; xfail on non-convergence The slow suite is deselected from the normal test workflow, so add a dedicated workflow that runs 'pytest -m slow tests/test_free_boundary_quasr.py' on every push to main and on manual dispatch (workflow_dispatch). Non-convergence is now marked xfail instead of failing, so the job stays green while remaining a convergence-diagnostic bed: a configuration that starts converging (e.g. after a solver improvement) surfaces as an xpass. Converged- but-unphysical results still fail hard. The Shafranov check is skipped (not failed) when only its vacuum reference is unavailable. * Add example: QUASR free-boundary cross-section plots + summary table examples/free_boundary_quasr_cross_sections.py solves the three profiles for each QUASR configuration (reusing the integration-test module's setup), saves a cross-section PNG per config overlaying each converged LCFS on the QUASR target boundary (magnetic axes marked), and prints/writes a summary table (status, volume, vol/target, magnetic-axis R, Shafranov shift, beta, ctor). Resolution, config IDs and OpenMP threads are CLI options so it can run cheaply.
proximafusion#538) * free_boundary: non-stellarator-symmetric (lasym) free-boundary support * output_quantities: correct lasym Nyquist normalization; tighten degenerate test * vmec_test: tighten lasym axisymmetric educational check Drop the accept-unconverged flag (the case converges to the 1e-11 force tolerance, identically across thread counts), tighten the scalar bound 1e-4 -> 1e-5, and assert the asymmetric physics directly: the b0 = rbtor0/Raxis on-axis-field identity (which guards the tmult = 0.5 Nyquist normalization), nonzero rmns/zmnc, and the off-midplane magnetic axis against the educational_VMEC reference. * vmec_indata: allocate asymmetric arrays in SetMpolNtor SetMpolNtor accesses the asymmetric coefficient arrays (raxis_s, zaxis_c, rbs, zbc) through std::optional::value() when lasym is set. An object assembled field-by-field via the Python bindings reaches it with those arrays still unset, so a lasym input threw "bad optional access", or, when mpol/ntor were unchanged, left them unallocated for the caller. Allocate them to their zero defaults at the current resolution when lasym is set and they are missing, before the early return and the resize, so both the resize and downstream indexing are well-defined regardless of how the object was built. * tests: add lasym exact-equivalence tests Drive a stellarator-symmetric configuration through transformations that need the asymmetric representation but reproduce the symmetric equilibrium's physics: lasym with zero asymmetric content (2D and 3D), a rigid z-shift, and a small toroidal rotation that mixes symmetric modes into asymmetric ones. Each compares volume, beta, and the iota profile against the symmetric run. * test_lasym: assert full rmnc/zmns geometry match in the reduces-to-symmetric tests * test_lasym: single-line the _assert_same_geometry docstring * outputs: fix lasym wout coefficients and save(); add rotation-law regressions Fixes proximafusion#675: the asymmetric rmns/zmnc/lmnc conversion now matches convert.f (m=0 sources and signs, m>0 cross-term signs, sm/sp lmnc half-mesh interpolation); the derived-coefficient transforms symmetrize the real-space fields before the cos-parity projection (symoutput); the bsubsmnc axis extrapolation and bsubsmnc_full are added; bsubvF now restores the pre-symforce blmn value, fixing bsubvmnc/bsubvmns and the derived currents against PARVMEC; the currvmns inner-point sqrt(s) denominator is fixed; save() handles the union-typed asymmetric annotations and writes atomically through a temp file so a failed save leaves no partial wout file. Regressions: off-grid boundary reconstruction of the rotated cth_like case, the rigid-rotation law across all wout Fourier pairs with PARVMEC-calibrated tolerances, a complete save/load round trip asserting every asymmetric array, and no-partial-file-on-failure. * ci: retrigger after runner loss in ubsan job
… intrinsic in README (proximafusion#677)
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
| // A physical inconsistency was detected deep in the MHD model (e.g. a | ||
| // degenerate flux-surface geometry or a free-boundary current mismatch) | ||
| // that the solver has no retry strategy for. | ||
| UNRECOVERABLE_ERROR = 5, |
There was a problem hiding this comment.
warning: invalid case style for enum constant 'UNRECOVERABLE_ERROR' [readability-identifier-naming]
Suggested change
| UNRECOVERABLE_ERROR = 5, | |
| kUnrecoverableError = 5, |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
Use VMEC++ as a differentiable component in an external optimizer: the
boundary-shape gradient
dJ/dx_Bby the implicit-function adjoint instead offinite-differencing over boundary DOFs.
examples/vmecpp_adjoint.py: partition the state into interior/boundary,converge the interior to force balance (
solve_interior), then one adjointsolve
H_II lambda = dJ/dx_I(GMRES preconditioned byM^-1) gives the fullboundary gradient. The interior Hessian is symmetric indefinite (the lambda
constraint is a saddle), so GMRES is used, not CG.
examples/simsopt_vmec_gradient.py: a SIMSOPTOptimizablewrapping it.This PR uses the finite-difference HVP. With the exact autodiff HVP (#23) the
same adjoint gets cheaper still (numbers below).
Verification (force evals counted in VMEC++, ns=11)
The adjoint computes the same gradient as finite-differencing over the boundary
but at a cost independent of the number of boundary DOFs: 7x fewer evals on
solovev (18 DOFs), 93x on cth_like (150 DOFs) with the FD HVP, and 25x / 263x
with the exact HVP (#23). Gradients agree with the FD reference to 3.9e-4 /
4.0e-2.
Stacked on #10 (HVP).