[latexpage]
Problem definition
The computation of magnetohydrodynamic (MHD) equilibrium of toroidal plasmas plays a central role in tokamak fusion studies, not only as a fundamental tool for the computation of transport and stability calculations but also for the development of advanced tokamak control systems. Under perfect axisymmetric conditions, the vector calculus formulation of the Grad-Shafranov equation is the one most commonly used:
\begin{equation}
-\nabla\cdot\left(\frac{1}{\mu r}\nabla\psi\right) = r\frac{\mathrm{d}p}{\mathrm{d}\psi} + \frac{1}{\mu_{0}r}f\frac{\mathrm{d}f}{\mathrm{d}\psi} := j_{t}\;. \label{eq:grad_shafranov}
\end{equation}
This equation can be written as a system of three equations:
\begin{equation}
\begin{cases}\nabla\psi = \vec{u} \\ \nabla\cdot\vec{h} = \tilde{j}_{t}\end{cases} \mathrm{and}\quad\mathbb{K}\vec{u} = \vec{h}\quad\mathrm{with}\quad\mathbb{K}:=\left[\begin{array}{cc}\frac{1}{\mu_{0}r} & 0 \\ 0 & \frac{1}{\mu_{0}r}\end{array}\right].\label{eq::gs_system}
\end{equation}
The two equations on the left define differential (topological) relations and the third one establishes a constitutive relation between $\vec{u}$ (the gradient of $\psi$) and the flux vector $\vec{h}$. The direct discretization of this formulation using standard techniques presents two challenges:
- Conservation of current at the discrete level ($\int_{\Omega}\nabla\cdot\vec{h}\,\mathrm{d}\Omega=\int_{\Omega}j_{t}\mathrm{d}\Omega$).
- Extension to curved meshes.
The first one is an important requirement for coupling with transport equations, the second enables optimal high order approximation in curved domains.
Differential Geometry Formulation
In this work, we propose to use the differential geometric formulation of this equation:
\begin{equation}
\mathrm{d}\star_{\mathbb{K}}\mathrm{d}\star\tilde{\psi}^{(2)} = \tilde{j}_{t}^{(2)}\,.
\end{equation}
This corresponds to two topological relations (left) and two constitutive relations (right):
\begin{equation}
\begin{cases}\mathrm{d}\psi^{(0)} = u^{(1)} \\ \mathrm{d}\tilde{h}^{(1)} = \tilde{j}_{t}^{(2)}\end{cases} \!\!\!\!\mathrm{and}\,\, \begin{cases}\star\tilde{\psi}^{(2)} = \psi^{(0)} \\ \mathbb{K}\star u^{(1)} = \tilde{h}^{(1)}\end{cases}\!\!\!\!\!,\!\!\!\label{eq::gs_system}
\end{equation}
where each physical quantity is either an inner, $\Lambda^{k}$, or outer, $\tilde{\Lambda}^{k}$, oriented a $k$-form:
\[
\psi^{(0)}\in\Lambda^{0},\quad u^{(1)}\in\Lambda^{1}, \quad \tilde{h}^{(1)}\in\tilde{\Lambda}^{1},\quad \tilde{j}_{t}^{(2)},\tilde{\psi}^{(2)}\in\tilde{\Lambda}^{2}.
\]
Here $\mathrm{d}$ is the exterior derivative and $\star$ is the Hodge-$\star$ operator that transforms $k$-forms into $(k+1)$-forms and changes their orientation from inner to outer and vice-versa:
\[
\star:\Lambda^{k}\rightarrow\tilde{\Lambda}^{k+1} \quad\mathrm{or}\quad\star:\tilde{\Lambda}^{k}\rightarrow\Lambda^{k+1}.
\]
With these three operators, $(\mathrm{d},\star,\mathbb{K})$, it is possible to construct the following double de Rham complex
\begin{equation}
\begin{matrix}
\Lambda^{0}(\Omega)
&\stackrel{\mathrm{d}}{\longrightarrow}& \Lambda^{1}(\Omega)
&\stackrel{\mathrm{d}}{\longrightarrow}\; &\Lambda^{2}(\Omega) \; \\
\star\updownarrow & & \mathbb{K}\star\downarrow\uparrow\star\mathbb{K}^{-1} &
&\star\updownarrow & \\
{\tilde{\Lambda}}^{2}(\Omega)
&\stackrel{\mathrm{d}}{\longleftarrow}& {\tilde{\Lambda}}^{1}(\Omega)
&\stackrel{\mathrm{d}}{\longleftarrow}\; &{\tilde{\Lambda}}^{0}(\Omega)
\end{matrix} \label{eq::diffGeom_deRham_Complex}
\end{equation}
These spaces of inner and outer oriented $k$-forms are intimately related to the fact that $k$-forms naturally integrate over $k$-dimensional domains, $\Omega_{k}$. For example, $h^{(1)}$ is a flux quantity and therefore it naturally integrates across lines (since these equations are defined in a two dimensional space).
This formal classification of physical quantities into inner and outer oriented $k$-forms, although seeming a mere mathematical artefact, has deep impacts in the formulation of the discretization of the equations. This association is important because a physical system generally preserves integral quantities and these quantities are topological in the sense that they do not depend on a particular coordinate system or metric employed. Integration of differential forms is independent of any metric notions in space; no Riemann metric, dot products or $k$-dimensional volume elements are required. Integrals are of paramount importance in mimetic discretizations, because we cannot exactly represent differential forms and the operations on and between differential forms in a finite dimensional setting. On the other hand we can represent integrals of differential forms and the integral relations between them exactly.
Using the de Rham complex, (\ref{eq::diffGeom_deRham_Complex}), the inner product between inner oriented $k$-forms ($\alpha^{(k)},\beta^{(k)}\in\Lambda^{k}$) becomes
\begin{align*}
\left(\alpha^{(i)},\beta^{(i)}\right) := \int_{\Omega}\alpha^{(i)}\wedge\star\beta^{(i)}, \quad i=\{0,2\}\quad\mathrm{and}\quad\left(\alpha^{(1)},\beta^{(1)}\right) := \int_{\Omega}\alpha^{(1)}\wedge\mathbb{K}\star\beta^{(1)}\,,
\end{align*}
and for outer oriented $k$-forms ($\tilde{\alpha}^{(k)},\tilde{\beta}^{(k)}\in\tilde{\Lambda}^{k}$) yields
\begin{align*}
\left(\tilde{\alpha}^{(i)},\tilde{\beta}^{(i)}\right) := \int_{\Omega}\tilde{\alpha}^{(i)}\wedge\star\tilde{\beta}^{(i)}, \quad i=\{0,2\} \quad\mathrm{and}\quad
\left(\tilde{\alpha}^{(1)},\tilde{\beta}^{(1)}\right) := \int_{\Omega}\tilde{\alpha}^{(1)}\wedge\star\mathbb{K}^{-1}\tilde{\beta}^{(1)}\,.
\end{align*}
This formulation particularly contrasts with the one obtained with traditional vector calculus for the inner product between 1-forms (vectors in the vector calculus formulation):
\[
\left(\vec{a},\vec{b}\right) := \int_{\Omega}\vec{a}\cdot\vec{b}\,\mathrm{d}\Omega\,.
\]
It is possible to note that neither $\mathbb{K}$ nor $\mathbb{K}^{-1}$ appear in the vector calculus inner product. This results in a different finite element formulation, as will be seen below.
Mimetic Discretization
The mimetic spectral element method used in this work follows directly from the framework introduced in [Palha et al. 2014]. The key ideas of this method are: (i) special basis functions capable of representing the integral character of $k$-forms and (ii) a weak formulation that takes into account the relations established in the de Rham complex, (\ref{eq::diffGeom_deRham_Complex}).
These basis functions, $\tilde{\epsilon}^{(k)}_{i}$, satisfy the following properties:
\[
\tilde{\epsilon}^{(0)}_{i}(\tilde{\tau}_{(0),j})=\delta_{ij}, \quad \int_{\tilde{\tau}_{(1),j}}\tilde{\epsilon}_{i}^{(1)} = \delta_{ij}\quad\mathrm{and}\quad \int_{\tilde{\tau}_{(2),j}}\tilde{\epsilon}_{i}^{(2)} = \delta_{ij},
\]
where $\tilde{\tau}_{(0),j}$, $\tilde{\tau}_{(1),j}$ and $\tilde{\tau}_{(2),j}$ are the points, edges and surfaces of the spectral element mesh.
The weak problem of (\ref{eq::gs_system}), then reduces to finding the approximate solution, $\tilde{\psi}_{h}\in\tilde{\Lambda}^{2}$, such that:
\[
\begin{cases}
\left(\mathbb{K}\star\mathrm{d}\star\tilde{\psi}^{(2)}_{h},\tilde{\epsilon}^{(1)}_{j}\right) = \left(\tilde{h}_{h}^{(1)},\tilde{\epsilon}_{j}^{(1)}\right) \quad \forall \tilde{\epsilon}_{j}^{(1)}\\
\left(\mathrm{d}\tilde{h}_{h}^{(1)},\tilde{\epsilon}_{j}^{(2)}\right) = \left(\tilde{j}_{t}^{(2)},\tilde{\epsilon}_{j}^{(2)}\right) \quad \forall \tilde{\epsilon}_{j}^{(2)}.
\end{cases}
\]
Using integration by parts we can rewrite this as:
\[
\begin{cases}
\int_{\partial\Omega}\tilde{\psi}_{h}^{(2)}\wedge\tilde{\epsilon}_{j}^{(1)}-\left(\tilde{\psi}_{h}^{(2)},\mathrm{d}\tilde{\epsilon}^{(1)}_{j}\right) = \left(\tilde{h}_{h}^{(1)},\tilde{\epsilon}_{j}^{(1)}\right) \quad \forall \tilde{\epsilon}_{j}^{(1)} \\
\left(\mathrm{d}\tilde{h}_{h}^{(1)},\tilde{\epsilon}_{j}^{(2)}\right) = \left(\tilde{j}_{t}^{(2)},\tilde{\epsilon}_{j}^{(2)}\right) \quad \forall \tilde{\epsilon}_{j}^{(2)}.
\end{cases}
\]
In doing this, the discrete analog of the flux operator ($\mathbb{K}\star\mathrm{d}\star$) is the negative adjoint of the discrete divergence ($\mathrm{d}$).
Results
This formulation was applied to a Soloviev solution of the Grad-Shafranov equation, [Howell et al. 2014], corresponding to:
\[
\psi(r,z) = \frac{f_{0}R_{0}^{2}a^{2}}{2}\left[1-\left(\frac{z}{a}\right)^{2} – \left(\frac{r-R_{0}}{a} + \frac{(r-R_{0})^{2}}{2aR_{0}}\right)^{2}\right],\quad p(\psi)=\frac{f_{0}}{\mu_{0}}\psi + p_{0},\quad f\frac{\mathrm{d}f}{\mathrm{d}\psi} = f_{0}R_{0}^{2},
\]
with $a = 0.5$, $f_{0} = 1.0$ and $R_{0}=1.1$. To assess this method on meshes of different curvature, the following space transformation that maps the space $(\xi,\eta)\in[-1,1]\times[-1,1]$ into a curved domain $(x,y)\in[x_{a},x_{b}]\times[y_{a},y_{b}]$ is used:
\[
\begin{cases}
x(\xi,\eta) = \left(\xi + c\sin(\pi\xi)\sin(\pi\eta) + 1\right)\frac{x_{b}-x_{a}}{2} \\
y(\xi,\eta) = \left(\eta + c\sin(\pi\xi)\sin(\pi\eta) + 1\right)\frac{y_{b}-y_{a}}{2}
\end{cases}.
\]
An example of these results, together with the error, are presented below.
With this formulation, optimal $h$- and $p$-convergence are obtained regardless of the curvature of the mesh, as can be observed in the following figures.
This shows the potential of this method for the accurate solution of the Grad-Shafranov equation on curved meshes.
[fruitful_sep]
References
Palha, A., Koren. B. and Felici, F. (2016). A mimetic spectral element solver for the Grad-Shafranov equation. Journal of Computational Physics, accepted. doi:10.1016/j.jcp.2016.04.002
Howell, E. C. and Sovinec, C. R. (2014). Solving the Grad–Shafranov equation with spectral elements. Computer Physics Communications, 185(5), 1415–1421. doi:10.1016/j.cpc.2014.02.008
Palha, A., Rebelo, P. P., Hiemstra, R., Kreeft, J. and Gerritsma, M. (2014). Physics-compatible discretization techniques on single and dual grids, with application to the Poisson equation of volume forms. Journal of Computational Physics, 257, 1394–1422. doi:10.1016/j.jcp.2013.08.005