Numerical methods for stellarator simulations in BOUT++

This paper presents recent enhancements to the BOUT++ simulation framework, including improvements to grid generation, physics models, and the library itself, to enable Flux Coordinate Independent (FCI) modeling of the Scrape-Off Layer in realistic stellarator geometries, specifically demonstrated with the Wendelstein 7-X device.

The Gist

Imagine trying to predict how heat and particles flow through a fusion reactor that looks less like a simple donut and more like a twisted, multi-lobed pretzel. This is the challenge of modeling the Stellarator, specifically the Wendelstein 7-X (W7-X), which is currently the most advanced fusion device on Earth.

The "Scrape-Off Layer" (SOL) is the messy, hot edge of this reactor where plasma touches the walls. In a simple donut-shaped reactor, the magnetic field lines are like neat, parallel train tracks. But in a stellarator, the tracks twist, turn, and sometimes even loop back on themselves in chaotic ways.

This paper describes a major software upgrade for BOUT++, a powerful computer program used to simulate these plasma flows. The authors, David Bold and Brendan Shanahan, explain how they fixed the software so it can handle the "pretzel" geometry of the W7-X without crashing or giving wrong answers.

Here is a breakdown of their work using everyday analogies:

1. The Problem: The "Train Track" vs. The "Maze"

The Old Way: Previously, simulation software tried to solve physics equations by laying a grid of lines perfectly along the magnetic field (like train tracks). This works great for simple shapes.
The Problem: In the W7-X, the magnetic field lines twist into islands and chaos. Trying to lay a straight grid on a twisting, chaotic path is like trying to draw a perfect grid on a crumpled piece of paper. The lines don't line up, and the math breaks down, especially near the walls where the plasma hits the target plates.

2. The Solution: The "Flux-Coordinate-Independent" (FCI) Method

Instead of forcing the grid to follow the magnetic field, the authors used a method called FCI.

• Analogy: Imagine you are trying to map the flow of water in a river with lots of whirlpools and islands. Instead of drawing your map lines along the water's path (which changes constantly), you draw a fixed, rigid grid over the whole river. You then calculate how the water moves across your fixed grid lines.
• The Benefit: This allows the computer to handle the chaotic, island-like magnetic fields of the stellarator without getting confused.

3. The Upgrades: Fixing the Engine

The authors didn't just use the FCI method; they had to rebuild the engine of the BOUT++ software to make it fast and accurate enough for real-world use.

A. The "Octagon" Diffusion (New Math)

• The Issue: When plasma diffuses (spreads out) near the sharp corners of the reactor walls, standard math tools fail because the grid has "corners" where the lines aren't smooth.
• The Fix: They created a new math tool based on Finite Volume.
• Analogy: Think of the grid not as a single point, but as a small, 8-sided room (an octagon) surrounding that point. Instead of guessing the flow at a single dot, they calculate how much "stuff" flows in and out of the walls of that 8-sided room. This ensures that nothing magically disappears or appears at the corners, keeping the physics honest.

B. The "Teamwork" Upgrade (Parallelization)

• The Issue: Simulating a whole stellarator requires massive computing power. The old software was like a team where everyone had to wait for one person to finish a task before the next could start. It was slow.
• The Fix: They rewrote the code so the computer could split the work up much better, using many processors at once (MPI).
• Analogy: Imagine a library. The old way was having one librarian check out books for the whole building. The new way is having 1,000 librarians, each handling a specific section of the library simultaneously. This made the simulation run much faster, allowing them to use thousands of computer cores at once.

C. The "Smart Boundary" (Handling the Walls)

• The Issue: When plasma hits the wall, it behaves differently than in the middle of the reactor. The old software struggled to handle these "edge cases," especially when the magnetic field lines hit the wall at weird angles or very short distances.
• The Fix: They implemented a "Leg-Value-Fill" method.
• Analogy: Imagine a game of "telephone" where you have to guess what a message is at the edge of the room. If the person is too far away, the message gets garbled. The new code uses a smart mathematical "guessing game" (Taylor expansion) to predict what happens right at the wall, even if the wall is very close. If the wall is too close, the code automatically switches to a simpler, safer method to avoid errors, acting like a safety net.

D. The "Grid Maker" (Zoidberg)

• The Issue: To run the simulation, you need a digital map (a grid) of the reactor. If the map has weird, squished cells, the simulation will be inaccurate.
• The Fix: They improved the grid generator (called Zoidberg).
• Analogy: Imagine trying to tile a floor that has a weird, curved shape. If you just use square tiles, you end up with tiny, jagged pieces at the edges that don't fit well. The new Zoidberg tool smooths out the curves and aligns the tiles so they fit perfectly, even around the complex "ports" and "islands" of the stellarator. It also uses data from the actual machine to make the map as realistic as possible.

4. The Result: A Stable Simulation

With these upgrades, the team successfully ran simulations of the W7-X that were stable and reached a "saturated phase" (meaning the turbulence settled into a realistic pattern).

Why does this matter?
Fusion energy is the "holy grail" of clean power. To build a working fusion power plant, we need to understand exactly how heat and particles behave at the edge of the reactor. If we get this wrong, the reactor walls could melt. This paper provides the "high-definition camera" and the "smart navigation system" needed to see and predict those edge conditions in the most complex fusion reactor ever built.

In short: They took a software tool that worked for simple shapes, gave it a brain transplant to handle chaotic 3D shapes, and tuned its engine to run fast enough to help us build the future of clean energy.

Technical Summary

1. Problem Statement

Modeling the Scrape-Off Layer (SOL) of stellarator fusion devices, specifically the Wendelstein 7-X (W7-X), presents significant numerical challenges due to their complex magnetic topology.

• Topological Complexity: Unlike tokamaks, stellarators like W7-X feature magnetic islands and chaotic field lines. The SOL is characterized by an "island divertor" where plasma intersects target plates at discrete toroidal locations.
• Limitations of Existing Methods: Traditional transport codes (e.g., EMC3-Eirene) are computationally cheap but rely on input diffusion coefficients ($D$ and $\chi$), lacking high-fidelity nonlinear physics. Higher-fidelity fluid turbulence models typically rely on field-aligned coordinates to efficiently calculate parallel gradients. However, in regions with magnetic islands and X-points (common in W7-X), field-aligned coordinates become singular or non-differentiable, rendering standard field-aligned approaches untenable.
• The Gap: While the Flux-Coordinate-Independent (FCI) method has been developed to handle complex topologies, implementing it efficiently in realistic stellarator geometries within the BOUT++ framework requires overcoming hurdles in grid generation, differential operator accuracy, parallelization, and boundary condition handling.

2. Methodology and Technical Improvements

The authors present a suite of numerical improvements to the BOUT++ framework and the associated grid generation tool (zoidberg) to enable realistic SOL simulations for W7-X.

A. Grid Generation and Geometry (zoidberg)

• Short Connection Lengths: Field lines entering and leaving the domain within a single slice create "short connection length" regions that are numerically unstable. The authors implemented a geometric algorithm using shapely (polygon intersection/union) to trace boundaries one slice in each direction, effectively removing these problematic regions from the computational domain.
• Poloidal Alignment: To improve mesh orthogonality, new alignment methods were introduced:
  • Minimizing Total Distance: Reduces the distance between inner and outer boundary points.
  • Orthogonal Projection with Smoothing: Projects outer points onto the inner boundary and smooths the poloidal position using a Gaussian filter. A user-defined parameter (mfi) controls the allowed variation in poloidal spacing, allowing for more orthogonal meshes.
• Elliptic Grid Generation: The governing equation for the grid was modified to include a relaxation coefficient ($\zeta$) for the $z$-direction derivatives. This balances smooth contours away from boundaries with uniform cell sizing, preventing cell clustering near the inner boundary.
• Boundary Smoothing: Implemented spline interpolation and spectral smoothing (Fourier series) to handle sharp corners and "pumping gaps" in the first wall geometry, preventing numerical oscillations.

B. Differential Operators

• Finite Volume Perpendicular Diffusion: The standard coordinate-transform approach for perpendicular diffusion ($\nabla \cdot a \nabla_\perp f$) fails at corners (X-points or physical targets) where the transformation is non-differentiable.
  • Solution: A finite volume approach was proposed. The cell is treated as an octagon. Gradients across the 8 surfaces are reconstructed using linear combinations of neighboring points.
  • Convergence: Method of Manufactured Solutions (MMS) tests show second-order convergence in smooth regions but only first-order convergence at corners. The total convergence is approximately 1.5 order, which is acceptable for the complex geometry.
• Pre-computation: The diffusion coefficients are pre-computed during grid generation, requiring only 5 parameters per cell at runtime, significantly reducing computational overhead.

C. Parallelization and Performance

• FCI Parallelization Challenge: In FCI, field lines intersect the next slice at arbitrary locations, requiring interpolation. The original BOUT++ implementation required entire parallel slices to reside on a single MPI task, limiting scalability.
• Linear Operator Reformulation: The FCI interpolation was re-implemented as a linear operator (sparse matrix-vector multiplication).
  • The interpolation stencil (originally complex) was reduced to a $4 \times 4$ stencil.
  • The operation is decomposed into matrix-matrix-vector multiplications. The matrix-matrix part is pre-computed; only the matrix-vector multiplication occurs at runtime.
  • Scalability: Using PETSc, the domain can now be split in the $x$-direction. This allows simulations on grids (e.g., $68 \times 36 \times 256$) to scale up to 1152 MPI processors.

D. Boundary Conditions

• Leg-Value-Fill (LVF): The LVF method uses Taylor expansions to extrapolate values to boundaries.
  • Analytical Solution: The linear systems for LVF are solved analytically using sympy code generation before compilation, ensuring efficiency and reducing bugs.
  • Robustness: A cutoff ($s/d_y < 0.01$) is applied to prevent divergence when the boundary is very close to the last grid point.
• Unified Interface: A new template-based API was developed to handle both Field-Aligned (FA) and FCI grids. This allows the same physics code (e.g., sheath boundary conditions) to run on both grid types, automatically falling back to lower-order schemes where high-order schemes are invalid (e.g., short connection lengths).

3. Key Results

• Convergence: MMS tests confirm the new differential operator achieves expected convergence rates (1.5 order globally due to corner singularities).
• Grid Quality: The new alignment methods significantly reduce cell shear and improve orthogonality compared to previous "original" implementations. The histogram of $\cos(\alpha)$ (angle between grid vectors) shows a strong peak near 1 (orthogonal) for the new methods.
• Scalability: The PETSc-based parallelization demonstrates that the FCI operator is independent of the number of processors regarding error norms. Wall-clock time per RHS evaluation scales well up to 1152 cores, though scaling is currently limited by the lack of $z$-direction domain decomposition in BOUT++.
• Simulation Stability: The authors successfully ran Hermes-2 turbulence simulations in W7-X geometry with reduced magnetic fields. The simulations reached a saturated turbulent state stably.
• Physics Validation: Preliminary results show electron temperature profiles in steady-state transport simulations that are comparable to EMC3-like setups, validating the implementation.

4. Significance and Outlook

• Enabling High-Fidelity Stellarator Modeling: This work removes the primary numerical barriers to applying high-fidelity fluid turbulence models (like BOUT++/Hermes) to realistic stellarator geometries. It moves beyond simplified models to capture nonlinear physics in the W7-X SOL.
• Unified Framework: The development of a unified interface for FA and FCI grids simplifies code maintenance and ensures that physics improvements (like sheath boundary conditions) benefit all simulation types.
• Future Work:
  • Hermes-3 Integration: Adapting the multi-component plasma code Hermes-3 to FCI to leverage these enhancements.
  • Drift Physics: Investigating the impact of drifts, which are suggested to be critical in W7-X.
  • Continuous Integration: The new features are being integrated into the BOUT++ and Hermes-3 CI pipelines with extensive unit tests and MMS tests to prevent regressions.

In conclusion, the paper establishes a robust, scalable, and accurate numerical foundation for simulating the complex SOL of stellarators, bridging the gap between simplified transport models and high-fidelity turbulence simulations for devices like Wendelstein 7-X.