Finding magnetic surfaces

Two-dimensional data $\Psi (R, Z)$ on a rectangular grids $(R, Z)$ is read from the G_EQDSK file (G-file) of EFIT code. Based on the 2D array data, I use 2D cubic spline interpolation to construct a interpolating function $\Psi = \Psi (R, Z)$. To construct a magnetic surface coordinate system, we need to find the contours of $\Psi$, i.e., magnetic surfaces. The values of $\Psi$ on the magnetic axis, $\Psi_0$, and the value of $\Psi$ on the last closed flux surface (LCFS), $\Psi_b$, are given in G-file. Using these two values, I construct a 1D array ``psival'' with value of elements changing uniform from $\Psi_0$ to $\Psi_b$. Then I try to find the contours of $\Psi$ with contour level value ranging from $\Psi_0$ to $\Psi_b$. This is done in the following way: construct a series of straight line (in the poloidal plane) that starts from the location of the magnetic axis and ends at one of the points on the LCFS. Combine the straight line equation, $Z = Z (R)$, with the interpolating function $\Psi (R, Z)$, we obtain a one variable function $h = \Psi (R, Z (R))$. Then finding the location where $\Psi$ is equal to a specified value $\Psi_i$, is reduced to finding the root of the equation $\Psi (R, Z (R)) - \Psi_i = 0$. Since this is a one variable equation, the root can be easily found by using simple root finding scheme, such as bisection and Newton's method (bisection method is used in GTAW code). After finding the the roots for each value in the array ``psival'' on each straight lines, the process of finding the contours of $\Psi$ is finished. The contours of $\Psi$ found this way are plotted in Fig. 8.

Figure 8: Verification of the numerical code that calculates the contours of the poloidal flux $\Psi$. The bold line in the figure indicates the LCFS. The contour lines (solid lines) given by the gnuplot program agrees well with the results I calculate by using interpolation and root-finding method (the two sets of contours are indistinguishable in this scale). My code only calculate the contour lines within the LCFS, while those given by gnuplot contains additional contour lines below the X point and on the left top in the figure. Eqdisk file of the equilibrium was provided by Dr. Guoqiang Li (filename: g013606.07104).

In the above, we mentioned that the point of magnetic axis and points on the LCFS are needed to construct the straight lines. In G-file, points on LCFS are given explicitly in an array. The location of magnetic axis is also explicitly given in G-file. It is obvious that some of the straight lines $Z = Z (R)$ that pass through the location of magnetic axis and points on the LCFS will have very large or even infinite slope. On these lines, finding the accurate root of the equation $\Psi (R, Z (R)) - \Psi_i = 0$ is difficult or even impossible. The way to avoid this situation is obvious: switch to use function $R = R (Z)$ instead of $Z = Z (R)$ when the slope of $Z = Z (R)$ is large (the switch condition I used is $\vert d Z / d R\vert > 1$).

In constructing the flux surface coordinate with desired Jacobian, we will need the absolute value of the gradient of $\Psi$, $\vert \nabla \Psi \vert$, on some specified spatial points. To achieve this, we need to construct a interpolating function for $\vert \nabla \Psi \vert$. The $\vert \nabla \Psi \vert$ can be written as

$$\vert \nabla \Psi \vert = \sqrt{\left( \frac{\partial \Psi}{\partial R} \right)^2 + \left( \frac{\partial \Psi}{\partial Z} \right)^2},$$ (206)

By using the center difference scheme to evaluate the partial derivatives with respect to $R$ and $Z$ in the above equation (using one side difference scheme for the points on the rectangular boundary), we can obtain an 2D array for the value of $\vert \nabla \Psi \vert$ on the rectangular $(R, Z)$ grids. Using this 2D array, we can construct an interpolating function for $\nabla \Psi$ by using the cubic spline interpolation scheme.

Figure 9: Grid points (the intersecting points of two curves in the figure) corresponding to uniform poloidal flux and uniform poloidal arc length for EAST equilibrium shot 13606 at 7.1s (left) (G-file name: g013606.07104) and shot 38300 at 3.9s (right) (G-file name: g038300.03900).

Figure 10: $\vert \nabla \Psi \vert$ as a function of the poloidal angle. The different lines corresponds to the values of $\vert \nabla \Psi \vert$ on different magnetic surfaces. The stars correspond to the values on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis). The equilibrium is for EAST shot 38300 at 3.9s.

Figure 11: The Poloidal magnetic field $B_p = \vert \nabla \Psi \vert / R$ (left) and toroidal magnetic field $B_{\phi } = g / R$ (right) as a function of the poloidal angle. The different lines corresponds to the values on different magnetic surfaces. The stars correspond to the values on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis). The equilibrium is for EAST shot 38300 at 3.9s.

Figure 12: Equal-arc Jacobian as a function of the poloidal angle on different magnetic surfaces. The dotted line corresponds to the values of Jacobian on the boundary magnetic surface. The equilibrium is for EAST shot 38300 at 3.9s.

Figure 13: $\vert \nabla \Psi \vert$ as a function of the poloidal angle. The different lines corresponds to values of $\vert \nabla \Psi \vert$ on different magnetic surfaces. The stars correspond to the values of $\vert \nabla \Psi \vert$ on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis). The equilibrium is a Solovev equilibrium.

Figure 14: Jacobian on different magnetic surfaces as a function of the poloidal angle. The equilibrium is a Solovev equilibrium and the Jacobian is an equal-arc Jacobian. The stars correspond to the values of Jacobian on the boundary magnetic surface while the plus signs correspond to the value on the innermost magnetic surface (the magnetic surface adjacent to the magnetic axis).