Support operator method

Motivation

The Support Operator Method (SOM) is a mimetic finite difference technique used for discretizing operators that satisfy integral identities by construction. This approach is not only advantageous for ensuring conservation properties, but it has also been shown to improve numerical accuracy, especially in highly anisotropic problems, by reducing spurious perpendicular numerical diffusion.

For a general introduction to the Support Operator Method, we recommend the book Conservative Finite-Difference Methods on General Grids by M. Shashkov (CRC Press, Boca Raton, 1996). The benefits of the SOM in reducing spurious numerical diffusion on non-aligned grids are discussed in S. Guenter et al., J. Comput. Phys. 209:354, (2005). The application of the SOM within the FCI framework was first introduced in A. Stegmeir et al., Comput Phys. Commun. 198:139, (2016) and has since been refined in several follow-up publications. We have arrived at the toroidally staggered version with a refined discrete parallel gradient, allowing for more accurate handling of strong map distortions.

Support operator method

In addition to the canonical (3D) grid, denoted by integer toroidal grid indices (\ k=0,1,2,\dots), referred to as $G$, a toroidally staggered dual grid $G^∗$ with half-integral toroidal grid indices $k=\frac{1}{2},\frac{3}{2},\frac{5}{2},\dots$ is introduced. We can come up with a discrete version for the parallel gradient as a sparse matrix $\mathbf{Q}$, which maps from the canonical grid to the staggered grid: $$\begin{align} \nabla_\parallel\Rightarrow\mathbf{Q}, \quad \mathbf{Q}: G\rightarrow G^*. \end{align}$$ The SOM provides a method for constructing the discrete analogue of the parallel divergence, represented by the matrix $\mathbf{P}$, which maps vice versa from the staggered grid to the canonical grid: $$\begin{align} \nabla\cdot\left(\mathbf{b}\circ\right)\Rightarrow\mathbf{P}, \quad \mathbf{P}: G^*\rightarrow G, \end{align}$$

To construct the discrete parallel divergence, we begin by utilizing the following integral identity: $$\begin{align} \int\limits_V u \nabla\cdot\left(\mathbf{b}\nabla_\parallel v\right) dV = -\int\limits_V\nabla_\parallel u \nabla_\parallel v dV + \int\limits_{\partial V} \dots, \end{align}$$ where we omit the surface term, as it is related to the boundary conditions, which are handled separately. Rather than directly discretizing the parallel divergence, we discretize the integral identity, with the matrix $\mathbf{P}$ still unknown: $$\begin{align} \sum\limits_{\alpha,\beta,\gamma} u_\alpha\mathbf{P}_{\alpha\beta}\mathbf{Q}_{\beta\gamma}v_\gamma \, \Omega_\alpha = -\sum\limits_{\mu,\nu,\sigma}\mathbf{Q}_{\mu\nu}u_\nu\mathbf{Q}_{\mu\sigma}v_\sigma\, \omega_\mu. \end{align}$$ Here, Greek indices indicate summation over the full 3D grid or dual grid. The left-hand side refers to a discrete integration over the canonical grid $G$, where $\Omega_\alpha$ represents finite, locally field-aligned flux box volumes for the grid points (these become essentially grid cells). The right-hand side corresponds to an integration over the staggered grid $G^*$, with the flux box volumes denoted by $\omega_\mu$. These flux box volumes can be accurately computed during the field line trace procedure using the toroidal flux expansion factor. For additional details on this and an illustrative depiction of the flux boxes, see the parallel gradient section.

Since the above expression holds for arbitrary fields $u$ and $v$, we can relabel the indices and arrive at the following result: $$\begin{align} \mathbf{P}_{\alpha\beta} = -\mathbf{Q}_{\beta\alpha}\frac{\omega_\beta}{\Omega_\alpha}, \end{align}$$ or, by interpreting $\Omega$ and $\omega$ as diagonal matrices, we obtain the compact form: $$\begin{align} \mathbf{P} = -\Omega^{-1}\mathbf{Q}^T\omega \end{align}$$

The integral equality above implies the adjointness relation $\nabla\cdot\left(\mathbf{b}\circ\right)=-\nabla_\parallel^\dagger$, which is preserved at the discrete level with the correspondence $^\dagger\rightarrow ^T$ up to volume correction factors. As a result, the discrete parallel diffusion operator is self-adjoint by construction, ensuring the conservation of the $L_1$​ norm (e.g. energy, particles) in time-dependent problems.