BOUT++: a framework for parallel plasma fluid simulations

B. D. Dudson; M. V. Umansky; X. Q. Xu; P. B. Snyder; H. R. Wilson

Source-linked AI summary

BOUT++: a framework for parallel plasma fluid simulations

B. D. Dudson, M. V. Umansky, X. Q. Xu, P. B. Snyder, H. R. Wilson

arXiv:0810.5757v2Oct 31, 2008physics.plasm-phphysics.comp-ph

TL;DR

The paper addresses the need for flexible, scalable simulation of tokamak ELMs and related fluid models. It presents BOUT++, a modular implicit framework in curvilinear coordinates, and demonstrates accurate benchmarks, efficient scaling, and close agreement with ELITE for linear ELM simulations.

Problem
ELM filament formation and particle and energy loss remain incompletely understood, motivating improved simulation tools for tokamak plasmas.
Method
BOUT++ uses modular 3D fluid equations in curvilinear coordinates, fully implicit time integration, WENO differencing, and reduced ideal-MHD ELM simulations.
Results
BOUT++ reproduces ELM mode structures and growth rates close to ELITE and demonstrates stable, accurate tests with efficient scaling to thousands of processors.
Takeaways & Limitations
BOUT++ provides a flexible platform for testing different algorithms and physics models, including dissipationless initial-value simulations of ELMs.
Takeaways & Limitations
Accuracy near shocks needs further improvement, and vacuum-region modelling is important for correct simulation results.

Abstract

from arXiv · show

A new modular code called BOUT++ is presented, which simulates 3D fluid equations in curvilinear coordinates. Although aimed at simulating Edge Localised Modes (ELMs) in tokamak X-point geometry, the code is able to simulate a wide range of fluid models (magnetised and unmagnetised) involving an arbitrary number of scalar and vector fields, in a wide range of geometries. Time evolution is fully implicit, and 3rd-order WENO schemes are implemented. Benchmarks are presented for linear and non-linear problems (the Orszag-Tang vortex) showing good agreement. Performance of the code is tested by scaling with problem size and processor number, showing efficient scaling to thousands of processors. Linear initial-value simulations of ELMs using reduced ideal MHD are presented, and the results compared to the ELITE linear MHD eigenvalue code. The resulting mode-structures and growth-rate are found to be in good agreement (BOUT++ = 0.245, ELITE = 0.239). To our knowledge, this is the first time dissipationless, initial-value simulations of ELMs have been successfully demonstrated.

1. Introduction

BOUT++ is a modular parallel framework for 3D plasma-fluid simulations that separates user physics from reusable numerical and computational infrastructure. Its general coordinate and model support is paired with efficient testing and scaling for tokamak-focused applications.

BOUT++ provides an object-oriented C++ framework for parallel 3D fluid simulations with arbitrary equations in curvilinear coordinates.
Metric-tensor support covers slab, sheared-slab, cylindrical, and non-orthogonal flux-coordinate geometries.
The framework separates differential geometry, communication, and input/output from user-specified physics equations.
Compared with broader frameworks, BOUT++ reduces required parallel-computing knowledge and development time for working simulations.
BOUT++ targets tokamak ELM simulations, with tests covering accuracy, flexibility, ELM comparisons, and scaling to thousands of processors.

2. Physics overview

ELMs are abrupt tokamak-edge eruptions whose particle and energy deposition can damage material surfaces. Their transition from peeling-ballooning instability to filamentary loss remains incompletely understood, motivating complementary simulation approaches.

ELMs release particles and energy within sub-millisecond timescales and produce filamentary structures from the plasma edge.
The released particles and energy reach material surfaces, creating potential damage in future fusion devices.
ELMs occur near the stability boundary of the ideal-MHD peeling-ballooning mode, driven by edge pressure gradients and currents.
How peeling-ballooning modes develop into observed filaments and ultimately lose particles and energy is not fully understood.
Existing nonlinear codes reproduce experimental observations, while a complementary strategy begins with ideal-MHD theory and adds complexity progressively.

3. Structure of the program

BOUT++ organizes time integration, data management, parallel communication, geometry, differencing, and physics into replaceable modules. CVODE advances states through an implicit black-box interface, while the physics module supplies model-specific derivatives.

The program comprises modules for time integration, I/O, data handling, MPI communication, coordinate operators, differencing, and physics equations.
Modules can be modified independently through common interfaces, with the physics module designed for especially easy replacement.
CVODE passes the system state to BOUT++, which returns time derivatives for implicit time advancement.
CVODE uses a Newton-Krylov BDF implicit method for stiff problems, while physics-based preconditioning remains a planned extension.
The solver is wrapped in C++, allowing replacement by another package without changing the remaining code.
Runtime options configure simulation settings, while the equations remain defined in a compiled physics module.

3.3. Data handling

BOUT++ encapsulates field storage, arithmetic, memory optimization, and debugging in data objects that simplify fluid-code development. These abstractions retain flexibility while adding hidden performance optimizations and configurable diagnostics.

Data classes manage memory and domain loops, allowing scalar and vector fields to be manipulated with concise overloaded operations.
The data layer supports 3D and axisymmetric scalar and vector fields, with arithmetic, dot-product, cross-product, and mathematical-function overloads.
3.3.1. Optimisation: BOUT++ prioritizes flexibility despite the usual speed advantage of specialized programs because flexibility accelerates high-level algorithm development.
3.3.1. Optimisation: Memory recycling and copy-on-change reduce allocation, freeing, and unnecessary data copying behind the user-facing interface.
3.3.2. Debugging support: Debugging support includes runtime checks, operation tracking, segmentation-fault handling, and a fast message stack for locating faults.
3.3.2. Debugging support: 10-15% is the reported runtime increase from enabling runtime checking and the message stack, while checks can be disabled for longer tested runs.

3.4. Communication and topology

BOUT++ separates parallel communication from field data objects to preserve optimization flexibility, while using two-dimensional regular domain decomposition and internal topology handling.

Parallel communication is handled by dedicated communicator objects so communication can be overlapped with calculations when useful.Fields are grouped into communicator objects that perform the communications.
Domain decomposition currently uses a regular two-dimensional grid over x and y, with equal point counts per dimension on each processor.
Topology is represented through branch-cuts in the grid file, allowing non-rectangular tokamak X-point domains to be processed as rectangular local subdomains.Branch-cuts must coincide with processor boundaries.

3.5. Coordinate system

BOUT++ represents curvilinear geometries through metric tensors and supports field-aligned tokamak coordinates that reduce resolution requirements along elongated magnetic structures, while requiring twist-shift treatment of magnetic shear.

Metric tensors and Christoffel symbols are stored as global scalar fields and used by differential operators to implement coordinate-system dependence.Metric quantities can be loaded from grid files or set in physics modules.
Moving meshes could in principle be implemented by evolving metric tensors, but this has not yet been attempted.
The coordinate system is restricted to one symmetry direction, z, so its metric tensor components are currently two-dimensional fields.Fully three-dimensional metric tensors are planned but are not currently necessary for tokamak simulations.
Field-aligned meshes reduce the number of grid points needed along magnetic field-lines when structures satisfy k|| << k⊥.This alignment correspondingly reduces computational cost.
In field-aligned coordinates, ψ labels flux surfaces, while z labels field lines on each surface, with local shear determining the coordinate mapping.The covariant basis vector connects neighboring grid points.
After one poloidal circuit, the mesh shifts toroidally by 2πq radians and is remapped using the twist-shift condition.The safety factor q depends on ψ.
The twist-shift condition preserves periodicity of field values but not radial derivatives because magnetic shear creates a discontinuity in the radial basis vector.This discontinuity can affect differencing stencils at the matching location.
BOUT++ uses quasi-ballooning coordinates to avoid a special poloidal location, accepting non-commutativity and its associated coordinate-torsion correction for vector equations.The correction is a differencing artifact rather than a physical effect and is expected to be small in most tokamak simulations.

3.6. Differencing methods

BOUT++ separates partial-derivative calculations from metric-based differential operators and supports runtime-selectable central and upwind schemes, including third-order WENO for advection.

BOUT++ uses Method of Lines finite differencing, calculating derivatives separately in each dimension before constructing metric-dependent operators.
Central differencing options include second- and fourth-order methods, CWENO, and FFT derivatives in the z dimension.Methods can be selected separately by dimension at runtime.
Third-order WENO is available for advection and is designed to remain well-behaved at steep gradients while maintaining high accuracy.The scheme was used in the reported tests and ELM simulations.
Differential operators are divided into coordinate-independent operators and operators specialized for Clebsch coordinates with the magnetic field aligned to y.The advection operation v·∇f has a dedicated function because upwinding requires special numerical treatment.
The perpendicular-gradient operator is inverted efficiently by Fourier transforming in z and solving band-diagonal problems in x for each Fourier mode.The resulting systems are tridiagonal for second-order differencing.

3.7. Physics module

BOUT++ isolates physical equations in physics modules, while input-configurable boundary conditions and centralized problem-specific code support verification and rapid implementation of new models.

Physics modules are the only BOUT++ component that assigns physical meaning to variables and determine the equations solved.They implement initialization and per-time-step derivative calculations.
Table 1 compares analytic ideal-MHD expressions with their corresponding BOUT++ code.
Initialization defines state and derivative fields, loads initial profiles, registers communicated fields, and can add output variables.Communication objects and extra output variables are optional.
Boundary conditions for each evolving variable can be selected in the input settings file, including zero-value, zero-gradient, Laplacian, and symmetric conditions for scalar fields.More complicated coupled boundary conditions are a possible future extension.
Concentrating problem-specific code in one place lets users verify the equations being solved and implement new physical models quickly.

4. Test problems

BOUT++ is benchmarked on drift-wave, interchange, and Orszag–Tang vortex problems, reproducing analytic growth behavior and maintaining accurate, stable ideal-MHD solutions around shocks.

Test objectives: Three tests assess wave-phase fidelity, interchange growth rates relevant to ELMs, and numerical stability and accuracy for shocks.The resistive drift instability tests wave propagation; the interchange mode tests growth-rate recovery; the Orszag–Tang vortex tests ideal-MHD shock handling.
Resistive drift-wave instability: BOUT++ recovers the resistive drift-wave growth-rate peak as parallel conductivity varies, indicating accurate wave-phase simulation.The analytic results are compared with BOUT++ simulations across conductivity values.
Interchange mode: Interchange-mode growth rates are well reproduced for curvature radii R = 1 and 10 metres.For R = 1 m, the mode amplitude increases by 8 orders of magnitude without significant numerical noise or instability.
Orszag–Tang vortex: At t = 0.5, |∇· B| is 1.5 × 10−12, 4.6 × 10−12, and 1.8 × 10−11 on 128 × 128, 256 × 256, and 512 × 512 grids.These values are small relative to large-scale B/L values of approximately 10−1, demonstrating accurate preservation of ∇· B = 0.
Orszag–Tang vortex: The Orszag–Tang simulation runs to t = 10 with only a very slow increase in ∇· B.The solution remains stable beyond the t = 0.5 snapshot despite shock formation.
Orszag–Tang vortex: BOUT++ solutions closely match specialised MHD codes, but small oscillations occur at shock tops and improve with increased grid resolution.The oscillations do not grow during the simulation and arise because many terms use fourth-order central differencing.

5. ELM simulations

BOUT++ performs linear ELM simulations with reduced ideal MHD and benchmarks growth rates and mode structures against ELITE. The results are close, while differences remain near the plasma boundary because the codes model the vacuum region differently.

ELM model and setup: BOUT++ simulates ELMs by evolving vorticity, pressure, and parallel vector potential with high-β reduced ideal MHD.The simulations use field-aligned coordinates, fourth-order central differencing, and third-order WENO advection.
ELM model and setup: These are the first reported ideal-MHD time-dependent ELM simulations without intentionally introduced dissipation.Only numerical dissipation is present in the simulations.
Linear benchmarking: γ = 0.245ωA for BOUT++ versus γ = 0.239ωA for ELITE in the n = 20 linear benchmark.The initial single-mode perturbation first undergoes a transient before settling to an eigenmode with one growth rate.
Mode structure: Poloidal-mode amplitudes peak near resonant magnetic surfaces, consistent with analytic expectations.The resonant surfaces are locations where q(ψ) is rational.
Mode structure: BOUT++ mode envelopes agree well with ELITE, although transient extra modes and slightly sharper individual harmonics remain.The main discrepancy occurs near normalized ψ = 1, where ELITE uses a vacuum treatment and BOUT++ uses low-pressure ideal plasma.

6. Performance

BOUT++ performance was evaluated across problem sizes and processor counts, including large 3D ELM and 2D Orszag–Tang simulations. The tests show efficient scaling, while Laplacian inversion and time advancement remain important performance limits.

Performance goals: 10^7–10^8-variable simulations motivate evaluating BOUT++ runtime scaling with problem size and processor number.The linear ELM problem is used as the default scaling case.
Performance goals: 80–90% of wall-clock time is spent evaluating physics-module time derivatives, with the remainder spent in the implicit CVODE time advance.Because the scheme is implicit, timestep selection depends on accuracy and simulation behavior rather than the fastest grid waves.
Scaling with problem size: Runtime scales approximately linearly with x and z domain sizes, while y-direction scaling lies between O(n) and a higher-order dependence because fast parallel dynamics affect the timestep.The tests indicate O(n) behavior when grid size does not affect the timestep, with worse scaling when it does.
Scaling with number of processors: 25,165,824 evolving variables on a 256×256×128 ELM mesh achieve approximately 50% efficiency at 2048 processors relative to 32 processors.Efficiency declines beyond 256 processors as the Laplacian inversion becomes significant; the parallel Thomas-based inversion loses efficiency when processor subdivision exceeds available poloidal points.
Performance limitations: The current Laplacian inversion uses a serial-style Thomas algorithm, providing a test baseline but remaining inefficient for parallel computers.The paper identifies faster parallel inversion algorithms as future work.
Scaling with number of processors: More than 50% efficiency is maintained at 4096 processors in another case, while CVODE limits the 2D Orszag–Tang scaling and communication is not a significant bottleneck.The reported cache effect may improve efficiency as each processor handles fewer grid points.

7. Conclusions

BOUT++ provides a modular platform for flexible plasma-fluid modeling with implicit time integration and WENO shock handling. Tests demonstrate efficient scaling and close agreement between reduced-ideal-MHD ELM simulations and ELITE, while faster parallel Laplacian inversion remains needed.

Contributions: BOUT++ supports arbitrary numbers of scalar and vector fields, enabling rapid testing of new algorithms and physics models.Its modular design allows the solved fluid model to be changed easily.
Numerical methods: The code combines a fully implicit CVODE solver with WENO schemes for shocks, which were tested on linear and nonlinear problems.The implementation is stable in the presence of shocks, but accuracy near shocks requires further improvement.
Performance: Runtime scaling is O(n) when the timestep is unaffected by grid size and worse when it is affected, with hard scaling demonstrated to thousands of processors.The conclusions identify faster parallel Laplacian inversion as a key improvement area.
ELM simulations: Reduced-ideal-MHD linear ELM simulations reproduce mode structures and growth rates close to those from the ELITE linear MHD code.The paper describes these as the first dissipationless, apart from numerical, ELM simulations.