The SPIDER Code - Axisymmetric Fixed Boundary Plasma Equilibrium Solver
Preprint, Inst. Appl. Math., the Russian Academy of Science

Вычислительный код SPIDER - расчет аксиально-симметричного равновесия плазмы с фиксированной границей

          Contents

          1.   Introduction    …………………………………………. 4

2.      Mathematical model   ………………………………… 6

3.      Finite-difference scheme   ……………………………. 10

4.      Iteration process   …………………………………….. 12

5.  Test simulations   ……………………………………… 16

6.  Examples of runs and Figures……………………… 17

7.   References    …………………………………………... 24

1.  Introduction

Tokamak plasma equilibrium computation is a fundamental problem of magnetic confinement studies. Many plasma processes, including linear and early nonlinear stages of magneto-hydrodynamic (MHD) instabilities, plasma evolution and transport, plasma flows, waves and turbulence, represent different kinds of deviations from MHD equilibrium. Thus they require accurate calculations of plasma equilibrium configurations (see, for example, Refs. [1-8]).

          In this paper we describe an axisymmetric fixed boundary plasma equilibrium solver - the SPIDER code - for solving the tokamak equilibrium problem based on the nonlinear Grad-Shafranov equilibrium equation with different sets of prescribed profiles, e.g.  and ,   and ,  and . Here  is the plasma pressure,  is the poloidal current,  is the safety factor,  is the averaged parallel to the magnetic field component of the current density,  is the poloidal flux. The SPIDER code output consists of the magnetic surfaces coordinates and other equilibrium magnetic field characteristics.

          Because of the nonlinearity of the Grad-Shafranov equation all numerical methods for equilibrium calculation are iterative. Conventionally they can be subdivided into two classes:

a.)               Eulerian, that use a prescribed (e.g., rectangular or conformal to the plasma boundary) mesh calculation as, for example, in Refs. [8-14];

b.)               Lagrangian, that use curvilinear flux coordinates and adaptive to magnetic surfaces mesh for equilibrium calculations as in Refs. [7], [15-20], [26-35].

          Eulerian methods have an advantage of being able to easily reproduce the two-dimensional geometry of complicated configurations as in Ref. [10]. They are widely used for simulations of equilibrium control in tokamaks as in Refs. [21], [36-38] and for the interpretation of experimental magnetic measurements (e.g. Refs. [22-24]). The disadvantage of those methods is in their limited range of formulations of equilibrium problems, restricted essentially to a prescribed right hand side of the Grad-Shafranov equation. Formulations requiring other input profiles are difficult to implement due to a necessity of frequent magnetic surface coordinate computing and averaging over them, as discussed in Ref. [25].

          Lagrangian methods have an advantage in permitting to solve  wide range of problems. Any set of two one-dimensional functions of the radial flux coordinate, that uniquely determines the current density and the pressure profile, can be used in flux coordinate methods. A disadvantage of flux coordinate methods consists in difficulties  treating free boundary plasma equilibrium, especially with a separatrix (see, e.g., Ref. [7]).

          Methods, that use flux coordinates, in turn, can be subdivided into two types:

a.)               variational (Refs. [16], [19], [28-30]) and inverse coordinate (Refs. [17], [31-32], [35]) methods;

b.)              adaptive grid methods (Refs. [15], [18], [33, 34]).

Variational and inverse coordinate codes solve equations in which unknowns are flux coordinates themselves. Adaptive grid codes solve equilibrium equation on given curvilinear grid and then use the numerical solution for advancing the computational grid and adjustment of the flux coordinate system.

The theory of perturbed equilibrium approach for solving the Grad-Shafranov equation for different formulations of the tokamak equilibrium problem has been described in Ref. [39]. The discrepancy between nonlinear equation and its linear analog in that work brings in the main contribution. Such approach is particularly adequate for the adaptive grid codes (such as Ref. [15]) and can use constraints, which are specific to the each equilibrium problem. However, it is more complicated to implement this approach in the variational methods or inverse variable codes due to the difficulty in the lineralization of equilibrium equations used in those codes (although some attempts have been made, for example, in Ref. [16] with successful results presented in Ref. [35]).

In this paper we describe the SPIDER code - an axisymmetric fixed boundary plasma equilibrium solver, which is based on an adaptive grid approach.

 The main restriction of the code consists in an assumption of nested magnetic surfaces with a single magnetic axis in plasma.

The main achievements of the code are in the increased computational speed and high accuracy of resolution of both differential plasma equilibrium problem and its discrete model.

The SPIDER code was thoroughly tested by means of both analytical tests and comparison with such well-known equilibrium solvers as POLAR-2D ( KIAM, Ref. [32, 33] ), CAXE ( KIAM, Ref. [34] ) and EFIT ( GA, Ref. [40] ).

At the present moment the SPIDER code is used in such ITER Central Team codes as PET ( KIAM, Ref. [8, 38] ), DINA ( TRINITI, Ref. [36,37] ) and ASTRA (Kurchatov Institute, Ref. [41] ).

Despite of large number of tokamak fixed boundary plasma equilibrium codes and long history of their development, there is still a demand for accurate, fast and robust code, which can be used for extensive calculations of complicated toroidal plasma equilibrium problems: such as plasma equilibrium with high beta, high  elongation of the plasma boundary, low aspect ratio and large Shafranov shift, with x-point at plasma boundary, plasma with reverse shear current density profile, and nonzero current density at plasma boundary etc. Large number of tests and real device numerical simulations show (in authors opinion) that the SPIDER code meets all the above-mentioned requirements.

2.  Mathematical model

In the case of an axisymmetric plasma configuration in a magnetic plasma confinement device (tokamak type) the following vector equations:

(where- gas kinetics isotropic plasma pressure,  - magnetic field vector,  - plasma current density vector) can be reduced into a well known scalar Grad-Shafranov equation given by (2.1), (2.2) for poloidal flux function .

The SPIDER code is designed for numerical solving of the Grad- Shafranov axisymmetric fixed boundary plasma equilibrium problem

 ,   ,                   (2.1)

 ,                                 

where - cylindrical coordinate system associated with plasma configuration symmetry axis, - plasma domain in plane,  is the poloidal magnetic flux function, is the toroidal plasma current density:

 ,                        (2.2)

 ,

 ,

Here  is the plasma pressure profile, is the toroidal component of the magnetic field vector potential,  is the toroidal component of the magnetic field, is the poloidal magnetic flux (measured from symmetry axis ),  is the poloidal current (measured from the symmetry axis ):

                    ,     ,

where  is the poloidal surface, enclosed by circular line in toroidal direction passing through the  point.

          Assuming presence of the nested surfaces with a single magnetic axis we can define plasma toroidal flux and current in the following manner:

                     ,

  ,

where  is toroidal cross-section  of magnetic surface  .

The following boundary condition at the plasma boundary curve  is used in the SPIDER code:

                             ,       .                  (2.3)

It is assumed that plasma boundary is prescribed. There are two possibilities to set plasma boundary  in the SPIDER code:

1.     by means of the following analytical formulas:

   ,                                     (2.4)

where  are coordinates of the geometric center of the plasma boundary, - poloidal angle in  plane,  is the plasma boundary aspect ratio,  is the plasma boundary elongation and  is the plasma boundary triangularity:

                             ,                            (2.5)

where  are the upper and lower elongations respectively,  are the upper and lower triangularities respectively. Of course it is possible to use any other appropriate analytical formulas.

2.     by means of  prescribed boundary points coordinates array:   .

There are following settings for  axisymmetric fixed boundary equilibrium problem in the SPIDER code:

1. Setting with prescribed  - right hand side (2.2) of the Grad-Shafranov equation (2.1) and toroidal plasma current value  , i.e. with prescribed profiles  as functions of  normalized poloidal flux  and with  toroidal plasma current value given by

                                                          .                                (2.6)

There are two possibilities to set profiles in the SPIDER code:

a.)               by means of the analytical formulas:

        ,                                     (2.7)

       ,

Of course it is possible to use any other appropriate analytical formulas.

b.)              by means of prescribed arrays.

          2. Setting with prescribed  profiles  and safety factor  as functions of  normalized poloidal flux  and with prescribed value of poloidal flux in plasma . Inside plasma domain the poloidal current function is connected with safety factor  as follows:

where integral is taken along the magnetic surface  toroidal cross-section line.

There are two possibilities to set  profile for the problem setting 2 with the prescribed safety factor profile:

a.)            by means of analytical formula:

 ;

Of course it is possible to use any other appropriate analytical formulas.

b.)           by means of prescribed arrays:  .

There is another possibility to set  the profile of  in the SPIDER code by means of prescribed arrays:  .

          In any considered setting of the equilibrium problem there is  possibility to prescribe the value of poloidal beta  by means of rescaling the pressure profile .

 In the SPIDER code the value of poloidal beta  is defined in the following manner:

 .

          The total beta  is defined as:

,

where  - vacuum toroidal field value at the plasma center and

, 

where  - the plasma volume.

3.   Finite-difference scheme

The computational domain  is covered by a computational grid, which is topologically equivalent to a radially-annular grid and it is used as initial guess for construction of final magnetic surface adaptive grid.

The SPIDER code discrete model is based on magnetic surface adaptive grid finite-difference method with 9-points difference scheme.

A priori unknown adaptive grid has radially-annular structure:

                   ,   ,                      (3.1)

where  is the adaptive radial direction variable,  is the poloidal angle direction variable. Grid lines  form the set of nested magnetic surfaces where  labels magnetic axis location,  labels magnetic surface  and  labels plasma boundary curve . Grid lines  are chosen to be straight.

The difference analog of the Grad-Shafranov operator  (2.1) is constructed on the basis of the conservative finite-difference approximation of the  operator by means of the operator-variational scheme [42].

Let us consider a vector potential  of the poloidal magnetic field  with only one non-trivial toroidal projection, connected with the poloidal magnetic flux  as follows:

    ,                                               (3.2)

where  - the unit vector in toroidal angle  direction. Then Grad-Shafranov operator  (2.1) can be written in the form of the non-trivial projection of the  operator on toroidal angle  direction:

            ,                                                         (3.3)

 or in a scalar form

                     .                                                         (3.4)

          For any two dimensional cell of our 2D radially-annular grid we can construct corresponding three dimensional cell, which is a body of rotation of this 2D cell in the  direction.

          The difference approximation of the vector potential  is determined by means of the orthogonal projections on the edges of the 3D cell. In our case we have only non-trivial projections on the edges in toroidal direction:

                      .                                                                            (3.5)

          The difference approximation of the  operator in terms of its orthogonal projections on the unit normals to the faces of the 3D cell is constructed on the basis of the invariant definition, that is valid for an arbitrary infinitesimal area  with a boundary curve :

                   ,                                                         (3.6)

where  is the unit vector normal to the  face. A direct approximation of  (3.6) on the faces of the 3D cell gives the following approximation of the orthogonal projections of the :

,         (3.7)

,       

,

,

 where , ,  and  are corresponding lengths of the 2D cell edges. These projections are related to the centers of the corresponding faces of the 3D cell.

          In turn the difference approximation of the  operator is determined by means of the orthogonal projections on the edges of the 3D cell and is related to the their centers. For construction of the operator the difference analog of the following integral theorem is used:

       .         (3.8)

4. Iteration process

Let us denote finite difference approximation of the Grad-Shafranov equation (2.1) on the adaptive grid (3.1) by means of subscript

.                                      (4.1)

          The final purpose of solving of the discrete problem (4.1) with corresponding boundary conditions approximation consists in the finding both unknown coordinates of magnetic surface adaptive grid (3.1) and desired unknown discrete function  in the nodes of this grid.

Let us assume that the adaptive grid  coordinates and the corresponding numerical solution  for iteration number  are known. Then the iteration loop of the SPIDER code consist of the following two parts.

The first part implies that we solve the linear problem with respect to  on known adaptive grid :

.                          (4.2)

The second part implies that we construct new  approximation of the unknown adaptive grid  assuming that prescribed  grid lines coincide with  lines of closed magnetic surfaces.

Let us note that the matrix of the linear algebraic equations (4.2) depends on adaptive grid coordinates. Therefore, we need to invert this matrix on each iteration step . As the matrix inversion on each iteration step is very expensive, a combination of direct and iterative methods  is employed to solve the equations (4.2) more efficiently in course of adaptive grid iterations. With this in mind let us write

                          ,         .                            (4.3)

Because matrix  changes very slowly from iteration to iteration, the inverse matrix  is close to the matrix . Therefore we can write

       ,                  (4.4)

where matrix  can be used as a "preconditioner" for some fast iterative method.

          Finally instead of the matrix inversion on each iteration step we use the following iterative procedure:

a.)               for iteration step  we compute the inverse matrix and solve equations (4.2);

b.)              for the next  iteration steps equations (4.2) are solved by means of some fast iterative method using as a "preconditioner" the matrix in the right hand side of the equation (4.4):

          .

          In the SPIDER code the following variants of the grid adaptation variable  may be used:

1.       ,

2.           ,

3.          ,

where  - normalized poloidal flux,  - normalized toroidal flux.

          Let us consider in more detail iteration process for equilibrium problem with prescribed safety factor  profile.

          On the basis of averaged 1D Kruscal-Kulsrud equilibrium equation

                                                                              (4.5)

and relations between currents and fluxes in axisymmetric configuration

           ,   ,   ,   ,           (4.6)

we can rewrite averaged 1D equilibrium equation in the following form

                     .                                 (4.7)

Here coefficients  are determined only by means of magnetic surfaces geometry. Assuming that the derivative =  in (4.6) , we define the coefficients   with the help of the following formulae

,  ,                                     (4.8)

where  is the magnetic surface volume.

On the basis of 1D relation (4.7) the following iteration process works very well for solving fixed boundary tokamak equilibrium problem with prescribed safety factor  profile:

(1.)           let us know magnetic surface geometry and all required magnetic surface functions  on iteration number n;

(2.)           we can solve 1D equation (4.7) on "frozen"  1D flux grid

with prescribed boundary values  and .

          and obtain   and   ;

(3.)           we can make the Picard iteration for solving of the 2D  Grad-Shafranov equilibrium equation (2.1) on "frozen" 2D grid  (3.1)

and on the basis of this solution fulfill the construction of new approximation of 2D adaptive grid (i.e. magnetic surfaces coordinates);

(4.)    after step (3.) we can compute all required magnetic surface functions  on iteration number    and go to the beginning of the iteration loop.

          Let us consider in more detail iteration process for equilibrium problem with prescribed  profile.

Taking into account that

                                                                 (4.9)

and relations (4.5), (4.6) we can obtain the following set of two 1D nonlinear equations for unknowns  and :

                         ,

                                                                                                                   (4.10)

                       .

with the following boundary conditions:

           and   ,

where  - vacuum poloidal current value – is prescribed.

We can solve system of 1D equations (4.10) on "frozen" 1D flux grid to obtain  ,  and  taking into account (4.6) relations between currents and fluxes, and get .

As in the case of prescribed safety factor  profile, we can make the Picard iteration for solving of the 2D Grad-Shafranov equilibrium equation (2.1) on "frozen" 2D grid  (3.1)

 ,

and on the basis of this solution fulfill the construction of new approximation of 2D adaptive grid (i.e. magnetic surfaces coordinates).

5.   Test simulations

          Validation of the code against exact solution on the set of the different grid sizes and plasma parameters has been carried out.  For this test the following formulas of the exact solution were taken from  [43], p.132 :

c_{0  ,}  where  ,

 ,    ,

 ,    ,  .

Parameter  was chosen to satisfy  .

The results of simulations for different values of plasma aspect ratio , elongation and triangularity  are presented at the table № 1. Calculated values of coordinate of magnetic axis  and corresponding poloidal flux value  are compared to the exact values ,  for different values of grid size . It is seen that difference between numerical and exact solution  is quite small. It is necessary to point out, that accuracy of numerical solution convergence from grid size is quadratic.

Table № 1.

6.   Examples of runs and Figures

          Several examples of the SPIDER code runs are presented in this section.

The first case corresponds to high beta poloidal  and prescribed reverse shear averaged current density profile  with circular plasma boundary shape. Results of simulation are shown in Figs. 1-3. It is seen that magnetic axes is shifted far away from the plasma geometrical center.

In the second case, plasma equilibrium with a high plasma boundary elongation  is calculated. This case demonstrates that code can reliably calculate quite exotic equilibrium configuration without problems. Results of simulation are shown in Figs. 4-6. 

As an example for the third case, plasma boundary with a single X-point has been chosen. This example shows, that code can be used to simulate realistic tokamak experiments with diverted plasma. Results are shown in the Figs. 7-10.

Fig.1  Magnetic surfaces adaptive grid: , given  .

Fig.2  Plasma profiles versus normalized  : safety factor , poloidal current , current density profile parameters  and

Fig.3. Plasma profiles versus normalized : safety factor , poloidal current , current density profile parameters  and averaged .

Fig.4. Magnetic surfaces adaptive grid: plasma elongation

Fig.5. Plasma profiles versus normalized  : safety factor , poloidal current , current density profile parameters  and .

Fig.6. Plasma profiles versus normalized  : safety factor , poloidal current , current density profile parameters  and averaged  .

          Fig.7.  Magnetic surfaces adaptive grid: diverted plasma.

Fig.8.  Plasma toroidal current density: diverted plasma.

Fig.9. Plasma profiles versus normalized  : safety factor , poloidal current , current density profile parameters  and .

Fig.10. Plasma profiles versus normalized : safety factor , averaged magnetic field , averaged toroidal current density , poloidal current .

7.   References

[ 1 ]   V.D.Shafranov, Reviews of Plasma Physics (Consultants Bureau, New York, 1966), Vol. 2, p. 103.

[ 2 ]   H.Grad and J.Hogan, Phys. Rev. Lett. 24, 1337 (1970).

[ 3 ]   H.Grad, P.N.Hu and D.C.Stevens, Proc. Natl. Acad. Sci. USA 72, 3789 (1975).

[ 4 ]   L.E.Zakharov and V.D.Shafranov, Reviews of Plasma Physics (Consultants Bureau, New York, 1981), Vol. 11, p. 153.

[ 5 ]   R.R.Khairutdinov and V.E.Lukash, J. Comput. Phys. 109, 193 (1993).

[ 6 ]   L.Degtyarev, A.Martynov, S.Medevedev, F.Troyon, L.Vilard and R.Gruber, Comput. Phys. Commun, 103, 10 (1997).

[ 7 ]   S.A.Galkin, V.V.Drozdov, A.A.Ivanov, S.Yu.Medvedev and Yu.Yu.Poshekhonov, General Atomic Report GA-A23045, May 1999.

[ 8 ]   S.A.Galkin, A.A.Ivanov, S.Yu.Medvedev and Yu.Yu.Poshekhonov, Nucl. Fusion 37, No. 10, 1455 (1997).

[ 9 ]   W.Feneberg and K.Lackner, Nucl. Fusion 13, 549 (1973).

[10]    F.J.Helton and T.S.Wang, Nucl. Fusion 18, 1523 (1978).

[11]    J.L.Lohnson, H.E.Dalhed, J.M.Green et al., J. Comput. Phys. 32, 212, (1979).

[12]    J.P.Goedbloed, Comput. Phys. Commun, 24, 311 (1981).

[13]    J.Blum, J.L.Foil and B.Thooris, Comput. Phys. Commun. 24, 235 (1981).

[14]    H.Lutjens, A.Bondesons and A.Roy, Comput. Phys. Commun. 69, 287 (1992).

[15]    J.DeLucia, S.C.Jardin and A.M.Todd, J. Comput. Phys. 37, 183, (1980).

[16]    K.M.Ling and S.C.Jardin, J. Comput. Phys. 58, 300, (1985).

[17]    R.Gruber, R.Iacono and F.Troyon, J. Comput. Phys. 73, 168, (1987).

[18]    G.T.A.Huysmans, J.P.Goedbloed and W.Kerner, in Proceedings CP90 Conference on Computational Physics Proceedings (Word Scientific, Singapore, 1991), p. 371.

[19]    B.Turkington, A.Lifschitz, A.Eydeland and J.Spruck, J. Comput. Phys. 106, 269 (1993).

[20]    G.O.Ludwig, Plasma Phys. Controlled Fusion 39, 2021 (1997).

[21]    R.Albanese and F.Villone, Nucl. Fusion 38, 723 (1998).

[22]    L.L.Lao, H.E.St.John, R.D.Stambaugh, D.Stambaugh, A.G.Kellman and W.Pfeiffer, Nucl. Fusion 25, 1611 (1985).

[23]    J.Blum, E.Lazzaro, J.O'Rourke, B.Keegan and Y.Stephan, Nucl. Fusion 30, 1475 (1990).

[24]    J.R.Ferron, M.L.Walker, L.L.Lao, H.E.St.John, D.A.Humphreys and J.A.Leuer, Nucl. Fusion 38, 1055 (1997).

[25]    L.LoDestro and L.D.Pearlstein, Phys. Plasmas 1, 90 (1994).

[26]    J.F.Clark and D.J.Sigmar, Phys. Rev. Lett. 38, 70 (1977).

[27]    S.P.Hirshman and S.C.Jardin, Phys. Fluids 22, 731, (1979).

[28]    V.D.Khait, Sov. J. Plasma Phys. 6, 476 (1980).

[29]    L.L.Lao, S.P.Hirshman and R.M.Wieland, Phys. Fluids 24, 1431 (1981).

[30]    L.L.Lao, Comput. Phys. Commun. 31, 201 (1984).

[31]    P.N.Vabishevich, L.M.Degtyarev and A.P.Favorski, Sov. J. Plasma Phys. 4, 554 (1978).

[32]    L.M.Degtyarev and V.V.Drozdov, Comput. Phys. Rep. 2, 341 (1985).

[33]    L.M.Degtyarev and V.V.Drozdov, Int. J. Mod. Phys. C 2, 30 (1991).

[34]    S.Yu.Medvedev, L.Villard, L.M.Degtyarev, A.A.Martynov, R.Gruber and F.Troyon, 20^th EPS Conf. on Controlled Fusion and Plasma Phys., Lisbon, Proc. Contrib. Papers, Vol. 17C, Part IV, 1279 (1993).

[35]    V.V.Drozdov, The POLAR-2D-NW code, Private coomunications.

[36]    R.R.Khayrutdinov, J.B.Lister, V.E.Lukash and J.P.Wainwright, Plasma Phys. Control. Fusion, 43, 321 (2001).

[37]    J.-Y.Favez, R.R.Khayrutdinov, J.B.Lister, and V.E.Lukash, Plasma Phys. Control. Fusion, 44, 171 (2002).

[38]    ITER Final Design Report. Control System Design and Assessment, Section 5.1.1 "Plasma Current, Position and Shape Control", G45 FDR 1 10-70-13 R1.0, 2001.

[39]    L.E.Zakharov and A.Pletzer, Phys. of Plasmas 6, No. 12, 4693 (1999).

[40]    Ref. on the EFIT code

[41]    IAEA-CN-94/CT/P-08, ITER Preprints, INTERNATIONAL ATOMIC ENERGY AGENCY NINETEENTH FUSION ENERGY CONFERENCE, Lyon, France, 14-19 October 2002.

[42]    Samarsky A.A., Tishkin V.F., Favorsky A.P., et al. \\ Sov. J. Differential Equations  7, 1171 (1981).

[43]   L.M.Degtyarev and V.V.Drozdov, S.Yu.Medvedev. Equilibrium and stability  numerical modeling of toroidal plasma, USSR Academy Science Book, M.V.Keldysh Institute of Applied Mathematics, Moscow 1989.