# Multi-Region OpenFOAM Solver Development for Compact Toroid Transport in Drift Tube
## 3 Program Verification  
### The mhdF Solver Verification  
#### 3.1.1 Orszag-Tang MHD Vortex Problem  
Here, the classical Orszag-Tang MHD vortex problem is selected as a benchmark case to verify the mhdF solver. This case is a classic two-dimensional numerical test problem in MHD, widely used to verify the accuracy of MHD solvers and their ability to simulate the coupling effects between magnetic fields and fluid dynamics. It demonstrates the nonlinear interaction between magnetic fields and fluid motion, especially complex phenomena such as magnetic field stretching, compression, shock formation, and energy dissipation. This problem does not have an analytical solution but can be verified by comparing the results with those of other standard MHD solvers. After non-dimensionalization, the computational domain is \(x,\ y\in[0\ 1]^{2}\), and the initial conditions for magnetic flux density and velocity are as follows:  
\[B_{0}=[-\sin(2\pi y),\ \sin(4\pi x)]/\gamma \tag{15}\]\[U_{0}=[-\sin(2\pi y),\ \sin(2\pi x)] \tag{16}\]  
where \(\gamma=5/3\) is the adiabatic index; the initial values of pressure and density are \(P_{0}=0.6\), \(\rho_{0}=1.0\); and magnetic permeability is \(\mu=1\). Here, we compare the results with those computed by the ideal MHD solver of the FLASH code [30]. The simulation results are as follows. Figure 1 displays the distribution of magnitude of the magnetic flux density and density calculated by the mhdF solver at \(t=0.4\). The results exhibit high structural consistency with those from the FLASH code [30], demonstrating identical complex vortex flow patterns. Figure 2 compares the results calculated by the mhdF solver and the FLASH code at \(t=0.4\) along the line \(y=0.25\), showing the variations in the magnetic flux density components (\(B_{x}\), \(B_{y}\)) and the density along the x-axis. The blue solid lines represent the results calculated by the mhdF solver, while the red dashed lines correspond to the results from the FLASH code. The overall computational results are largely consistent, with some discrepancies at peak values. The PIMPLE algorithm in the mhdF solver explicitly treats the pressure term, resulting in insufficient resolution for transonic shock discontinuities. In contrast, FLASH uses a high-order Godunov method, which better captures shock waves. Nevertheless, for supersonic problems, the calculation has already met the requirements for simulating the dynamic processes of CTs in the drift tube to a certain extent. Figure 1: Displays the magnitude of the magnetic flux density and density distribution of the Orszag-Tang MHD vortex calculated by the mhdF solver at \(t=0.4\). Figure 2: The variations in three physical quantities along the \(x\)-axis at \(y=0.25\) and \(t=0.4\). (**a**) \(B_{x}\), (**b**) \(B_{y}\), (**c**) density. Here, the convergence order of the solver mhdF is discussed through relative errors. For the numerical results of the Orszag-Tang vortex simulation on an \(N\times N\) grid at time \(t\), the relative error of each variable is defined as follows:  
\[\delta_{N}(W)=\frac{\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{N}\left|W_{i,k}^{N}- W_{i,k}^{\max}\right|}{\sum\limits_{i=1}^{N}\sum\limits_{k=1}^{N}\left|W_{i,k}^{ \max}\right|} \tag{17}\]  
Here, \(W_{i,k}^{\max}\) represents the simulated value at the coordinate \((i,k)\) on the high-resolution grid. The value of \(\max\) is taken as 512. The combined relative error for multiple variables is defined as follows:  
\[\delta_{N}=\frac{1}{4}\left(\delta_{N}(U_{x})+\delta_{N}(U_{y})+\delta_{N}(B_{ x})+\delta_{N}(B_{y})\right) \tag{18}\]  
The corresponding formula for calculating the order of convergence is  
\[R_{N}=-\frac{\lg(\delta_{N}/\delta_{N_{0}})}{\lg(N/N_{0})} \tag{19}\]  
Here, let \(N_{0}=50\), which refers to the number of grid points in the coarsest grid, and \(N_{0}<N\).