Structure preserving hybrid Finite Volume Finite Element method for compressible MHD

TL;DR

This work tackles robust simulation of compressible viscous and resistive MHD across all Mach numbers by introducing a structure-preserving semi-implicit method that merges a conservative finite-volume treatment for convective terms with a compatible finite-element exterior calculus discretization for magneto-acoustic terms. The spline-like FEEC framework enforces discrete energy stability, magnetic-helicity conservation, and a divergence-free magnetic field, while the time integration splits into an ADI-like three-operator scheme solving a sequence of linear, symmetric positive-definite systems via matrix-free conjugate gradients. The method is validated on a broad suite of 1D, 2D, and 3D tests, including Brio–Wu, Orszag–Tang vortex, rotor problems, and long-time stationary vortices, demonstrating robustness, accuracy, and low dissipation even at large CFL numbers. The results indicate a flexible and scalable approach that can be extended to higher-order accuracy and adaptivity, with potential applications in astrophysical and fusion-plasma contexts where multi-scale MHD dynamics are essential.

Abstract

In this manuscript we present a novel and efficient numerical method for the compressible viscous and resistive MHD equations for all Mach number regimes. The time-integration strategy is a semi-implicit splitting, combined with a hybrid finite-volume and finite-element (FE) discretization in space. The non-linear convection is solved by a robust explicit FV scheme, while the magneto-acoustic terms are treated implicitly in time. The resulting CFL stability condition depends only on the fluid velocity, and not on the Alfvénic and acoustic modes. The magneto-acoustic terms are discretized by compatible FE based on a continuous and a discrete de Rham complexes designed using Finite Element Exterior Calculus (FEEC). Thanks to the use of FEEC, energy stability, magnetic-helicity conservation and the divergence-free conditions can be preserved also at the discrete level. A very efficient splitting approach is used to separate the acoustic and the Alfvénic modes in such a fashion that the original symmetries of the PDE governing equations are preserved. In this way, the algorithm relies on the solution of linear, symmetric and positive-definite algebraic systems, that are very efficiently handled by the simple matrix-free conjugate-gradient method. The resulting algorithm showed to be robust and accurate in low and high Mach regimes even at large Courant numbers. Non-trivial tests are solved in one-, two- and three- space dimensions to confirm the robustness, accuracy, and the low-dissipative and conserving properties of the final algorithm. While the formulation of the method is very general, numerical results for a second-order accurate FV-FE scheme will be presented.

Figures (18)

• Figure 1: Sketch of the location of the FEEC degrees of freedom: $V^0$ scalar variables are described by nodal valued functions (bottom-left); $V^1$ vector variables by edge-valued functions (bottom-right figures for the three spatial components, respectively); $V^2$ vector variables by face-valued functions (top-right figures for the three components, respectively); $V^3$ scalar variables by barycentric-valued functions (top-left figure). Note and remember that the nodal-point corresponds to the dual barycenters of the corresponding dual-grid.
• Figure 2: Sketch of the main (blue) and the barycentric dual (green) grids.
• Figure 3: Exact and numerical solutions for the Riemann problem RP0 ($\Delta x= 1/100$). The low Mach isolated, stationary contact wave is exactly preserved even at long time $t_f=1000$ for our hybrid FV-FEEC method ($s(\lambda^v)$). The FV-FEEC solution obtained by setting a higher numerical diffusion, i.e. $s=s(\lambda^b)$ or $s=s(\lambda^{\text{MHD}})$, is also plotted at time $t=0.01$. (see colored version online)
• Figure 4: Exact and numerical solution for the MHD Riemann problem RP1 at $t=0.1$. Density (left) and horizontal velocity $u_x$ (right) are plotted, comparing the numerical solutions obtained after choosing different time-step scales, i.e. the convection scale $\Delta t(\lambda^{v})$ (red squares), and the MHD scale $\Delta t(\lambda^{\text{MHD}})$ (blue deltas). (see colored version online)
• Figure 5: Exact and numerical solution for the Riemann problem RP2 at $t=0.2$. Density (left column) and magnetic field $B_y$ and $|B|$ (right column) are plotted, comparing the numerical solutions obtained after choosing different time-step scales, i.e. the convection scale $\Delta t(\lambda^{v})$ with $\theta_b=1$ (red squares), $\theta_b=0.55$ (green squares), and the MHD scale $\Delta t(\lambda^{\text{MHD}})$ (blue deltas). The vertical component $B_y$ is plotted next to the magnitude $|B|$ to show the correct approximation of the two rotational waves. (see colored version online)