{"id":"7a873b39-8ca8-4524-8bf6-38d0c4dbe057","arxiv_id":"2608.01750","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A differentiable, GPU-accelerated solver computes finite-toroidal-mode ideal MHD instabilities and their gradients, with a benchmark agreement against NIMSTELL on a stellarator case.","lead":"AGNI is a new computer code that accurately finds plasma instabilities in fusion devices and computes how they change as the machine's shape and settings change. It is the first differentiable solver of its kind for global (finite-n) instabilities, which makes it useful for optimizing stellarators and tokamaks.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Radial mapping Jacobian inconsistently applied in discretized δW: Appendix B radial-derivative blocks omit W_s^{-1} required by Eq. (26) and Eq. (51), so the assembled matrix may not be the stated variational discretization.","rationale":"The reader's weakest assumption is that the pseudospectral discretization yields a variational, pollution-free eigenvalue problem. My concern targets a specific, internal prerequisite: even the basic radial mapping factors in the discretized matrix are inconsistent between Section 3.2 and Appendix B. If the code implements the formulas as written, the variational claim is false, and the agreement with NIMSTELL would be coincidental or driven by unstated cancellation. This is load-bearing because the central contribution is precisely the reliability of the finite-n solver. I do not claim the code is wrong—only that the manuscript as written does not establish the correctness of the assembled matrix. The recommended test is a direct numerical comparison of the two assemblies, which would settle the issue. The verdict remains CONDITIONAL: the authors should correct or clarify the radial mapping in the discretization and demonstrate the corrected assembly on an analytic equilibrium before full acceptance. This is a partial agreement with the reader because I identify a more specific facet of the same weak assumption, rather than a wholly new concern.","tokens_in":20297,"tokens_out":10791,"duration_ms":122993,"concrete_test":"Assemble the potential-energy matrix A for the §5 equilibrium twice: (i) exactly as written in Appendix B with a single W in radial-derivative blocks, and (ii) with the radial blocks corrected to D_{ρ_s}=W_s^{-1}D_ρ and the required W_s^{-1} factor per Eq. (26) and Eq. (51). Compare the dominant eigenvalue from the two assemblies. If the eigenvalues differ by more than the 6% NIMSTELL discrepancy, the published discretization is not the one that produced the benchmark.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"Section 3.2 (Eqs. 25-26) introduces the radial mapping ρ_s = f(ρ) and states that every discretized term carries the quadrature weight W_s W, with W_s = diag(f'(ρ))⊗I. For radial derivatives, the chain rule gives D_{ρ_s} = W_s^{-1}D_ρ (Appendix A.2, Eq. 51). Substituting into Eq. (26) yields a factor W_s^{-1} in the radial-derivative blocks—since (D_{ρ_s})^T (W_s W) (D_{ρ_s}) = D_ρ^T W_s^{-1} W D_ρ—not a bare W. However, the explicit discretized expressions in Appendix B for the field-line bending terms (Eqs. 54-55) and the mixed term (Eq. 59) show only a single W in the radial-derivative terms (e.g., [(D_ρ ιψ'^2)^† ψ'/√g W gθθ/ψ' (D_ρ ιψ'^2)]), with no W_s^{-1}. Moreover, Eq. (59) uses D_{ρ_s} explicitly while Eqs. (54)-(55) use D_ρ, making the notation inconsistent. If the code follows the displayed formulas, the radial weighting of the Q² and compressibility terms is not the discretization of δW stated in Eq. (26), and the computed eigenvalue would not be the physical growth rate. The NIMSTELL benchmark is the only external check, but it cannot resolve a systematic error in the radial discretization because the code assembly is not independently verified; the paper also provides no analytic test case (e.g., cylindrical equilibrium with known eigenvalue) and no sensitivity scan over the mapping parameters for the benchmark case.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents AGNI, a GPU-accelerated, JAX-based finite-n ideal MHD stability solver and optimizer. The authors discretize the Bernstein energy principle using pseudospectral differentiation matrices in real-space coordinates obtained from DESC equilibria, form a Hermitian generalized eigenvalue problem for the plasma displacement, and solve for the most unstable mode with a shift-invert Lanczos method. They benchmark the growth rate and eigenfunction against the initial-value code NIMSTELL for a modified Landreman–Buller–Drevlak quasi-helically symmetric stellarator, reporting agreement for the dominant m = n = 4 interchange mode. They also verify reverse-mode automatic differentiation gradients against central finite differences, measure CPU/GPU timings for eigenvalue and gradient evaluation, discuss finite-precision limits for near-marginal eigenvalues, and propose two ways of enforcing incompressibility. The stated contribution is a differentiable finite-n MHD stability solver that can be used for gradient-based optimization of stellarators, tokamaks, and mirrors.","tokens_in":20816,"tokens_out":11397,"duration_ms":136008,"significance":"If the discretization is correct, this is a timely and useful contribution: it extends differentiable MHD stability optimization from the infinite-n ballooning limit to finite-n global modes, and it brings GPU acceleration and JAX-based reverse-mode gradients to a class of problems previously handled by non-differentiable finite-element or initial-value codes. The paper contains concrete numerical comparisons with an external code, finite-difference gradient checks, resolution scans, and hardware timing measurements, all of which are valuable for a methods paper. The finite-precision resolution criterion in Section 4.3 is a sensible practical addition. However, the significance is conditional: the central validation rests on a single equilibrium/mode comparison, and the radial-mapping discretization contains an inconsistency that must be resolved before the variational and benchmark claims can be fully accepted.","major_comments":[{"comment":"The radial mapping Jacobian is applied inconsistently. Equation (26) states that each discretized term carries the quadrature weight W_s W, and Eq. (51) gives D_{ρ_s} = W_s^{-1} D_ρ. For a term with two radial derivatives in the mapped coordinate, substituting (51) into (26) yields factors W_s^{-1} in the radial-derivative blocks, not a single W. However, the explicit discretized field-line bending terms in Eqs. (54)-(55) display D_ρ with only W, omitting the required W_s^{-1} factors, while Eq. (59) uses D_{ρ_s} explicitly and Eq. (61) uses C_{ρ_s} containing D_{ρ_s}. If the code assembles the displayed D_ρ forms, the radial weighting of Q^2 is not the discretization of δW stated in Eq. (26), and the computed eigenvalue is not demonstrably the physical growth rate. Please reconcile the notation, correct the formulas if needed, and add an analytic test case (e.g., a cylindrical equilibri","section":"Section 3.2, Eq. (26); Appendix A.2, Eq. (51); Appendix B, Eqs. (54)-(55), (59)"},{"comment":"The external validation is limited to one equilibrium and one mode, and the benchmark is not self-contained: the details of the NIMSTELL setup, reduction to ideal MHD, and convergence criteria are deferred to the under-review paper [31]. Table 3 shows the eigenvalue changing by more than an order of magnitude between the coarsest and finest resolutions (3.97×10^-6 at 8×24×4 versus 5.85×10^-5 at 40×48×16), yet no NIMSTELL value or uncertainty is given in the table, and the 6% growth-rate difference in Fig. 4 is presented without error bars. Because the radial-mapping issue above could represent a systematic discretization error, a single external match is insufficient to establish the variational, pollution-free claim. Please add an analytic test, a resolution study tied quantitatively to NIMSTELL, and a scan over the radial mapping parameters for the benchmark case.","section":"Section 5 and Table 3"},{"comment":"The gradient formula (46) is the standard expression for the derivative of an eigenvalue of a Hermitian problem at fixed matrix coefficients. The paper claims reverse-mode gradients with respect to boundary-shape and profile parameters 'without re-solving the equilibrium,' but the finite-difference comparison in Fig. 5 varies the boundary coefficient R_{b,10} and states that force balance is enforced for each equilibrium. It is not specified whether the DESC equilibrium solve is included inside the differentiated program and whether Eq. (46) is chained through the equilibrium solve, or whether the computed gradient is only at fixed equilibrium quantities. This distinction is central to the optimizer claim. Please state the chain-rule path explicitly and, where possible, verify the total derivative against finite differences.","section":"Section 5.2, Eq. (46)"}],"minor_comments":[{"comment":"The symbol ρ is used both for the physical normalized toroidal flux in [0,1] and for the computational Legendre variable in [-1,1] used with the mapping ρ_s=f(ρ). This makes equations such as (22), (25), and (26) hard to parse. Consider using a distinct symbol (e.g., u) for the computational coordinate.","section":"Section 3.1 and 3.2"},{"comment":"The expression for |δV| omits the Alfvén-speed and length normalization factors that appear in the physical growth rate γ = (√λ/a_N)(B_N/√(μ0 n_i M_i)). This is harmless for the normalized eigenfunction plots but should be stated if Eq. (45) is used standalone.","section":"Section 5, Eq. (45)"},{"comment":"The text says both codes solve for the 'most dominant n = 0 family' and then describes a 'm = n = 4' mode. The relation between the field-periodic n=0 family and the n=4 mode should be clarified for readers unfamiliar with stellarator mode classification.","section":"Section 5"},{"comment":"The claim that without the instability drive all eigenvalues lie below 4×10^-10 is stated verbally; a quantitative axis label or a table entry for the largest stable eigenvalue would make the finite-precision argument easier to verify.","section":"Section 4.3 and Figure 2"},{"comment":"The code is described as 'publicly available as a pull request' in the DESC repository. For a methods paper, a standalone repository with a DOI or versioned release would improve reproducibility and long-term access.","section":"Data availability"}],"recommendation":"major_revision","confidential_remarks":"The manuscript leans heavily on companion or under-review references for two central ingredients: the NIMSTELL benchmark details ([31]) and the adjoint gradient derivation ([35]). I would ask the editor to confirm that these references are available and that the benchmark description is self-contained. The radial-mapping inconsistency in Appendix B is the main technical risk; it is fixable, but it needs to be addressed with corrected formulas and an analytic test before the variational claim can be accepted."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Worth a look if you care about stellarator stability optimization. AGNI is the first solver I know that makes finite-n ideal MHD growth rates differentiable end-to-end in JAX, and the reverse-mode gradient test against finite differences is solid. The CPU/GPU timing comparison is useful, and the paper gives a careful description of the eigenvalue handling (shift-invert, incompressibility projection vs. increased Gamma). That part is real work and I'd credit it.\n\nThe soft spot is the radial mapping. Eq. (26) defines the quadrature with an extra weight W_s = diag(f'(rho)), but only for derivatives not in rho. For radial derivatives the chain rule (Eq. 51) inserts W_s^{-1} on one side. The discretized formulas in Appendix B for the field-line bending terms (Eqs. 54-55) show a bare W with D_rho, while the mixed term (Eq. 59) uses D_{rho_s} explicitly. Those two notations can't both be the discretization stated in Eq. (26). If the code follows the displayed formulas, the radial weighting of the Q^2 and compressibility terms is not the variational form, and the eigenvalue could be systematically off. The NIMSTELL comparison shows a 6% difference, which is in the right ballpark, but it's one equilibrium and one mode; you can't rule out a systematic radial-weighting error on the basis of that alone. There is no analytic test case, no sensitivity scan over the mapping parameters or epsilon, and the NIMSTELL details live in an unpublished companion paper. The code is only a PR at the moment.\n\nNone of this kills the paper: the concept is sound, the gradient checks are convincing, and the resolution scan shows a converged-looking eigenvalue. But the radial mapping has to be clarified or fixed, and the verification needs at least one analytic or cylindrical benchmark plus a mapping-parameter sweep. I'd send it to peer review with a request for those changes. The core idea is worth referee time; the current evidence just doesn't certify the numerics yet.","headline":"Novel finite-n MHD stability optimizer, but the radial mapping discretization is suspect and the validation is too thin to certify it.","tokens_in":21216,"tokens_out":9740,"would_cite":false,"duration_ms":102269,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":["52.35.Py","52.55.Hc"],"model":"deepseek-v4-flash","headline":"AGNI is a differentiable, GPU-accelerated finite-n ideal MHD stability solver that recovers a stellarator's dominant m=n=4 interchange mode and returns accurate growth-rate gradients.","keywords":["ideal MHD stability","finite toroidal mode number","stellarator","differentiable optimization","generalized eigenvalue problem","spectral differentiation","automatic differentiation","interchange mode"],"falsifier":"Select an equilibrium with a documented low-n kink or peeling instability and compute its growth rate with AGNI and with an independent finite-element stability code at matched resolution; if the AGNI eigenvalue converges to a different branch or disagrees beyond the grid-convergence tolerance, the variational discretization is not capturing the physical spectrum.","tokens_in":20195,"feed_emoji":"🧲","tokens_out":9646,"duration_ms":93250,"temperature":0.7,"pith_summary":"The paper presents AGNI, a stability solver built to make finite toroidal-mode-number ideal MHD instabilities directly optimizable. AGNI discretizes the ideal MHD energy principle in real space with spectral differentiation matrices, using geometry from an equilibrium calculation, and turns the stability problem into a Hermitian eigenvalue problem whose largest positive eigenvalue is the squared growth rate. The authors benchmark AGNI on a modified quasi-helically symmetric stellarator, recovering the dominant m=n=4 interchange mode with growth rate and eigenfunction in reasonable agreement with an established initial-value code. They also demonstrate reverse-mode automatic differentiation gradients of the growth rate with respect to equilibrium shape and profile parameters, agreeing with central finite differences. If correct, this closes a gap: previously only infinite-n ballooning stability could be optimized, while finite-n global modes had to be checked afterwards.","feed_headline":"Solver computes finite-n MHD stability and its gradients in one pass","feed_subtitle":"It recovers the dominant interchange mode and returns accurate growth-rate gradients for design optimization.","key_machinery":"The carrying object is the discretized ideal-MHD energy principle δWp = λ δK, expressed as a generalized Hermitian eigenvalue problem Aξ = λBξ. Derivatives are taken by spectral differentiation matrices—Fourier in poloidal and toroidal angles, Legendre-Lobatto in radius with a mapping that clusters nodes where the mode peaks—and all geometry and metric coefficients come from the equilibrium. Two structural facts make the method fast and differentiable: the kinetic energy matrix is block-diagonal after node-major reordering, so its Cholesky factorization is O(N) rather than O(N³), and the shift-invert Lanczos iteration needs only the most unstable eigenpair. Reverse-mode automatic differentia","core_discovery":"AGNI's central claim is that a variational, pollution-free eigenvalue problem for finite-n ideal MHD stability can be assembled directly in real space from spectral differentiation matrices and equilibrium metric coefficients, and that the whole pipeline can be made differentiable. Discretizing δWp = λ δK with Fourier bases in angle and a Legendre-Lobatto radial grid gives a Hermitian problem Aξ = λBξ; a shifted rescaled radial component regularizes the axis, and a shift-invert Lanczos iteration isolates the most unstable mode. On a modified quasi-helically symmetric stellarator, AGNI recovers the dominant m = n = 4 interchange mode with growth rate and eigenfunction in agreement with an ini","pith_inferences":["The paper leaves implicit that the same gradient machinery could be combined with quasisymmetry and transport objectives in one multi-objective stellarator design loop; this is a direct next application of the reverse-mode gradients.","Because the near-marginal noise floor sits near |λ| ≈ 10⁻¹⁰ in the stated normalization, an optimizer pushing a configuration to marginal stability would need higher precision or a target well above that floor.","A natural stress test not reported here is a low-n kink or peeling-like case where the spectrum is not interchange-dominated; that would separate the discretization's capacity from the specific benchmark equilibrium.","The scalability of AGNI ultimately hinges on the preconditioned matrix-free solver the authors list as future work; without it, the dense-LU shift-invert path limits resolution to what fits in GPU or CPU memory."],"forward_implications":["Finite-n ideal MHD stability becomes a gradient-optimizable objective, closing the gap with the existing infinite-n ballooning optimization capability.","Growth-rate gradients with respect to boundary shape and profile parameters are computed without re-solving the equilibrium, which is what makes high-dimensional design optimization tractable.","GPU evaluation is roughly an order of magnitude faster than CPU for both eigenvalues and gradients, limited mainly by memory at high resolution.","Incompressibility can be imposed inside the differentiable pipeline by raising the adiabatic index Γ, letting an optimizer target incompressible modes without dense projection matrices.","The solver extends in principle to tokamaks and mirrors once vacuum (external) mode coupling is added, as the authors state."],"supporting_citations":[{"why":"Supplies the ideal MHD energy principle that AGNI discretizes into the eigenvalue problem.","marker":"[5]"},{"why":"Provides the 3D stability variational formulation and notation that AGNI's δW expansion follows.","marker":"[14]"},{"why":"DESC equilibrium code supplies the flux-surface geometry and metric coefficients AGNI differentiates.","marker":"[21]"},{"why":"JAX provides automatic differentiation, GPU execution, and the linear algebra primitives.","marker":"[24]"},{"why":"NIMSTELL is the initial-value benchmark code that AGNI's growth rate and eigenfunction are compared against.","marker":"[25]"},{"why":"Defines the modified quasi-helically symmetric equilibrium used as the benchmark case.","marker":"[29]"},{"why":"Supplies the Lanczos implementation used for the shift-invert eigensolver.","marker":"[30]"},{"why":"Derivation of the reverse-mode gradient formula for eigenvalues used in the AD gradient calculation.","marker":"[35]"}],"fun_headline_variants":["AGNI: differentiable finite-n MHD stability solver for fusion devices","Finite-n MHD stability with gradients in one differentiable solver","AGNI computes finite-n MHD stability gradients without re-solving","AGNI: finite-n MHD stability and gradients, no re-solves"],"cache_read_input_tokens":2688,"weakest_assumption_plain":"The load-bearing premise is that the pseudospectral discretization of the MHD energy integral produces a variational eigenvalue problem whose most unstable eigenvalue is the true physical growth rate; this is supported by only one resolution scan and a single external benchmark.","fun_headline_variants_meta":{"raw":{"variants":["AGNI: differentiable finite-n MHD stability solver for fusion devices","Finite-n MHD stability with gradients in one differentiable solver","AGNI computes finite-n MHD stability gradients without re-solving","AGNI: finite-n MHD stability and gradients, no re-solves"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001639,"raw_usage":{"total_tokens":6379,"prompt_tokens":802,"completion_tokens":5577,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":546,"completion_tokens_details":{"reasoning_tokens":5502}},"tokens_in":546,"tokens_out":5577,"duration_ms":37476,"temperature":1.0,"reasoning_tokens":5502,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-04T21:40:50.095749+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Select an equilibrium with a documented low-n kink or peeling instability and compute its growth rate with AGNI and with an independent finite-element stability code at matched resolution; if the AGNI eigenvalue converges to a different branch or disagrees beyond the grid-convergence tolerance, the variational discretization is not capturing the physical spectrum.","supporting_citations":[{"cited_title":"An energy principle for hydromagnetic stability problems","cited_arxiv_id":null,"evidence_quote":"Supplies the ideal MHD energy principle that AGNI discretizes into the eigenvalue problem."},{"cited_title":"TERPSICHORE: A Three-Dimensional Ideal Magnetohydrodynamic Stability Program","cited_arxiv_id":null,"evidence_quote":"Provides the 3D stability variational formulation and notation that AGNI's δW expansion follows."},{"cited_title":"DESC: A stellarator equilibrium solver","cited_arxiv_id":null,"evidence_quote":"DESC equilibrium code supplies the flux-surface geometry and metric coefficients AGNI differentiates."},{"cited_title":"JAX: composable transformations of Python+NumPy programs","cited_arxiv_id":null,"evidence_quote":"JAX provides automatic differentiation, GPU execution, and the linear algebra primitives."},{"cited_title":"Semi-Implicit Stellarator Magnetohydrodynamics with Nodal Spectral Elements","cited_arxiv_id":"2606.28613","evidence_quote":"NIMSTELL is the initial-value benchmark code that AGNI's growth rate and eigenfunction are compared against."},{"cited_title":"Gradients of functions of large matrices","cited_arxiv_id":null,"evidence_quote":"Supplies the Lanczos implementation used for the shift-invert eigensolver."},{"cited_title":"An adjoint-based method for optimising MHD equilibria against the infinite-n, ideal ballooning mode","cited_arxiv_id":null,"evidence_quote":"Derivation of the reverse-mode gradient formula for eigenvalues used in the AD gradient calculation."}],"review_version":1}