Overview

The Magnetohydrodynamic Algorithm outside a Sphere (MAS) code integrates the time-dependent resistive thermodynamic magnetohydrodynamic (MHD) equations in three-dimentional spherical coordinates and is used extensively in models of coronal structure [Mikić and Linker (1996), Mikić et al. (1999), Linker et al. (1999), Lionello et al. (2008), Downs et al. (2013) , Mikić et al. (2018)], coronal dynamics [Lionello et al. (2003), Lionello et al. (2006), Linker et al. (2011)], and coronal mass ejections out to the Earth [Linker et al. (2003), Lionello et al. (2013), Török et al. (2018)]. It is the primary MHD model in CORHEL (Corona-Heliosphere), a suite of models for describing the solar corona and inner heliosphere [Riley et al. (2012)] and CORHEL-CME, an interface-based system for modeling CME eruptions from the Sun to Earth, both of which are available for public runs-on-demand at NASA's Community Coordinated Modeling Center.

To develop a realistic model, photospheric magnetic field observations are used to specify the boundary condition on the radial magnetic field (Br). This data is taken from measurements at the Wilcox Solar Observatory, the Mount Wilson Observatory, the GONG and SOLIS instruments from the National Solar Observatory, the SOI/MDI instrument aboard the SOHO spacecraft, the HMI instrument aboard the SDO spacecraft, or the PHI instrument aboard the SolO spacecraft. The resulting solution provides a three-dimensional description of the coronal and/or heliospheric plasma and magnetic fields.

To compare the results with observations, we calculate synthetic observables that can be quantitatively compared with ground and space-based observations. For example, from the simulated electron density, we compute the total and polarized brightness of the K-Corona in white light. The white light corona is best observed on the ground during total solar eclipses, but it is also observed by coronographs, such as the HAO Mauna Loa Coronameter, and in space by the LASCO coronagraph aboard SOHO, the SECCHI COR1, and COR2 coronagraphs and HI heliospheric imager aboard STEREO, and the WISPR instrument aboard Parker Solar Probe. Using the model density and temperature, we also synthesize EUV and soft X-ray emission as observed in space, from instruments like EIT aboard SOHO, EUVI aboard STEREO, XRT aboard Hinode, and AIA aboard SDO. Direct comparisons of plasma and magnetic field properties can be made with in-situ measurements from spacecraft such as OMNI, STEREO, and Parker Solar Probe. We also visualize structural and topological features of the magnetic field, such as open field regions and the squashing factor "Q" [Titov (2007)].

The MAS Model

The MHD equations used in MAS are \[ \begin{alignat}{2} \frac{\partial {\bf A}}{\partial t} &= {\bf v} \times {\bf B} -{\frac{c^2}{4\,\pi}\,{\textcolor{lightblue}{\eta}}\,\nabla \times \nabla \times {\bf A}} \\ \frac{\partial \rho}{\partial t} &=-\nabla\cdot(\rho\,{\bf v}) \\ \frac{\partial T}{\partial t} &=-\nabla\cdot(T\,{\bf v})- (\gamma - 2)\left(T\,\nabla\cdot{\bf v}\right)+\left.\frac{(\gamma-1)}{2\,k}\frac{m_p}{\rho}\right[\left.- \nabla\cdot\left({\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{q_1}}}} + {\bf q_2}\right) -\frac{\rho^2}{m_p^2}\,{\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{Q}}}} + {\bf H}\right] \notag \\ \rho\,\frac{\partial{\bf v}}{\partial t}& = -\rho\,{\bf v}\cdot \nabla\,{\bf v} + \frac{1}{c}\,{\bf J} \times {\bf B} - \nabla (p+p_w) + \rho\,{\bf g} + {\bf F_{c}} + \nabla\cdot({\textcolor{lightblue}{\nu}}\,\rho\,\nabla {\bf v}) {\textcolor{lightblue}{+ \nabla\cdot\left({\bf S}\,\rho\,\nabla{\frac{\partial \bf v}{\partial t}}\right)}} \notag \end{alignat} \]
where ${\bf A}$ is the magnetic vector potential, $\rho$ is the plasma density, $T$ is the temperature, ${\bf v}$ is the plasma velocity, ${\bf B}=\nabla\times {\bf A}$ is the magnetic field which satisfies the $\nabla \cdot {\bf B}=0$ condition, ${\bf J}=(c/4\pi)\,\nabla\times {\bf B}$ is the current density, $p=2\,k\,T\,\rho/m_p$ is the plasma pressure, $p_w$ is the Alfvén wave pressure (see below), ${\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{q_1}}}}$ is the modified heat flux due to parallel Spitzer collisional thermal conduction (see below), ${\bf q_2}$ is the heat flux due to collisionless thermal conduction (see below), ${\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{Q}}}}$ is the modified optically thin radiative loss function (see below), ${\bf H}$ is the coronal heating (see below), $\gamma$ is the adiabatic index, ${\bf g}=-g_0\,R_{\odot}^2\,{\bf \hat r}/r^2$ is the gravitational force, ${\bf F_c}=\rho\,\left(2\,\Omega \times {\bf v} + \Omega \times \left(\Omega \times {\bf r}\right)\right)$ combines the Coriolis and centrifugal terms used when computing in the co-rotating frame, $c$ is the speed of light in a vacuum, $m_p$ is the proton mass, $k$ is Boltzman's constant, $\eta$ is the resistivity coefficient, $\nu$ is the kinematic viscosity coefficient, and ${\bf S}$ is a semi-implicit factor (see below). Terms colored light blue are added/modified for numerical efficiency and/or stability (see below).

The modified parallel Spitzer collisonal heat flux ${\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{q_1}}}}$ [Spitzer and Harm (1953)] and collisionless heat flux ${\bf q_2}$ [Hollweg (1976)] are defined as \[ {\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{q_1}}}} =f_{\mbox{c}}(r)\,\textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}\,\kappa_0\,T^{5/2}\, {\bf \hat b}{\bf \hat b}\cdot\nabla T, \qquad {\bf q_2}=(1-f_{\mbox{c}}(r))\,\frac{\alpha_{\mbox{nc}}\,k\,\rho}{(\gamma-1)\,m_p}\,T\,({\bf v}-\Omega\times {\bf r})\,{\bf \hat b}{\bf \hat b} \] where ${\bf \hat b}=|{\bf B}|/{\bf B}$ is the normalized direction of the magnetic field, $\kappa_0=9\times10^{-7}$ is the Spitzer parallel electron thermal conductivity constant, $\alpha_{\mbox{nc}}\sim O(1)$ is based on the deviation from Maxwellian of the solar wind's electron distribution, $f_{\mbox{c}}(r)=1/2-1/2\,\mbox{tanh}\left((r-r_{\mbox{tc}})/\sigma_{\mbox{tc}}\right)$ is a radial profile used to transition between collisional and collisionless thermal conduction, and $\textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}$ is a modification function used for resolving the transition region (see below).

The radiative loss function in MAS is \[ {\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{Q}}}}={\bf Q}/\textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}, \] where ${\bf Q}$ is a selected temperature-dependent radiative loss function and $\textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}$ is a modification function used for resolving the transition region (see below). MAS contains multiple options for specifying the radiative loss function including [Rosner et al. (1978), Mok et al. (2004), Athay (1986)], as well as some directly based on the CHIANTI database. The CHIANTI options currently use version 7.1 [Landi et al. (2013)] through lookup tables based on coronal, photospheric, or hybrid abundances.

MAS contains numerous coronal heating models, both empirical and physics-based, that can be combined together. Empirical heating can be set by a combination of analytic functions [e.g. Lionello et al. (2009)], while the primary physics-based approach is the wave-turbulence-driven (WTD) model [Lionello et al. (2014), Downs et al. (2016), Mikić et al. (2018)]. In the WTD approach, auxiliary equations are solved to describe the large-scale evolution of turbulent Alfvénic fluctuations in the corona and solar wind. This takes the form of two coupled PDEs governing the evolution of the envelope of the fluctuation amplitudes ($z_{\pm}$) via propagation, reflection, and dissipation: \[ \frac{\partial z_{\pm}}{\partial t} =-\left({\bf v}\pm {v_A}\right)\cdot\nabla\,z_{\pm} + \frac{z_{\pm}}{4}\left({\bf v} \mp {v_A}\right)\cdot\nabla\left(\mbox{ln}\,\rho\right) + \frac{z_{\mp}}{2}\left({\bf v} \mp {v_A}\right)\cdot\nabla\left(\mbox{ln}\left|{v_A}\right|\right) - \frac{z_{\pm}\left|z_{\mp}\right|}{2\,\lambda_{\perp}}, \] where the Alfvén wave speed is defined as $v_A^2=|{\bf B}|^2/(4\,\pi\,\rho)$ and the Poynting flux of waves entering the domain is chosen by setting a uniform outward wave amplitude $z_0$ at the lower boundary ($\left|z_{\pm}(r=R_{\odot})\right|=z_0$).
The WTD heating term takes the form of: \[ H_{\mbox{wtd}}=\frac{\rho}{4\,\lambda_{\perp}}\,\left(|z_{-}|\,z_{+}^2 + |z_{+}|\,z_{-}^2\right) \] where $\lambda_{\perp}=\lambda_0\,\sqrt{{B_w}/{|{\bf B}|}}$ is the correlation length, and $B_w=6.09\,\mbox{G}$ is a chosen reference magnetic field and $\lambda_0$ is a tunable parameter.

MAS can also solve for the Alfvén wave pressure using a more simplified approach than WTD. In this case, it is found by averaging the forward and backward Alfvén wave pressures as $p_w=\left(\epsilon_{+}+\epsilon_{-}\right)/2$ where $\epsilon{\pm}$ are computed by integrating the Wentzel-Kramers-Brillouin (WKB) approximation [Mikić et al. (1999)] for the wave energies: \[ \frac{\partial \epsilon_{\pm}}{\partial t} =-\nabla\cdot \left( \epsilon_{\pm}\left[{\bf v}\pm v_A\,{\bf \hat b}\right]\right) - \frac{1}{2}\,\epsilon_{\pm}\,\nabla\cdot {\bf v} \]

The MAS code uses an adaptive time-step that conforms to the plasma flow CFL condition. As this is often a much slower time-scale than the magnetosonic wave travel time, a “semi-implicit ” term is added to the velocity advance to provide unconditional linear stability for time steps larger than the fast magnetosonic wave limit.
The term “semi-implicit ” is used in this context to refer to a term that modifies the dispersion relations for MHD waves in a manner similar to a fully implicit algorithm [Schnack et al. 1987] and has also been described as the approximate implicit approach [Carmana et. al. (1991)]. The semi-implicit term is computationally advantageous compared to a fully-implicit approach because it results in a symmetric positive definite operator (easier to solve) and can be solved separately (operator split) from other terms in the MHD equations.
The semi-implicit term is computed locally on the grid based on the Alfvén and sound speeds, and in MAS is set as [Lionello et al. (1999)]: \[ \textcolor{lightblue}{{\bf S(r,\theta,\phi)}=\frac{1}{\Delta t^2\,\tilde k^2}\,\left(\frac{C_w^2}{(1-C_f)^2} - 1\right)} \] where $\Delta t$ is the time-step, $\tilde k = 2\,\left(\Delta r^{-2} + (r\,\Delta\theta)^{-2}+(r\,\Delta \phi\,\sin\theta)^{-2}\right)^{1/2}$, $C_f=\Delta t\,\tilde k \cdot {\bf v}$ is the flow CFL, and $C_w=(1/2)\,\Delta t\,\tilde k\,(v_c^2+v_A^2)^{1/2}$ is the wave CFL, where $v_c=(\gamma\,p/\rho)^{1/2}$ is the sound speed.
The semi-implicit term goes to zero in regions that are explicitly stable for the selected time step, as well as in the limit of a steady-state solution ($d/dt\rightarrow 0$).

Although stable, the wave-like terms do not guarantee monotonicity, resulting in the potential emergence of oscillatory behavior on the order of the grid size. In order to avoid such behavior (and to prevent physical structures smaller than the grid from trying to form), the dissipative resistivity $\textcolor{lightblue}{\eta}$ and kinematic viscosity $\textcolor{lightblue}{\nu}$ coefficients are set artificially higher than their physical values. However, to maintain accuracy, they are chosen such that the resulting diffusion time scales remain much larger than the typical Alfvén travel time. For flexibility, these coefficients can be constant or spatially-dependent.

Lastly, resolving the transition region between the lower chromosphere and the corona poses a difficult challenge for numerical models. Here, the temperature and density change by orders of magnitude over very little distance. Even with a non-uniform grid, resolving this region requires an exorbitant number of grid points. To alleviate this problem, we apply a smooth modification function to the Spitzer thermal conduction coefficient and its inverse to the radiative loss function [Lionello et al. (2009), Mikić et al. (2013)]: \[ {\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{q_1}}}} = \textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}\,{\bf q_1}, \] \[ {\bf \textcolor{lightblue}{\widetilde{\textcolor{white}{Q}}}} = {\bf Q}/\textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}, \] where \[ \textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T) = \left(1+(T/T_{\mbox{cut}})^{-10}\right)^{1/4}} \] such that $\textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}\rightarrow 1$ for $T \gt T_{\mbox{cut}}$ and $\textcolor{lightblue}{f_{\mbox{$\kappa$-mod}}(T)}\rightarrow (T/T_{\mbox{cut}})^{-5/2}$ for $T \lt T_{\mbox{cut}}$ (acting as a temperature cut-off function for ${\bf q_1}$ at temperatures below $T_{\mbox{cut}}$). This causes the transition region to broaden substantially while preserving the coronal solution above.

Numerical Methods

The MAS code is written in Fortran and parallelized for CPUs and GPUs using Message Passing Interface (MPI) , OpenACC, and Fortran's `do concurrent' standard parallelism [Caplan et al. (2023)]. It solves the MHD equations on a non-uniform logically-rectangular staggered grid using finite differences. The non-uniformity of the grid helps MAS to efficiently resolve small-scale structures such as the transition region and active regions, while allowing for coarser grid points over the global scale. Each component of the field vectors are staggered which, combined with the use of the vector potential as a primitive, guarantees that the magnetic field is divergence-free. Advective terms are spatially differenced using first-order upwinding, while parabolic and gradient terms are centrally differenced using `Method A' from [Veldman and Rinzema (1992)].

Several time integration schemes are used in the model. Advective terms are integrated with an explicit predictor-corrector upwinding scheme [Lionello et al. (1999)] where any “reaction ” terms (such as heating, radiative loss, gravitational forces etc.) are computed in the corrector step. The semi-implicit term is solved in the predictor-corrector framework implicitly using two preconditioned conjugate gradient (PCG) solves, while diffusive-like terms (such as thermal conduction and viscosity) are operator-split and either solved implicitly with backward-Euler+PCG or explicitly with a second-order Runge-Kutta Gegenbauer (RKG2) super time-stepping scheme combined with a PTL time step limiting method [Caplan et al. (2024)]. The WKB Alfvén wave pressure advance is operator split and integrated with the forward-Euler scheme sub-cycled at the Alfvén wave CFL limit. The wave turbulence advance is also operator split from the main MHD equations and sub-cycled, but uses a 2nd-order TVD Runge-Kutta time-stepping scheme [Gottlieb and Shu (1998)]. Within each step, we use Strang splitting of the various operators. A OSPRE flux-limiting scheme [Waterson and Deconinck (2007)] is used for the advective operators, while a point-implicit scheme is used for the remaining operators [Tóth et al. (2012)].

The PCG solver in MAS has two communication-free preconditioning options. Depending on the problem, we either use a non-overlapping domain decomposition with zero-fill incomplete LU factorization [Saad (2003)] or point-Jacobi (diagonal scaling). The operator matrices are stored in a modified diagonal (DIA) sparse storage format for internal grid points, while the boundary conditions are applied matrix-free. The preconditioner matrix is stored in a compressed sparse row (CSR) storage format optimized for vectorized memory access [Smith and Zhang].

The code exhibits parallel strong scaling up to thousands of CPU cores (MPI ranks) on large HPC systems. A 3D decomposition is used to distribute to domain across MPI ranks in a topologically-aware communicator. Inter-rank communication is performed asynchronously whenever possible. A hybrid MPI+STDPAR+OpenACC approach is used to scale the code to multiple GPU compute nodes [Caplan et al. (2023)].