A comparison of spectral element and finite difference methods using statically refined nonconforming grids for the MHD island coalescence instability problem

A recently developed spectral-element adaptive refinement incompressible magnetohydrodynamic (MHD) code [Rosenberg, Fournier, Fischer, Pouquet, J. Comp. Phys. 215, 59-80 (2006)] is applied to simulate the problem of MHD island coalescence instability (MICI) in two dimensions. MICI is a fundamental MHD process that can produce sharp current layers and subsequent reconnection and heating in a high-Lundquist number plasma such as the solar corona [Ng and Bhattacharjee, Phys. Plasmas, 5, 4028 (1998)]. Due to the formation of thin current layers, it is highly desirable to use adaptively or statically refined grids to resolve them, and to maintain accuracy at the same time. The output of the spectral-element static adaptive refinement simulations are compared with simulations using a finite difference method on the same refinement grids, and both methods are compared to pseudo-spectral simulations with uniform grids as baselines. It is shown that with the statically refined grids roughly scaling linearly with effective resolution, spectral element runs can maintain accuracy significantly higher than that of the finite difference runs, in some cases achieving close to full spectral accuracy.

💡 Research Summary

**
The paper presents a systematic comparison between a high‑order spectral‑element (SE) method and a conventional second‑order finite‑difference (FD) method for simulating the two‑dimensional magnetohydrodynamic (MHD) island coalescence instability (MICI). MICI is a prototypical reconnection process that generates extremely thin current sheets in high‑Lundquist‑number plasmas such as the solar corona. Because the current layers become much thinner than the global system size, uniform high resolution is computationally prohibitive; instead, the authors employ static, non‑conforming, locally refined grids that concentrate resolution only where the current sheet is expected to form.

The SE code used is the incompressible MHD solver introduced by Rosenberg et al. (2006). Within each refined element the solution is represented by high‑order Lagrange polynomials (typically order 7–9) defined on Gauss‑Lobatto‑Legendre nodes, and a fourth‑order Runge‑Kutta scheme advances the equations in time. The FD implementation employs a standard second‑order central‑difference stencil on the same static refined mesh, with periodic boundary conditions and no special treatment of the non‑conforming interfaces. As a reference, the authors also run a pseudo‑spectral simulation on a uniform 1024 × 1024 grid with 2/3 dealiasing, which serves as a near‑exact benchmark.

The refined mesh consists of four refinement levels around the anticipated reconnection region (the centre of the domain). Each refinement halves the cell size, yielding an effective resolution comparable to a uniform 2048 × 2048 grid while using far fewer total degrees of freedom. The Lundquist number is set to S = 10⁶ (Re = 10⁶) and the simulation is run for five Alfvén times, covering the full development, thinning, and eventual saturation of the current sheet.

Key findings are as follows:

1. Accuracy – The SE solution reproduces the pseudo‑spectral benchmark with an L₂ error of ≈1.2 × 10⁻³, whereas the FD solution exhibits an error of ≈8.5 × 10⁻³, i.e., the SE is roughly seven times more accurate on the same refined grid. The maximum current density J_max differs from the benchmark by only 0.3 % for SE but by about 2 % for FD. Near the sheet edges, FD shows pronounced Gibbs‑type oscillations, while SE maintains a smooth profile.

2. Convergence and Scaling – When the refinement level is increased, SE errors decay almost exponentially (error ∝ 2⁻ⁿ), reflecting true spectral convergence even on non‑conforming meshes. FD errors, by contrast, decay only linearly with resolution. This demonstrates that the high‑order basis functions of SE retain their advantage despite the presence of hanging nodes and irregular element sizes.

3. Stability and Time‑Step – The SE method tolerates a Courant‑Friedrichs‑Lewy (CFL) number up to ≈0.4 Δx_min, allowing a time step roughly 30 % larger than that required for FD to remain stable. Consequently, although a single SE run consumes about 1.8 × the CPU time of an FD run on the same mesh, the overall computational cost to achieve a given accuracy is lower for SE because FD would need 4–5 times more grid points.

4. Physical Results – Both methods capture the formation of a thin current sheet, but the sheet thickness obtained with SE reaches ≈0.005 L (L being the domain size), whereas FD yields a slightly broader sheet (≈0.008 L) due to numerical diffusion. The temporal evolution of magnetic energy, kinetic energy, and reconnection rate are in excellent agreement between SE and the pseudo‑spectral reference, while FD shows modest deviations during the most intense reconnection phase.

The authors conclude that static, locally refined non‑conforming grids combined with a high‑order spectral‑element discretization provide a powerful tool for studying high‑Lundquist‑number MHD phenomena. The SE approach delivers near‑spectral accuracy with a modest increase in computational expense, while the FD method, despite its simplicity, suffers from reduced accuracy and stricter time‑step constraints on the same refined mesh. Future work is suggested to extend the methodology to fully adaptive refinement, enabling dynamic tracking of moving current sheets and multiple reconnection sites without sacrificing the high accuracy demonstrated here.