pith. sign in

arxiv: 2506.16676 · v7 · pith:GZ6C5BI6new · submitted 2025-06-20 · ⚛️ physics.plasm-ph · cs.PF

Fast solvers for Tokamak fluid models with PETSC

Mark F. Adams , Jin Chen , Benjamin Sturdevant This is my paper

Pith reviewed 2026-05-22 01:01 UTC · model grok-4.3

classification ⚛️ physics.plasm-ph cs.PF
keywords multigridtokamakMHDpreconditionerPETScstellaratordisruptionplasma simulation
0
0 comments X

The pith

A toroidal semi-coarsening multigrid solver offers robust performance for implicit MHD solves in tokamak fluid models.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper shows how to build an effective multigrid solver for the linear systems that arise when advancing the momentum equation implicitly in the M3D-C1 fusion simulation code. Standard block Jacobi preconditioning slows down or fails as more toroidal planes are used because it ignores long-range couplings along the torus. The authors exploit the regular grid structure to apply geometric multigrid coarsening only in the toroidal direction, keeping fine resolution where the magnetic field alignment matters. Tests on a runaway-electron disruption in the SPARC device and on a stellarator configuration demonstrate that the new approach matches or exceeds the speed of the old solver and succeeds in cases where block Jacobi diverges.

Core claim

By replacing the block Jacobi preconditioner with a geometric multigrid method that semi-coarsens in the toroidal direction, the implicit momentum solve in M3D-C1 becomes more robust and scalable for both axisymmetric tokamak disruptions and non-axisymmetric stellarator equilibria.

What carries the argument

semi-coarsening geometric multigrid applied only in the toroidal direction on the partially field-aligned grid

If this is right

  • Convergence rates remain stable as toroidal resolution increases.
  • The solver succeeds on stellarator models that defeat the production block Jacobi method.
  • Performance stays competitive on SPARC-scale disruption simulations.
  • The approach fits naturally into the existing PETSc-based linear solver infrastructure of M3D-C1.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Similar semi-coarsening ideas might apply to other toroidal devices or to the radial and poloidal directions with suitable adaptations.
  • Faster and more reliable solvers could shorten the time needed for parameter scans in fusion reactor design.
  • The method may generalize to other semi-implicit plasma fluid codes that share a toroidal grid structure.

Load-bearing premise

The grid remains regular and partially aligned with the magnetic field in the toroidal direction so that simple geometric coarsening captures the dominant error modes.

What would settle it

A simulation on a grid without regular toroidal spacing or without partial field alignment where the multigrid solver fails to converge faster or more reliably than block Jacobi.

Figures

Figures reproduced from arXiv: 2506.16676 by Benjamin Sturdevant, Jin Chen, Mark F. Adams.

Figure 1
Figure 1. Figure 1: Snapshots of fields during thermal quench and runaway electron generation (courtesy R. Datta) [PITH_FULL_IMAGE:figures/full_fig_p004_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: 4 plane case, iteration counts (left axis/bar) and solve times with setup (right axis/bar) [PITH_FULL_IMAGE:figures/full_fig_p005_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: 8 plane case, iteration counts (left axis/bar), solve times with setup (right axis/bar) vs time [PITH_FULL_IMAGE:figures/full_fig_p006_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: 16 planes case, iteration counts (left axis/bar), solve times with setup (right axis/bar) vs time [PITH_FULL_IMAGE:figures/full_fig_p006_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: 32 plane case, iteration counts (left axis/bar), solve times with setup (right axis/bar) vs time [PITH_FULL_IMAGE:figures/full_fig_p006_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: 64 plane case with GPUs, iteration counts (left axis/bar) and solve times with setup (right [PITH_FULL_IMAGE:figures/full_fig_p007_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: 64 plane, Iteration counts (left axis/bar) and solve times with setup (right axis/bar) vs time [PITH_FULL_IMAGE:figures/full_fig_p009_7.png] view at source ↗
read the original abstract

Multigrid (MG) is widely recognized as a highly effective solver for the model problem, the Laplacian, but textbook MG fails on most problems of interest. MG methods have been applied to complex, real-world applications with careful consideration of the physical model and discretization. This work develops the first step in applying MG methods to science and engineering relevant magnetohydrodynamics (MHD) tokamak models in the \textit{M3D-C1} https://m3dc1.pppl.gov fusion energy science code. The semi-implicit time integrator in \textit{M3D-C1} is composed of many linear solves. The implicit advance of the momentum equation is the most challenging and is the focus of this work. The current production solver in \textit{M3D-C1} is a block Jacobi (BJ) preconditioner within a Krylov solver, where blocks group degrees of freedom on planes of constant toroidal coordinate. BJ convergence degrades as the number of planes increases due to the spectral properties of the matrix preconditioned with BJ. The partially magnetic field-aligned, regular toroidal grid structure in \textit{M3D-C1} is amenable to semi-coarsening geometric MG in the toroidal direction. This paper develops such a solver and demonstrates competitive performance on a runaway electron model of a SPARC https://cfs.energy/technology/sparc disruption, and superior robustness on a stellarator model on which the BJ solver fails to converge.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

1 major / 2 minor

Summary. The manuscript develops a geometric multigrid solver with semi-coarsening in the toroidal direction, implemented via PETSc, for the implicit momentum equation in the semi-implicit time integrator of the M3D-C1 MHD code. It contrasts this approach with the existing block Jacobi (BJ) preconditioner, whose convergence degrades with increasing toroidal plane count, and reports competitive performance on a runaway-electron SPARC disruption test case together with superior robustness on a stellarator case where BJ fails to converge.

Significance. If the reported robustness holds under the stated grid assumptions, the work would provide a practical advance in scalable linear solvers for large-scale tokamak and stellarator fluid models, directly addressing a known bottleneck in disruption and equilibrium simulations. The reliance on standard PETSc infrastructure is a strength that supports reproducibility and adoption.

major comments (1)
  1. [MG solver construction and numerical results] The central claim that the partially magnetic field-aligned regular toroidal grid is amenable to effective semi-coarsening geometric MG rests on the unverified assumption that the chosen smoother damps high-toroidal-wavenumber error modes on the discrete momentum operator. No Fourier/symbol analysis, smoothing-factor computation, or spectral-radius verification (target <0.5) is supplied to confirm this damping; without it the superiority over BJ as plane count increases remains untested and load-bearing for the performance claims on both the SPARC and stellarator cases.
minor comments (2)
  1. [Abstract] The abstract states 'competitive performance' and 'superior robustness' without any quantitative metrics (iteration counts, residual histories, wall-clock times, or scaling data); adding at least one table or figure summarizing these quantities would strengthen the presentation.
  2. [Implementation details] Explicit specification of the PETSc MG components (smoother type, restriction/prolongation operators, coarsening ratio in the toroidal direction, and stopping tolerances) is needed for full reproducibility.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for the constructive review and the recommendation for major revision. The comment raises a valid point about the lack of theoretical smoothing analysis, which we address below by proposing specific additions to strengthen the presentation of the multigrid solver.

read point-by-point responses

  1. Referee: [MG solver construction and numerical results] The central claim that the partially magnetic field-aligned regular toroidal grid is amenable to effective semi-coarsening geometric MG rests on the unverified assumption that the chosen smoother damps high-toroidal-wavenumber error modes on the discrete momentum operator. No Fourier/symbol analysis, smoothing-factor computation, or spectral-radius verification (target <0.5) is supplied to confirm this damping; without it the superiority over BJ as plane count increases remains untested and load-bearing for the performance claims on both the SPARC and stellarator cases.

    Authors: We agree that a Fourier or symbol analysis of the smoother would provide additional theoretical support for the semi-coarsening strategy. However, the discrete momentum operator in the M3D-C1 MHD model incorporates variable coefficients arising from the partially field-aligned grid, density, and magnetic field, making a complete smoothing-factor computation or spectral-radius verification substantially more involved than for constant-coefficient model problems. The manuscript instead demonstrates the practical effectiveness of the approach through direct numerical comparisons: the geometric multigrid preconditioner maintains iteration counts and wall-clock performance as the number of toroidal planes is increased in the SPARC disruption case, while the block Jacobi preconditioner degrades; on the stellarator case the multigrid solver converges robustly where block Jacobi fails. To respond to this comment we will add a short subsection describing the chosen smoother (Chebyshev or SOR relaxation) together with supplementary numerical diagnostics that quantify residual reduction on representative high-toroidal-wavenumber modes extracted from the test problems. These additions will make the empirical evidence for damping more explicit without altering the core claims or results. revision: yes

Circularity Check

0 steps flagged

No circularity: standard MG application to existing grid and model

full rationale

The paper presents the construction of a semi-coarsening geometric multigrid solver in the toroidal direction for the M3D-C1 momentum equation, motivated by the known partially field-aligned regular toroidal grid. This is a direct engineering application of established PETSc MG components to a given discretization, followed by numerical demonstrations on SPARC and stellarator cases. No equations reduce a claimed result to its own inputs by construction, no parameters are fitted and then relabeled as predictions, and no load-bearing premise rests on a self-citation chain. The central performance claims are independently falsifiable through the reported timings and convergence behavior.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The paper relies on standard multigrid theory and the existing grid structure of M3D-C1 without introducing new free parameters or postulated entities.

axioms (1)
  • domain assumption Multigrid methods can be effectively applied to the discretized MHD equations with the given grid structure.
    Invoked when developing the semi-coarsening approach for the toroidal direction.

pith-pipeline@v0.9.0 · 5789 in / 1217 out tokens · 86034 ms · 2026-05-22T01:01:15.469628+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

16 extracted references · 16 canonical work pages

  1. [1]

    M. F. Adams. Multigrid Equation Solvers for Large Scale Nonlinear Finite Element Simulations . Ph.D. dissertation, University of California, Berkeley, 1998. Tech. report UCB//CSD-99-1033

    work page 1998

  2. [2]

    M. F. Adams and M. K. Knepley. A bespoke multigrid approach for magnetohydrodynamics models of magnetized plasmas in petsc, 2023

    work page 2023

  3. [3]

    M. F. Adams, R. Samtaney, and A. Brandt. Toward textbook multigrid efficiency for fully implicit resistive magnetohydrodynamics. Journal of Computational Physics , 229(18), 9 2010

    work page 2010

  4. [4]

    M. F. Adams, P. Wang, J. Merson, K. Huck, and M. G. Knepley. A performance portable, fully implicit landau collision operator with batched linear solvers. SIAM Journal on Scientific Computing , 47(2):B360–B381, 2025

    work page 2025

  5. [5]

    J. H. Adler, T. R. Benson, E. C. Cyr, P. E. Farrell, S. P. MacLachlan, and R. S. Tuminaro. Monolithic multigrid methods for magnetohydrodynamics. SIAM Journal on Scientific Computing , 43(5):S70– S91, 2021

    work page 2021

  6. [6]

    K. Bell. A refined triangular plate bending finite element. International Journal for Numerical Methods in Engineering, 1(1):101–122, 1969

    work page 1969

  7. [7]

    Chac´ on

    L. Chac´ on. An optimal, parallel, fully implicit newton–krylov solver for three-dimensional viscoresis- tive magnetohydrodynamicsa). Physics of Plasmas , 15(5):056103, 02 2008

    work page 2008

  8. [8]

    J. P. Freidberg. Ideal magnetohydrodynamics. Plenum Press,New York, NY, 12 1986

    work page 1986

  9. [9]

    J. P. H. Goedbloed and S. Poedts. Principles of Magnetohydrodynamics: With Applications to Labo- ratory and Astrophysical Plasmas. Cambridge University Press, 2004

    work page 2004

  10. [10]

    D. S. Harned and D. D. Schnack. Semi-implicit method for long time scale magnetohydrodynamic computations in three dimensions. Journal of Computational Physics , 65(1):57–70, 1986

    work page 1986

  11. [11]

    S. C. Jardin. A triangular finite element with first-derivative continuity applied to fusion MHD applications. J. Comp. Phys. , 200:133–152, 2004

    work page 2004

  12. [12]

    S. C. Jardin. Review of implicit methods for the magnetohydrodynamic description of magnetically confined plasmas. Journal of Computational Physics , 231(3):822–838, 2012. Special Issue: Computa- tional Plasma Physics

    work page 2012

  13. [13]

    Trottenberg, C

    U. Trottenberg, C. W. Oosterlee, and A. Sch¨ uller.Multigrid. Academic Press, London, 2000. A Appendix: The Resistive MHD Equations This appendix sketches the numerical methods in the evolution of the momentum equations with a focus on the algebraic form of the diagonal blocks of the Jacobian matrix in M3D-C1, which is of the most relevance to algebraic s...

    work page 2000

  14. [14]

    the velocity solve ... vn+1

    work page

  15. [15]

    the magnetic field solve ... Bn+1

    work page

  16. [16]

    the pressure solve ... pn+1 where vn+1 is solved as nR2(νi, ˙U)−n[νi, ˙χ] =··· nR2νi ˙ω=··· n R4 (νi, ˙χ) +n[νi, ˙U] =··· (28) The inner product is defined as (a,b )≡∇⊥a·∇⊥b =aRbR +aZbZ and the Poisson bracket as [a,b ]≡[∇a×∇b·∇ϕ] =R−1(aZbR−aRbZ) After parabolization: nR2(νi, ˙U)−n[νi, ˙χ] −θ2δt2{δW0 11(νi, ˙U) +δW1 11(νi, ˙U′) +δW2 11(νi, ˙U′′)+ δW1 12(ν...

    work page