WARNING: still work in progress
Giacomo is working on this
Spatial Discretization in JOREK
The discretization of JOREK is well described in the papers O. Czarny, G. Huysmans, J.Comput.Phys 227, 7423 (2008) and S. Pamela et al. J. Comput. Phys. 464, 111101 (2022). Here we present a summary of the most important aspects and the code implementation.
As extensively explained here, JOREK uses the cylindrical coordinates $(R,Z,\phi)$. The discretization on the poloidal plane (variables $R$ and $Z$) is discretized using 2D Bezier finite elements, discussed the following section whereas the discretization along the toroidal direction (variable $\phi$) is performed with a truncated real Fourier series, as explained in this section.
The two discretizations combined read as follows:
\[\mathbf{X}(R,Z,\phi) = \left(\sum_{l=1}^{N_{tor}} Z_l(\phi) \left (\sum_{i = 1}^{N_{vert}}\sum_{j = 1}^{\text{dof_per_vertex}} h^{ij} \vec{u}^{ij} b^*_{i,j}(R,Z)\right) \right ) \tag{1}\]The inner parenthesis is the poloidal discretization, which uses the basis functions $b^*_i$, while $Z_l$ are the basis functions for the toroidal discretization.
2D Bezier finite elements in the poloidal plane
JOREK uses isoparametric Gm-continuous Bezier finite elements for spatial discretization on the poloidal plane $RZ$. Instead of using control points as "direct" degrees of freedom, a different formulation, called nodal formulation (see here), is used, which greatly facilitate imposing Gm-continuity across multiple finite elements. What isoparametric and Gm-continuous (with $m\in \mathbb{N}$) means will be explained soon.
For an overview of Bezier curves, Bezier surfaces and how to pass from these to the concept of Bezier finite elements see the Appendix A
We use the parameters $s,t\in [0,1]$ as local coordinates of a Bezier finite element, so the discretization on a single Bezier finite element of degree $n$ reads as
\[X_\text{pol}(s,t) = \sum_{i=0}^{n}\sum_{j=0}^{n} \mathbf{P}_{ij} B_i^n(s)B_j^n(t) \tag{2}\]with $B_i^n$ and $B_j^n$ Bernstein polynomials and $P_{ij}$ control points of the surface.
As explained in the following section, Equation $(2)$ is rewritten in a more convenient way, that is:
\[X_\text{pol}(s,t) = \sum_{i=0}^{N_\text{vert}-1} \sum_{j=0}^{\text{dof_per_vertex}} h^{ij} \vec{u}^{ij} b_i(s,t)\]It is important to note that this formulation holds only locally to one finite element.
Formally it can be extended to multiple elements, creating a mapping between $(R,Z)$ and $(s,t)$, which is needed anyway. … (finish here)
Bezier finite element nodal representation
In JOREK a particular formulation of the finite elements is used. In particular, the control points are divided into $4$ groups, based on which interpolation point they are close to. When the degree $n$ is odd, each group has $(n+1)^2/4$ points. To each interpolation point (each node), a basis of vectors (the actual degrees of freedom!) is assigned and the control points are rewritten as a combination of this vectors.
In the following of this section we provide an example with bi-cubic Bezier finite element, which is the most used finite element in JOREK. Next we discuss about the generalized nodal representation to any order, not only bi-cubic, that is implemented in JOREK and we show that $G_m$ continuity is simpler to be imposed on such a representation. Finally we provide an insight on why this particular formulation is chosen and is meaningful.
Nodal representation of bi-cubic Bezier finite element
With bi-cubic Bezier finite element the Bernstein polynomials' basis is a tensor product of the cubic Bernstein polynomials' basis.
The control points are divided in 4 groups as shown in the following illustration.
\[\begin{matrix} 30 &&& 31 &&& 21 &&& 20 \\ \bullet &&& \bullet &&& \bullet &&& \bullet \\ \\ 32 &&& 33 &&& 23 &&& 22 \\ \bullet &&& \bullet &&& \bullet &&& \bullet \\ \\ 02 &&& 03 &&& 13 &&& 12 \\ \bullet &&& \bullet &&& \bullet &&& \bullet \\ \\ 00 &&& 01 &&& 11 &&& 10 \\ \bullet &&& \bullet &&& \bullet &&& \bullet \\ \\ \end{matrix}\]Points $x0$, with $x\in {0,1,2,3}$ are the interpolation points and the points in their group are $xy$ with $y\in {1,2,3}$.
Then for each group $i$ the vectors $p_i, v_i, u_i$ and $w_i$ are introduced along with certain pre-defined quantities $d_{u,i}$ and $d_{v,i}$
\[\begin{align} \mathbf{P}_{0,0} &= p_1 \\ \mathbf{P}_{0,1} &= p_1 + v_1\, d_{v,1} \\ \mathbf{P}_{1,0} &= p_1 + u_1\, d_{u,1} \\ \mathbf{P}_{1,1} &= P_{0,1} + P_{1,0} - p_1 + w_1\, d_{u,1}\, d_{v,1} \\ \\ \mathbf{P}_{3,0} &= p_2 \\ \mathbf{P}_{3,1} &= p_2 + v_2\, d_{v,2} \\ \mathbf{P}_{2,0} &= p_2 + u_2\, d_{u,2} \\ \mathbf{P}_{2,1} &= P_{3,1} + P_{2,0} - p_2 + w_2\, d_{u,2}\, d_{v,2} \\ \\ \mathbf{P}_{3,3} &= p_3 \\ \mathbf{P}_{3,2} &= p_3 + v_3\, d_{v,3} \\ \mathbf{P}_{2,3} &= p_3 + u_3\, d_{u,3} \\ \mathbf{P}_{2,2} &= P_{2,3} + P_{3,2} - p_3 + w_3\, d_{u,3}\, d_{v,3} \\ \\ \mathbf{P}_{0,3} &= p_4 \\ \mathbf{P}_{0,2} &= p_4 + v_4\, d_{v,4} \\ \mathbf{P}_{1,3} &= p_4 + u_4\, d_{u,4} \\ \mathbf{P}_{1,2} &= P_{0,2} + P_{1,3} - p_4 + w_4\, d_{u,4}\, d_{v,4} \end{align} \tag{4}\]This formulation is based on the following work: O. Czarny, G. Huysmans, J.Comput.Phys 227, 7423 (2008).
Generalized nodal representation of Bezier finite element
The work S. Pamela et al. J. Comput. Phys. 464, 111101 (2022) generalizes what has been shown in the previous paragraph for bi-cubic elements to any degree $n$ and usually imposing $G_m$ continuity with $m=(n-1)/2$. Instead of using $p_i,u_i,v_i,w_i$ vectors, the notation is $u_{ij}$ and $d_{u,i}$ with $h^{ij}$ For example, in the bi-cubic case we have $u_{i0}:=p_i, u_{i1}:=u_i, …$. Then the generalized formulation, only for group of node $\mathbf{P}_{00}$, reads as follows:
\[\mathbf{P}_{ij} = h^{ij}u^{ij} + \sum_{k=0}^i\sum_{l=0}^i (-1)^{1+i+j+k+l}(1-\delta_{ki}\delta_{lj}) \begin{pmatrix} i \\ k \end{pmatrix} \begin{pmatrix} j \\ l \end{pmatrix} \mathbf{P}_{kl}\]with $0\leq i,j \leq (n+1)/2$
Imposing the continuity on generalized nodal representation
In JOREK one node, i.e. $\mathbf{P}_{00}$ is always shared by $4$ finite elements. Let $\xi_{1,1}, \xi_{-1,1}, \xi_{-1,-1}, \xi_{1,-1}$ be the finite elements that share the node, listed counter clockwise. Then each of these elements has a different set of vectors $u^{ij}$ (and values $h^{ij}$) as nodal representation used in tha shared node $\mathbf{P}_{00}$. We say that $u^{ij}$ and $h^{ij}$ are the nodal representation used by element $\xi_{1,1}$, $u^{-ij}$ and $h^{-ij}$ for element $\xi_{-1,1}$, $u^{i-j}$ and $h^{i-j}$ for element $\xi_{1,-1}$ and finally $u^{-i-j}$ and $h^{-i-j}$ for element $\xi_{-1,-1}$.
With the introduced notation, $G_m$ continuity on the shared nodes is obtained by imposing the following constrains.
$\forall j, \ \ h^{-ij}$ is constrained by:
\[h^{-ij} = \begin{cases} -\alpha h^{ij} & \text{for } i = 1 \text{ and } \alpha > 0 \\ h^{ij} & for i \neq 1 \end{cases}\]$\forall i, \ \ h^{i-j}$ is constrained by:
\[h^{i-j} = \begin{cases} -\beta h^{ij} & \text{for } j = 1 \text{ and } \beta > 0 \\ h^{ij} & for j \neq 1 \end{cases}\]Brief intuition on the formulation
As demonstrated in Corollary 1 of S. Pamela et al. J. Comput. Phys. 464, 111101 (2022), in the particular nodal formulation used, the vectors $u^{ij}$ are equal in direction and sign (but not absolute value) with the derivative
\[\frac{\partial^{i+j} \mathbf{P}}{\partial^{i} s \partial^{j} t}(s=0,t=0)\]This is extremely helpful, for instance, when building the mesh, since derivatives of $R$ and $Z$ are used to obtain a flux aligned grid.
Discretization with nodal formulation
With the nodal representation,
Implementation in the code
The following is still to be adapted
Two-dimensional third-order Bernstein polynomials $B_{i,j}^{3}(s,t)$ defined by
\[B_{i,j}^{3}(s,t) = B_{i}^{3}(s)B_{j}^{3}(t)\qquad i,j = 0\ldots 3 \tag{1}\]with
\[B_{i}^{3}(s) = \frac{3!}{i!(3 - i)!} s^{i}(1 - s)^{3 - i} \tag{2}\]are used for the discretization in the poloidal plane. The coordinates $s = 0\ldots 1$ and $t = 0\ldots 1$ form the element-local coordinate system. A quantity $X$ which may be a coordinate or a physical variable $(R,Z,\Psi,T,\ldots$; iso-parametric discretization) can be expressed by
\[X(s,t) = \sum_{i = 0}^{3}P_{i,j}B_{i,j}^{3}(s,t). \tag{3}\]As first order continuity of the finite elements is demanded, not all combinations of control points $P_{i,j}$ are valid. Effectively, four free parameters $p_k$, $u_k$, $v_k$, and $w_k$ remain per node $k$ and quantity $X$ of the element. The control points $P_{i,j}$ can be reconstructed from these free parameters as given in the following. Nodes $k = 1\ldots 4$ correspond to the positions $(i,j) = (0,0)$, $(3,0)$, $(3,3)$, and $(0,3)$.
\[\begin{array}{r l} P_{0,0} &= p_1 \\ P_{0,1} &= p_1 + v_1\, d_{v,1} \\ P_{1,0} &= p_1 + u_1\, d_{u,1} \\ P_{1,1} &= P_{0,1} + P_{1,0} - p_1 + w_1\, d_{u,1}\, d_{v,1} \\ P_{3,0} &= p_2 \\ P_{3,1} &= p_2 + v_2\, d_{v,2} \\ P_{2,0} &= p_2 + u_2\, d_{u,2} \\ P_{2,1} &= P_{3,1} + P_{2,0} - p_2 + w_2\, d_{u,2}\, d_{v,2} \\ P_{3,3} &= p_3 \\ P_{3,2} &= p_3 + v_3\, d_{v,3} \\ P_{2,3} &= p_3 + u_3\, d_{u,3} \\ P_{2,2} &= P_{2,3} + P_{3,2} - p_3 + w_3\, d_{u,3}\, d_{v,3} \\ P_{0,3} &= p_4 \\ P_{0,2} &= p_4 + v_4\, d_{v,4} \\ P_{1,3} &= p_4 + u_4\, d_{u,4} \\ P_{1,2} &= P_{0,2} + P_{1,3} - p_4 + w_4\, d_{u,4}\, d_{v,4} \end{array} \tag{4}\]Here, the $d_{u,k}$ and $d_{v,k}$ are measures for the distances of the control points from the element nodes. These are element properties, i.e., for a node that belongs to several elements, the value of $d_{u,k}$ depends on the element considered. In JOREK, they are represented by
\[\begin{align} &\texttt{element%size(k,1)} \rightarrow 1 \\ &\texttt{element%size(k,2)} \rightarrow d_{u,k} \\ &\texttt{element%size(k,3)} \rightarrow d_{v,k} \\ &\texttt{element%size(k,4)} \rightarrow d_{u,k} d_{v,k} \end{align}\]Projecting on the bezier elements
The distribution $X(s,t)$ of quantity $X$ within a certain element can be expressed by Equation $(3)$. Inserting Equations $(1)$-$(2)$ and Equation $(4)$ leads to the following expression
\[X(s,t) = \sum_{k = 1}^{4} \tilde{p}_k(s,t) \tag{5}\]with the contributions
\[\tilde{p}_k = b_{k,1}\, p_k + b_{k,2}\, u_k\, d_{u,k} + b_{k,3}\, v_k\, d_{v,k} + b_{k,4}\, w_k\, d_{u,k}\, d_{v,k} \tag{6}\]associated to the four nodes. Here,
\[\begin{array}{r l} b_{1,1} &= (1 - s)^{2} (1 - t)^{2} (1 + 2 s) (1 + 2 t) \\ b_{1,2} &= 3 (1 - s)^{2} (1 - t)^{2} s (1 + 2 t) \\ b_{1,3} &= 3 (1 - s)^{2} (1 - t)^{2} (1 + 2 s) t \\ b_{1,4} &= 9 (1 - s)^{2} (1 - t)^{2} s t \\ b_{2,1} &= s^{2} (1 - t)^{2} (3 - 2 s) (1 + 2 t) \\ b_{2,2} &= 3 s^{2} (1 - t)^{2} (1 - s) (1 + 2 t) \\ b_{2,3} &= 3 s^{2} (1 - t)^{2} (3 - 2 s) t \\ b_{2,4} &= 9 s^{2} (1 - t)^{2} (1 - s) t \\ b_{3,1} &= s^{2} t^{2} (3 - 2 s) (3 - 2 t) \\ b_{3,2} &= 3 s^{2} t^{2} (1 - s) (3 - 2 t) \\ b_{3,3} &= 3 s^{2} t^{2} (3 - 2 s) (1 - t) \\ b_{3,4} &= 9 s^{2} t^{2} (1 - s) (1 - t) \\ b_{4,1} &= (1 - s)^{2} t^{2} (1 + 2 s) (3 - 2 t) \\ b_{4,2} &= 3 (1 - s)^{2} t^{2} s (3 - 2 t) \\ b_{4,3} &= 3 (1 - s)^{2} t^{2} (1 + 2 s) (1 - t) \\ b_{4,4} &= 9 (1 - s)^{2} t^{2} s (1 - t) . \end{array} \tag{7}\]In JOREK, the quantities $b_{k,j}$ are required at the position of the Gaussian points with indices $(i_G, j_G)$. They are denoted
$\texttt{H(k,j,i_G,j_G)} b_{k,j}$ at Gaussian point $(i_G,j_G) \tag{8}$
and are initialized in the subroutine ''basisfunctions2''.
Notation in the code
The degrees of freedom of each node corresponding to the $l$-th toroidal harmonic of the $\nu$-th physical quantity are denoted in the following way in the code:
\[\begin{align} &\texttt{node%values(l,1,\nu)} \rightarrow p_k \\ &\texttt{node%values(l,2,\nu)} \rightarrow u_k \\ &\texttt{node%values(l,3,\nu)} \rightarrow v_k \\ &\texttt{node%values(l,4,\nu)} \rightarrow w_k \end{align}\]The degrees of freedom of the coordinates $R$ ($\mu = 1$) and $Z$ ($\mu = 2$) are called,
\[\begin{align} &\texttt{node%x(1,\mu)} \rightarrow p_k \\ &\texttt{node%x(2,\mu)} \rightarrow u_k \\ &\texttt{node%x(3,\mu)} \rightarrow v_k \\ &\texttt{node%x(4,\mu)} \rightarrow w_k \end{align}\]For all $\nu = 1\ldots N_{var}$ and $i_G, j_G = 1\ldots N_{Gauss}$ and $p = 1\ldots N_{plane}$, the variable values at a given Gaussian point $(s_{i_G}, t_{j_G})$ at the toroidal position $\phi_p$ in a given finite element can be expressed by
\[\begin{aligned} X_{\nu}(s_{i_G},t_{j_G},\phi_{p}) &= \sum_{k = 1}^{N_{vert}}\sum_{j = 1}^{N_{ord}}\sum_{l = 1}^{N_{tor}} \mathrm{nodes}(i)\%\mathrm{values}(1,j,\nu) \\ &\qquad \cdot \mathrm{H}(k,j,i_G,j_G) \cdot \mathrm{element}\%\mathrm{size}(k,j) \cdot \mathrm{HZ}(1,p), \end{aligned} \tag{9}\]where i=element%vertex(k). Here, $N_{var}$ denotes the number of physical variables, $N_{Gauss} = 4$ the number of Gaussian points used for Gauss quadrature in $s$ and $t$ directions each, $N_{plane} = 2(N_{tor} - 1)$ the number of toroidal planes located at $\phi_{p} = 2\pi (p - 1) / N_{plane}$, $N_{vert} = 4$ the number of vertices in each element, $N_{ord} = 4$ the number of degrees of freedom per vertex, and $N_{tor}$ the number of different toroidal Fourier modes. The quantity ''HZ(l,p)'' corresponds to the value of the $l$-th Fourier mode at the toroidal position $\phi_{p}$ and is denoted $Z_{l}(\phi_{p})$ in the following. The table below lists which Fourier modes correspond to the different mode indices $l$.
| l | 1 | 2 | 3 | 4 | 5 | … |
|---|---|---|---|---|---|---|
$Z_l(\phi_p) \equiv$ HZ(l,p) | 1 | $\cos \phi_p$ | $\sin \phi_p$ | $\cos 2\phi_p$ | $\sin 2\phi_p$ | … |
Real Fourier series in toroidal direction
- For the toroidal direction, a real Fourier series ($\cos$, $\sin$) is used:
| JOREK harmonic | Toroidal mode number | Toroidal basis function |
|---|---|---|
| 1 | 0 | 1 |
| 2 | n_period | cos(n_period*phi) |
| 3 | n_period | sin(n_period*phi) |
| 4 | 2*n_period | cos(2n_periodphi) |
| 5 | 2*n_period | sin(2n_periodphi) |
| … | … | … |
| n_tor-1 | (n_tor-1)/2*n_period | cos((n_tor-1)/2n_periodphi) |
| n_tor | (n_tor-1)/2*n_period | sin((n_tor-1)/2n_periodphi) |
- The hard-coded_parameters
n_torandn_periodare used to select the harmonics included in the simulation.n_tor: Total number of real Fourier modes (odd integer number)n_period: Toroidal periodicity (positive integer number)
- The maximum toroidal mode number included in a simulation is
n_max=(n_tor-1)/2*n_period. - A sufficient number of toroidal planes is required to avoid aliasing. The minimal requirement is typically
n_plane >= 2 * n_tor. It is always a good idea to scann_planefor a simulation to ensure convergence.n_planemust be a power of 2, if FFTW is not used. If FFTW is used (see here for more),n_planecan be an arbitrary positive integer number, but not for all values FFTW is similarly efficient.
Appendix A: Introduction on Bezier curves and Bezier surfaces
We start with a brief overview of Bezier curves and then move from this background to the particular formalism used in JOREK.
Bezier curves
Bezier curves are parametric curves $\mathbf{P}(t)$ which are controlled by a set of points, intuitively called control points. The parameter $t$ is in the range $[0,1]$. The first and the last points, namely $\mathbf{P}_0$ and $\mathbf{P}_{n - 1}$ on the sequence of control points are interpolation points, that is $\mathbf{P}(0) = \mathbf{P}_0$ and $\mathbf{P}(1) = \mathbf{P}_1$. The influence of the control points on the shape of the curve is determined by the a base ${B_0(t), \dots, B_n(t)}$, called Bernstein polynomials. Point i is associated to the polynomial $\mathbf{P}_i$, indeed the mathematical expression of the Bezier curve is
Bernstein polynomial are defined by the following formula:
\[B_{i}^{n}(s) = \frac{n!}{i!(n - i)!} s^{i}(1 - s)^{n - i}\qquad i = 0\ldots n\]Bernstein polynomials' basis is fundamental in the Bezier finite element framework. In JOREK, it is not the basis local to the single finite element, but the latter is derived from the Bernstein basis. More on this in the section Bezier finite element nodal representation.
When referring to Bezier curves of degree $n$ we mean that we have $n + 1$ points and $n + 1$ polynomials of degree $n$.
Intuition on Bezier curves: Bezier curves were invented by the French engineer Pierre Bézier and the intuitive way of building them is through the de Casteljau's algorithm. You can see Bezier curves as Linear intERPolation (in graphics vocabulary, LERP) applied recursively. First, you linearly interpolate with parameter
tevery two consecutivive control points, so you generatenpoints. Then you do the same with this new "derived" control point, again with parametert, until you reach only one point. This is $\mathbf{P}(t)$ The following image depicts what just described:
Important note on dimensionality: It is important to underline that points $\mathbf{P}_i$ can belong to a space $\mathbb{R}^m$ with any $m$, given that $m > 1$. Bezier curves are generally though with $\mathbf{P}_i \in \mathbb{R}^2$ but it my be as well $\mathbf{P}_i \in \mathbb{R}^{100}$! Indeed in JOREK the smallest space is $\mathbb{R}^8$ since at least $6$ variables are used.
Bezier surfaces
Bezier surfaces can be viewed as the tensor product of Bezier curves. The basis has now the following form:
\[\{B_{i}(s)B_{j}(t)\}_{i = 0\ldots n,j = 0\ldots n}\]and to each polynomial there is the corrspective $\mathbf{P}_{i,j}$ element. The Bezier surface then writes:
\[\mathbf{P}(s,t) = \sum_{i = 0}^{n}\sum_{j = 0}^{n}B_{i}^{n}(s)B_{j}^{n}(t)\]
From Bezier surfaces to Bezier finite elements
In JOREK Bezier surfaces are used for the discretization on the poloidal plane. The points $\mathbf{P}_i$ has as first two coordinates the domain coordinates $R$ and $Z$, whereas the remaining $6+$ coordinates correspond to the variables (functions) of interest (see here for the list of variables used).
In its most elementary version a Bezier finite element is nothing but a Bezier surface and the degrees of freedom are the points $\mathbf{P}_i$. To be precise, each entry of a point $\mathbf{P}_i \in \mathbb{R}^m$ is a degree of freedom, so the total number of degrees of freedom is $n \cdot m$.
The formulation gets more complex when patching together multiple finite elements to create the global representation. Indeed constraints on the continuity has to be added. We will give now an example with $G_0$ continuity, for a more in depth explanation of $G_0$ and $G_1$ look at O. Czarny, G. Huysmans, J.Comput.Phys 227, 7423 (2008) and for $G_2$ and a generalization for every $G_m$ look at S. Pamela et al. J. Comput. Phys. 464, 111101 (2022).
Imposing $G_0$ continuity
In JOREK all the finite elements have $4$ interpolation points and they have $4$ neighbouring finite elements (except on the boundary or on particular points such as the center of the grid). Imposing $G_0$ continuity between two neighbouring finite elements means that all the control points on the "interfacing" edge, that is the 2 shared interpolation points plus the control points along the edge that connects those two interpolation points, must be equal.
For example, with bi-cubic Bezier finite elements, $4$ control points (2 interpolation points and 2 non interpolation points) reside on the edge and has to be imposed equal to the $4$ points of the other finite element.
From imposing continuity to the nodal representation
Imposing high continuity involves a great number of control points and so there is not a clear vision of what are the degrees of freedom and what is instead fixed.
The nodal representation, discussed in this section, is meant to greatly simplify this, bringing clarity on what are the degrees of freedom. See the aforementioned section for more details on this.