Large-Scale High-Lundquist Number Reduced MHD Simulations of the Solar Corona Using GPU Accelerated Machines

We have recently carried out a computational campaign to investigate a model of coronal heating in three-dimensions using reduced magnetohydrodynamics (RMHD). Our code is built on a conventional scheme using the pseudo-spectral method, and is parallelized using MPI. The current investigation requires very long time integrations using high Lundquist numbers, where the formation of very fine current layers challenge the resolutions achievable even on massively parallel machines. We present here results of a port to Nvidia CUDA (Compute Unified Device Architecture) for hardware acceleration using graphics processing units (GPUs). In addition to a brief discussion of our general strategy, we will report code performance on several machines which span a variety of hardware configurations and capabilities. These include a desktop workstation with commodity hardware, a dedicated research workstation equipped with four Nvidia C2050 GPUs, as well as several large-scale GPU accelerated distributed memory machines: Lincoln/NCSA, Dirac/NERSC, and Keeneland/NICS.

💡 Research Summary

The paper presents a comprehensive effort to accelerate a three‑dimensional reduced magnetohydrodynamics (RMHD) code—originally written in Fortran 90 and parallelized with MPI—by porting it to NVIDIA’s CUDA platform. The scientific motivation is to study Parker’s model of coronal heating at very high Lundquist numbers (S ≈ 10⁶–10⁸), which requires extremely fine current sheets and long integration times to obtain statistically steady states. In the original implementation, the pseudo‑spectral scheme’s two‑dimensional Fast Fourier Transforms (FFTs) dominate the runtime, consuming more than 80 % of the computational budget, and the memory footprint becomes prohibitive on conventional CPU clusters when the resolution is increased to 1024² × 128 grid points.

To overcome these limitations, the authors adopt a four‑step CUDA porting strategy: (1) replace all FFTs with the CUFFT library, which is an order of magnitude faster than the hand‑coded CPU version and dramatically reduces host‑device data transfers; (2) recycle intermediate arrays to stay within the limited GPU memory; (3) implement point‑wise arithmetic (Poisson brackets, linear terms, etc.) as simple CUDA kernels, ensuring that the non‑FFT portion of the code does not become a bottleneck according to Amdahl’s law; (4) preserve the existing MPI domain decomposition in the axial (z) direction, assigning one MPI rank to one GPU, thereby achieving a hybrid MPI‑GPU parallelism that scales both across nodes and within each node.

Performance is evaluated on a spectrum of hardware: a traditional CPU cluster (Carver at NERSC), a GPU‑enabled workstation (UAF with four C2050 cards), and three large‑scale GPU‑accelerated clusters—Lincoln/NCSA (44 nodes × C2050), Dirac/NERSC (44 nodes × C2050), and Keeneland/NICS (120 nodes × 3 C2070 each). Relative to Carver, the CUDA version runs roughly 10× faster on Lincoln, 10× on Dirac, and up to 30× faster on Keeneland. Strong‑scaling tests show near‑linear speedup as the number of GPUs increases, confirming that the hybrid MPI‑GPU model efficiently utilizes both inter‑node and intra‑node parallelism.

Beyond performance, the scientific results are noteworthy. By extending the resistivity η down to ≈10⁻⁴ (far below the η ≈ 10⁻³ limit of earlier 2‑D studies), the authors observe a clear deviation from the η⁻¹⁄³ scaling of the average Ohmic dissipation ⟨W_η⟩ and the perpendicular magnetic field strength ⟨B_⊥⟩ reported by Longcope & Sudan (1994). The 3‑D simulations reveal that Ohmic dissipation dominates viscous dissipation (⟨W_ν⟩ ≪ ⟨W_η⟩) and that the time‑averaged Poynting flux matches the total dissipation within statistical error, indicating a well‑balanced energy budget. These findings suggest that three‑dimensional effects and long‑time statistical averaging modify the scaling laws derived from 2‑D models, providing new insight into the role of fine current sheets and magnetic reconnection in coronal heating.

In summary, the work demonstrates that GPU acceleration via CUDA can transform a computationally intensive RMHD code into a tool capable of tackling the extreme parameter regimes required for realistic coronal heating studies. The porting methodology—leveraging CUFFT, memory reuse, lightweight kernels, and a preserved MPI decomposition—offers a template for other plasma‑physics and astrophysics codes. Future directions include exploiting newer GPU architectures (e.g., NVIDIA Ampere with NVLink), extending the code to include additional physics (e.g., compressibility, Hall effects), and performing direct comparisons with solar observations to further validate Parker’s heating scenario.