Time-Adaptive PIROCK Method with Error Control for Multi-Fluid and Single-Fluid MHD Systems

TL;DR

This paper tackles numerical stiffness in multi-fluid, multi-species MHD models of the solar atmosphere by introducing PIROCK, a second-order partitioned Runge-Kutta method that unifies ROCK2 for diffusion, a third-order explicit Heun method for advection, and a two-stage SDIRK for stiff reactive terms, all under adaptive timestep control. By avoiding splitting errors and leveraging embedded error estimators, PIROCK delivers substantial gains in stability, accuracy, and computational efficiency over traditional explicit and Lie-splitting approaches, as demonstrated on Sod-Shock and 2.5D magnetic reconnection problems with multiple species. The results show PIROCK achieving orders-of-magnitude reductions in steps and wall time while maintaining comparable density and flux accuracy, and enabling chemical fractionation studies relevant to solar physics, such as the First-Ionization-Potential (FIP) effect. This work provides a versatile, efficient temporal integrator for MFMS and related stiff diffusion–reaction–advection problems, with broad applicability to single-fluid and multi-fluid MHD systems in solar and stellar contexts.

Abstract

The solar atmosphere is a complex environment with diverse species and varying ionization states, especially in the chromosphere, where significant ionization variations occur. This region transitions from highly collisional to weakly collisional states, leading to complex plasma state transitions influenced by magnetic strengths and collisional properties. These processes introduce numerical stiffness in multi-fluid models, imposing severe timestep restrictions on standard time integration methods. New numerical methods are essential to address these computational challenges, effectively managing the diverse timescales in multi-fluid and multi-physics models. The widely used time operator splitting technique offers a straightforward approach but requires careful timestep management to avoid stability issues and errors. Despite some studies on splitting errors, their impact on solar and stellar astrophysics is often overlooked. We focus on a Multi-Fluid Multi-Species (MFMS) model, which presents significant challenges for time integration. We propose a second-order Partitioned Implicit-Explicit Runge-Kutta (PIROCK) method that combines efficient explicit and implicit integration techniques with variable time-stepping and error control. Compared to a standard third-order explicit method and a first-order Lie splitting approach, the PIROCK method shows robust advantages in accuracy, stability, and computational efficiency. Our results reveal PIROCK's capability to solve multi-fluid problems with unprecedented efficiency. Preliminary results on chemical fractionation represent a significant step toward understanding the First-Ionization-Potential (FIP) effect in the solar atmosphere.

Figures (6)

• Figure 1: Comparison of the STS implementation in Eq. (\ref{['eq:defnaiv']}) (top) and the stable implementation based on Chebyshev recurrence in Eq. (\ref{['eq:defrec']}) (bottom) of the first-order Chebyshev method with optimal stability polynomial given in Eq. (\ref{['eq:defRs']}). The stability polynomial $R_s(z)$ defined in Eq. (\ref{['eq:defRs']}) for $s=10$ (see black curves) and the stability polynomials of internal stages (related to $k_j,j=1\ldots,s-1$, see color curves) are plotted as a function of $z$ which is assumed purely real here. The STS implementation exhibits poor internal stability, with internal stability polynomials oscillating with large amplitude resulting in instabilities with respect to roundoff errors, in contrast to the Chebyshev recurrence implementation defined in Eq. \ref{['eq:defrec']}, with favorable internal stability where the amplitude of the internal stage stability polynomials remains bounded by $1$.
• Figure 2: Stability domains as defined in Eq. \ref{['eq:defS']} of the schemes involved in PIROCK: ROCK2 (in gray. here $s=13$ internal stages) for diffusion terms, an explicit method of order $3$ for advection terms (dark gray), and the L-stable diagonally implicit Runge-Kutta method of order $2$ for stiff reaction terms (light gray). Right picture: zoom near the origin to illustrate that the advection method (dark gray) includes the portion of the imaginary axis. The stability domains are plotted for $z$ in the complex plane (horizontal axis: real part $\mathrm{Re}(z)$, vertical axis: imaginary part $\mathrm{Im}(z)$).
• Figure 3: Density distribution across the entire domain for PIROCK with different tolerances, the explicit approach, and reference solutions (top left). The top-middle panel depicts the distribution of the logarithm of the absolute difference, $\log_{10}(|\rho-\rho^{ref}|)$, between the different scenarios and the reference solution. The top-right panel shows the distribution of the logarithm of the absolute difference between analytical and reference solution $\log_{10}(|\rho-\rho^{ref}|)$. The bottom panels mirror the top ones but focus specifically on the shock region between x=0.82 and 0.88 (indicated by the two vertical blue lines in the top-left panel). The red dashed line represents the reference solution, and the black dashed line corresponds to the explicit. The remaining colors represent different PIROCK cases with varying tolerances.
• Figure 4: From left to right: 1) Distribution of the reference solution of the velocity $u_{z,H}$ in a restricted domain along the y direction ($y\in[0.2,0.8]$ Mm) in km.s$^{-1}$, 2D figures represent the logarithm of the difference between the reference and the numerical solution $\log_{10}(u^{ref}_{z,H}-u_{z,H})$ for PIROCK with tolerance $10^{-5}$, $10^{-2}$, splitting, and explicit. The last plot at the bottom represents the 1D distribution $\log_{10}(u^{ref}_{z,H}-u_{z,H})$ of all these cases at $y=0.5$ Mm and $z \in[1.5, 3]$ Mm at $t=50$ s.
• Figure 5: Same as Fig. \ref{['fig:comparison']} but for $\rho_{H}$.